text:
using JSON2, JSON
using BenchmarkTools

const suite2 = BenchmarkGroup()

suite2["print"] = BenchmarkGroup(["serialize"])
# suite["pretty-print"] = BenchmarkGroup(["serialize"])

struct CustomListType
    x::Int
    y::Float64
    z::Union{CustomListType, Void}
end

struct CustomTreeType
    x::String
    y::Union{CustomTreeType, Void}
    z::Union{CustomTreeType, Void}
end

list(x) = x == 0 ? nothing : CustomListType(1, 1.0, list(x - 1))
tree(x) = x == 0 ? nothing : CustomTreeType("!!!", tree(x - 1), tree(x - 1))

const micros = Dict(
    "integer" => 88,
    "float" => -88.8,
    "ascii" => "Hello World!",
    "ascii-1024" => "x" ^ 1024,
    "unicode" => "ສະບາຍດີຊາວໂລກ!",
    "unicode-1024" => "ℜ" ^ 1024,
    "bool" => true,
    "null" => nothing,
    "flat-homogenous-array-16" => collect(1:16),
    "flat-homogenous-array-1024" => collect(1:1024),
    "heterogenous-array" => [
        1, 2, 3, 7, "A", "C", "E", "N", "Q", "R", "Shuttle to Grand Central"],
    "nested-array-16^2" => [collect(1:16) for _ in 1:16],
    "nested-array-16^3" => [[collect(1:16) for _ in 1:16] for _ in 1:16],
    "small-dict" => Dict(
        :a => :b, :c => "💙💙💙💙💙💙", :e => 10, :f => Dict(:a => :b)),
    "flat-dict-128" => Dict(zip(collect(1:128), collect(1:128))),
    "date" => Date(2016, 08, 09),
    "matrix-16" => eye(16),
    "custom-list-128" => list(128),
    "custom-tree-8" => tree(8))

@benchmark JSON2.write($(IOBuffer()), $(micros["integer"]))
@benchmark JSON.print($(IOBuffer()), $(micros["integer"]))

@benchmark JSON2.write($(IOBuffer()), $(micros["float"]))
@benchmark JSON.print($(IOBuffer()), $(micros["float"]))

@benchmark JSON2.write($(IOBuffer()), $(micros["ascii"]))
@benchmark JSON.print($(IOBuffer()), $(micros["ascii"]))

@benchmark JSON2.write($(IOBuffer()), $(micros["ascii-1024"]))
@benchmark JSON.print($(IOBuffer()), $(micros["ascii-1024"]))

@benchmark JSON2.write($(IOBuffer()), $(micros["bool"]))
@benchmark JSON.print($(IOBuffer()), $(micros["bool"]))

@benchmark JSON2.write($(IOBuffer()), $(micros["null"]))
@benchmark JSON.print($(IOBuffer()), $(micros["null"]))

@benchmark JSON2.write($(IOBuffer()), $(micros["flat-homogenous-array-16"])) # should be better
@benchmark JSON.print($(IOBuffer()), $(micros["flat-homogenous-array-16"]))

@benchmark JSON2.write($(IOBuffer()), $(micros["nested-array-16^2"]))
@benchmark JSON.print($(IOBuffer()), $(micros["nested-array-16^2"]))

@benchmark JSON2.write($(IOBuffer()), $(micros["small-dict"]))  #TODO
@benchmark JSON.print($(IOBuffer()), $(micros["small-dict"]))

@benchmark JSON2.write($(IOBuffer()), $(micros["date"]))
@benchmark JSON.print($(IOBuffer()), $(micros["date"]))

@benchmark JSON2.write($(IOBuffer()), $(micros["custom-list-128"]))
@benchmark JSON.print($(IOBuffer()), $(micros["custom-list-128"]))

@benchmark JSON2.write($(IOBuffer()), $(micros["custom-tree-8"]))
@benchmark JSON.print($(IOBuffer()), $(micros["custom-tree-8"]))

b = @benchmarkable JSON2.read("\"id\"", Symbol); tune!(b); run(b)
b = @benchmarkable JSON2.read("\"id\"", String); tune!(b); run(b)

b = @benchmarkable JSON2.read("{\"id\": 1, \"name\": \"sally\"}"); tune!(b); run(b)
b = @benchmarkable JSON.parse("{\"id\": 1, \"name\": \"sally\"}"); tune!(b); run(b)

@time JSON2.read(IOBuffer(Mmap.mmap("/Users/jacobquinn/Downloads/100k_twitter_users.json")));
@time JSON.parsefile("/Users/jacobquinn/Downloads/100k_twitter_users.json");


text:
using Iterators, Base.Test
using Compat

# count
# -----

i = 0
c0 = countfrom(0, 2)

@test eltype(c0) == Int

for j in c0
	@test j == i*2
	i += 1
	i <= 10 || break
end

# take
# ----

t = take(0:2:8, 10)

@test length(collect(t)) == 5
@test eltype(t) == Int

i = 0
for j = t
	@test j == i*2
	i += 1
end
@test i == 5

i = 0
for j = take(0:2:100, 10)
	@test j == i*2
	i += 1
end
@test i == 10

# drop
# ----

d0 = drop(0:2:10, 2)

@test eltype(d0) == Int

i = 0
for j in d0
	@test j == (i+2)*2
	i += 1
end
@test i == 4

# cycle
# -----

cy1 = cycle(0:3)

@test eltype(cy1) == Int

i = 0
for j in cy1
	@test j == i % 4
	i += 1
	i <= 10 || break
end

# repeated
# --------

r0 = repeated(1, 10)

@test eltype(r0) == Int

i = 0
for j in r0
	@test j == 1
	i += 1
end
@test i == 10

r1 = repeated(1)

@test eltype(r1) == Int

i = 0
for j in r1
	@test j == 1
	i += 1
	i <= 10 || break
end

# repeatedly
# ----------

i = 0
for j = repeatedly(() -> 1, 10)
	@test j == 1
	i += 1
end
@test i == 10
for j = repeatedly(() -> 1)
	@test j == 1
	i += 1
	i <= 10 || break
end

# chain
# -----

ch1 = chain(1:2:5, 0.2:0.1:1.6)

@test eltype(ch1) == typejoin(Int, Float64)
@test collect(ch1) == [1:2:5; 0.2:0.1:1.6]

ch2 = chain(1:0, 1:2:5, 0.2:0.1:1.6)

@test eltype(ch2) == typejoin(Int, Float64)
@test collect(ch2) == [1:2:5; 0.2:0.1:1.6]

# product
# -------

x1 = 1:2:10
x2 = 1:5

p0 = product(x1, x2)

@test eltype(p0) == Tuple{Int, Int}
@test collect(p0) == vec([(y1, y2) for y1 in x1, y2 in x2])

p1 = product()

@test eltype(p1) == Tuple{}
@test length(p1) == 1
@test collect(p1) == [()]

# distinct
# --------

x = [5, 2, 2, 1, 2, 1, 1, 2, 4, 2]
di0 = distinct(x)
@test eltype(di0) == Int
@test collect(di0) == unique(x)

# partition
# ---------

pa0 = partition(take(countfrom(1), 6), 2)
@test eltype(pa0) == Tuple{Int, Int}
@test collect(pa0) == [(1,2), (3,4), (5,6)]

pa1 = partition(take(countfrom(1), 4), 2, 1)
@test eltype(pa1) == Tuple{Int, Int}
@test collect(pa1) == [(1,2), (2,3), (3,4)]

pa2 = partition(take(countfrom(1), 8), 2, 3)
@test eltype(pa2) == Tuple{Int, Int}
@test collect(pa2) == [(1,2), (4,5), (7,8)]

pa3 = partition(take(countfrom(1), 0), 2, 1)
@test eltype(pa3) == Tuple{Int, Int}
@test collect(pa3) == []

@test_throws ArgumentError partition(take(countfrom(1), 8), 2, 0)

# imap
# ----

function test_imap(expected, input...)
    result = collect(imap(+, input...))
    @test result == expected
end

# Empty arrays
test_imap(
    Any[],
    []
)

test_imap(
    Any[],
    Union{}[]
)

# Simple operation
test_imap(
    Any[1,2,3],
    [1,2,3]
)

# Multiple arguments
test_imap(
    Any[5,7,9],
    [1,2,3],
    [4,5,6]
)

# Different-length arguments
test_imap(
    Any[2,4,6],
    [1,2,3],
    countfrom(1)
)

# groupby
# -------

function test_groupby(input, expected)
    result = collect(groupby(x -> x[1], input))
    @test result == expected
end

# Empty arrays
test_groupby(
    [],
    Any[]
)

test_groupby(
    Union{}[],
    Any[]
)

# Singletons
test_groupby(
    ["xxx"],
    Any[["xxx"]]
)

# Typical operation
test_groupby(
    ["face", "foo", "bar", "book", "baz"],
    Any[["face", "foo"], ["bar", "book", "baz"]]
)

# Trailing singletons
test_groupby(
    ["face", "foo", "bar", "book", "baz", "xxx"],
    Any[["face", "foo"], ["bar", "book", "baz"], ["xxx"]]
)

# Leading singletons
test_groupby(
    ["xxx", "face", "foo", "bar", "book", "baz"],
    Any[["xxx"], ["face", "foo"], ["bar", "book", "baz"]]
)

# Middle singletons
test_groupby(
    ["face", "foo", "xxx", "bar", "book", "baz"],
    Any[["face", "foo"], ["xxx"], ["bar", "book", "baz"]]
)

# subsets
# -------

s0 = subsets(Any[])
@test eltype(eltype(s0)) == Any
@test collect(s0) == Vector{Any}[Any[]]

s1 = subsets([:a])
@test eltype(eltype(s1)) == Symbol
@test collect(s1) == Vector{Symbol}[Symbol[], Symbol[:a]]

s2 = subsets([:a, :b, :c])
@test eltype(eltype(s2)) == Symbol
@test collect(s2) == Vector{Symbol}[
    Symbol[], Symbol[:a], Symbol[:b], Symbol[:a, :b], Symbol[:c],
    Symbol[:a, :c], Symbol[:b, :c], Symbol[:a, :b, :c],
]

# subsets of size k
# -----------------

sk0 = subsets(Any[],0)
@test eltype(eltype(sk0)) == Any
@test collect(sk0) == Vector{Any}[Any[]]

sk1 = subsets([:a, :b, :c], 1)
@test eltype(eltype(sk1)) == Symbol
@test collect(sk1) == Vector{Symbol}[Symbol[:a], Symbol[:b], Symbol[:c]]

sk2 = subsets([:a, :b, :c], 2)
@test eltype(eltype(sk2)) == Symbol
@test collect(sk2) == Vector{Symbol}[Symbol[:a,:b], Symbol[:a,:c], Symbol[:b,:c]]

sk3 = subsets([:a, :b, :c], 3)
@test eltype(eltype(sk3)) == Symbol
@test collect(sk3) == Vector{Symbol}[Symbol[:a,:b,:c]]

sk4 = subsets([:a, :b, :c], 4)
@test eltype(eltype(sk4)) == Symbol
@test collect(sk4) == Vector{Symbol}[]

sk5 = subsets(collect(1:4), 1)
@test eltype(eltype(sk5)) == Int
@test length(collect(sk5)) == binomial(4,1)

sk6 = subsets(collect(1:4), 2)
@test eltype(eltype(sk6)) == Int
@test length(collect(sk6)) == binomial(4,2)

sk7 = subsets(collect(1:4), 3)
@test eltype(eltype(sk7)) == Int
@test length(collect(sk7)) == binomial(4,3)

sk8 = subsets(collect(1:4), 4)
@test eltype(eltype(sk8)) == Int
@test length(collect(sk8)) == binomial(4,4)

sk9 = subsets(collect(1:4), 5)
@test eltype(eltype(sk9)) == Int
@test length(collect(sk9)) == binomial(4,5)

# nth element
# -----------
for xs in Any[[1, 2, 3], 1:3, reshape(1:3, 3, 1)]
    @test nth(xs, 3) == 3
    @test_throws BoundsError nth(xs, 0)
    @test_throws BoundsError nth(xs, 4)
end

for xs in Any[take(1:3, 3), drop(-1:3, 2)]
    @test nth(xs, 3) == 3
    @test_throws BoundsError nth(xs, 0)
end

s = subsets([1, 2, 3])
@test_throws BoundsError nth(s, 0)
@test_throws BoundsError nth(s, length(s) + 1)

# Every nth
# ---------

tn0 = takenth(Any[], 10)
@test eltype(tn0) == Any
@test collect(tn0) == Any[]

tn1 = takenth(Int[], 10)
@test eltype(tn1) == Int
@test collect(tn1) == Int[]

@test_throws ArgumentError takenth([], 0)

tn2 = takenth(10:20, 3)
@test eltype(tn2) == Int
@test collect(tn2) == [12,15,18]

tn3 = takenth(10:20, 1)
@test eltype(tn3) == Int
@test collect(tn3) == collect(10:20)

# peekiter
# --------
result = Int[]
for x in peekiter(1:10) push!(result, x) end
@test result == collect(1:10)

result = Int[]
for x in peekiter([]) push!(result, x) end
@test result == []

result = Int[]
for x in peekiter(1:10) push!(result, x) end
@test result == collect(1:10)

it = peekiter([:a, :b, :c])
s = start(it)
@test get(peek(it, s)) == :a

it = peekiter([])
s = start(it)
@test isnull(peek(it, s))

it = peekiter(1:10)
s = start(it)
x, s = next(it, s)
@test get(peek(it, s)) == 2

## @itr
## ====

# @zip
# ----

macro test_zip(input...)
    n = length(input)
    x = Expr(:tuple, ntuple(i->gensym(), n)...)
    v = Expr(:tuple, map(esc, input)...)
    w = :(zip($(map(esc, input)...)))
    quote
	br = Any[]
	for $x in zip($v...)
	    push!(br, $x)
	end
	mr = Any[]
	@itr for $x in $w
	    push!(mr, $x)
	end
	@test br == mr
    end
end

@test_zip [1,2,3] [:a, :b, :c] ['x', 'y', 'z']
@test_zip [1,2,3] [:a, :b] ['w', 'x', 'y', 'z']
@test_zip [1,2,3] Any[] ['w', 'x', 'y', 'z']
@test_zip [1,2,3] Union{}[] ['w', 'x', 'y', 'z']

# @enumerate
# ----------

macro test_enumerate(input)
    i = gensym()
    x = gensym()
    v = esc(input)
    quote
	br = Any[]
	for ($i,$x) in enumerate($v)
	    push!(br, ($i,$x))
	end
	mr = Any[]
	@itr for ($i,$x) in enumerate($v)
	    push!(mr, ($i,$x))
	end
	@test br == mr
    end
end

@test_enumerate [:a, :b, :c]
@test_enumerate Union{}[]
@test_enumerate Any[]

# @take
# -----

macro test_take(input, n)
    x = gensym()
    v = esc(input)
    quote
	br = Any[]
	for $x in take($v, $n)
	    push!(br, $x)
	end
	mr = Any[]
	@itr for $x in take($v, $n)
	    push!(mr, $x)
	end
	@test br == mr
    end
end

@test_take [:a, :b, :c] 2
@test_take [:a, :b, :c] 5
@test_take [:a, :b, :c] 0
@test_take Any[] 2
@test_take Union{}[] 2
@test_take Any[] 0
@test_take [(:a,1), (:b,2), (:c,3)] 2

# @takestrict
# -----

macro test_takestrict(input, n)
    x = gensym()
    v = esc(input)
    quote
	br = Any[]
	bfailed = false
	try
	    for $x in takestrict($v, $n)
		push!(br, $x)
	    end
	catch
	    bfailed = true
	end

	mr = Any[]
	mfailed = false
	try
	    @itr for $x in takestrict($v, $n)
		push!(mr, $x)
	    end
	catch
	    mfailed = true
	end
	@test br == mr
	@test bfailed == mfailed
    end
end

@test_takestrict [:a, :b, :c] 2
@test_takestrict [:a, :b, :c] 3
@test_takestrict [:a, :b, :c] 5
@test_takestrict [:a, :b, :c] 0
@test_takestrict Any[] 2
@test_takestrict Union{}[] 2
@test_takestrict Any[] 0
@test_takestrict [(:a,1), (:b,2), (:c,3)] 2
@test_takestrict [(:a,1), (:b,2), (:c,3)] 3
@test_takestrict [(:a,1), (:b,2), (:c,3)] 4

# @drop
# -----

macro test_drop(input, n)
    x = gensym()
    v = esc(input)
    quote
	br = Any[]
	for $x in drop($v, $n)
	    push!(br, $x)
	end
	mr = Any[]
	@itr for $x in drop($v, $n)
	    push!(mr, $x)
	end
	@test br == mr
    end
end

@test_drop [:a, :b, :c] 2
@test_drop [:a, :b, :c] 5
@test_drop [:a, :b, :c] 0
@test_drop Any[] 2
@test_drop Union{}[] 2
@test_drop Any[] 0
@test_drop [(:a,1), (:b,2), (:c,3)] 2

# @chain
# -----

macro test_chain(input...)
    x = gensym()
    v = Expr(:tuple, map(esc, input)...)
    w = :(chain($(map(esc, input)...)))
    quote
	br = Any[]
	for $x in chain($v...)
	    push!(br, $x)
	end
	mr = Any[]
	@itr for $x in $w
	    push!(mr, $x)
	end
	@test br == mr
    end
end

@test_chain [1,2,3] [:a, :b, :c] ['x', 'y', 'z']
@test_chain [1,2,3] [:a, :b] ['w', 'x', 'y', 'z']
@test_chain [1,2,3] Any[] ['w', 'x', 'y', 'z']
@test_chain [1,2,3] Union{}[] ['w', 'x', 'y', 'z']
@test_chain [1,2,3] 4 [('w',3), ('x',2), ('y',1), ('z',0)]

ext: jl


text:
# test for comparisons and hashing

let
    x = HTMLText("test")
    y = HTMLText("test")
    x1 = HTMLText("test1")
    @test x == y
    @test hash(x) == hash(y)
    @test x1 != y
    @test hash(x1) != hash(y)
end

let
    x = HTMLElement(:div)
    y = HTMLElement(:div)
    @test x == y
    @test hash(x) == hash(y)
    push!(x, HTMLElement(:p))
    @test x != y
    @test hash(x) != hash(y)
    push!(y, HTMLElement(:p))
    @test x == y
    @test hash(x) == hash(y)
    setattr!(x,"class","test")
    @test x != y
    @test hash(x) != hash(y)
    setattr!(y,"class","test")
    @test x == y
    @test hash(x) == hash(y)
end

let
    x = HTMLDocument("html", HTMLElement(:html))
    y = HTMLDocument("html", HTMLElement(:html))
    @test x == y
    @test hash(x) == hash(y)
    x.doctype = ""
    @test x != y
    @test hash(x) != hash(y)
    y.doctype = ""
    @test x == y
    @test hash(x) == hash(y)
end


text:
function calc_shapes(ncell::Int64, sol::Array{Float64})

    E = zeros(ncell+2); B = zeros(ncell+2); F = zeros(ncell+2); S = zeros(ncell+2)
    dfde = zeros(ncell+2); del = zeros(ncell+2)

    for i = 1:ncell+2
        del[i] = sol[1,i]
        E[i] = sol[2,i]/del[i] - 1.
        F[i] = 4.8274*E[i]^4 - 5.9816*E[i]^3 + 4.0274*E[i]^2 + 0.23247*E[i] + 0.15174
        if E[i] < -0.0616
            B[i] = -225.86*E[i]^3 - 3016.6*E[i]^2 - 208.68*E[i] - 17.915
        elseif E[i] > -0.0395
            B[i] = 131.9*E[i]^3 - 167.32*E[i]^2 + 76.642*E[i] - 11.068
        else
            B[i] = 0.5*(-225.86*E[i]^3 - 3016.6*E[i]^2 - 208.68*E[i] - 17.915
                        + 131.9*E[i]^3 - 167.32*E[i]^2 + 76.642*E[i] - 11.068)
        end
        if E[i] < -0.0582
            S[i] = 451.55*E[i]^3 + 2010*E[i]^2 + 138.96*E[i] + 11.296
        elseif E[i] > -0.042
            S[i] = -96.739*E[i]^3 + 117.74*E[i]^2 - 46.432*E[i] + 6.8074
        else
            S[i] = 0.5*(451.55*E[i]^3 + 2010*E[i]^2 + 138.96*E[i] + 11.296
                        - 96.739*E[i]^3 + 117.74*E[i]^2 - 46.432*E[i] + 6.8074)
        end
        dfde[i] = 4*4.8274*E[i]^3 - 3*5.9816*E[i]^2 + 2*4.0274*E[i] + 0.23247
    end

    return E, F, B, S, dfde, del
end

function calc_eigen(ncell::Int64, E::Array{Float64}, F::Array{Float64},
                    dfde::Array{Float64}, ue::Array{Float64})

    lamb1 = zeros(ncell+2); lamb2 = zeros(ncell+2)

    for i = 1:ncell+2
        a_q = 1.
        b_q = -ue[i] * (dfde[i] - 1)
        c_q = ue[i] * ue[i] * (E[i]*dfde[i] - F[i])
        lamb1[i] = (-b_q + sqrt(b_q*b_q - 4*a_q*c_q))/(2*a_q)
        lamb2[i] = (-b_q  -sqrt(b_q*b_q - 4*a_q*c_q))/(2*a_q)

        #Always have lamb1 > lamb2
        if lamb2[i]>lamb1[i]
            temp  = lamb2[i]
            lamb2[i] = lamb1[i]
            lamb1[i] = temp
        end
    end

    return lamb1, lamb2
end

function calc_dt(ncell::Int64, cfl::Float64, dx::Float64, lamb1::Array{Float64},
                 lamb2::Array{Float64})

    dt = 10000
    for i = 2:ncell+1
        dti = cfl*(dx/(abs(lamb1[i]) + abs(lamb2[i])))
        if dti < dt
            dt = dti
        end
    end

    return dt
end

function calc_fluxes(ncell::Int64, E::Array{Float64}, F::Array{Float64},
                     del::Array{Float64}, ue::Array{Float64}, lamb1::Array{Float64}, lamb2::Array{Float64})

    Apos = zeros(2,2); Aneg = zeros(2,2)
    flux = zeros(2,2,ncell+2)

    for i = 1:ncell+2
        if lamb1[i] >= 0. && lamb2[i] >= 0.
            flux[1,1,i] = ue[i]*E[i]*del[i]
            flux[1,2,i] = ue[i]*F[i]*del[i]
        elseif lamb1[i] <= 0. && lamb2[i] <= 0.
            flux[1,:,i] = 0.
        else
            Apos[1,1] = ue[i]*lamb1[i]/(lamb1[i] - lamb2[i])*(-1. - lamb2[i]/ue[i])
            Apos[1,2] = ue[i]*lamb1[i]/(lamb1[i] - lamb2[i])
            Apos[2,1] = -ue[i]*lamb1[i]/(lamb1[i] - lamb2[i])*(1. + lamb1[i]/ue[i])*(1. + lamb2[i]/ue[i])
            Apos[2,2] = ue[i]*lamb1[i]/(lamb1[i] - lamb2[i])*(1. + lamb1[i]/ue[i])
            flux[1,1,i] = Apos[1,1]*del[i] + Apos[1,2]*(E[i] + 1.)*del[i]
            flux[1,2,i] = Apos[2,1]*del[i] + Apos[2,2]*(E[i] + 1.)*del[i]
        end
    end

    for i = 1:ncell+2
        if lamb1[i] >= 0. && lamb2[i] >= 0.
            flux[2,1,i] = 0
            flux[2,2,i] = 0
        elseif lamb1[i] < 0. && lamb2[i] < 0.
            flux[2,1,i]= ue[i]*E[i]*del[i]
            flux[2,2,i]= ue[i]*F[i]*del[i]
        else
            Aneg[1,1] = ue[i]*lamb2[i]/(lamb1[i]-lamb2[i])*(1. + lamb1[i]/ue[i])
            Aneg[1,2] = -ue[i]*lamb2[i]/(lamb1[i]-lamb2[i])
            Aneg[2,1] = ue[i]*lamb2[i]/(lamb1[i]-lamb2[i])*(1. + lamb1[i]/ue[i])*(1. + lamb2[i]/ue[i])
            Aneg[2,2] = ue[i]*lamb2[i]/(lamb1[i]-lamb2[i])*(-1. - lamb2[i]/ue[i])
            flux[2,1,i] = Aneg[1,1]*del[i] + Aneg[1,2]*(E[i]+1)*del[i]
            flux[2,2,i] = Aneg[2,1]*del[i] + Aneg[2,2]*(E[i]+1)*del[i]
        end
    end

    return flux
end

function matuCyl(tTot::Float64, ncell::Int64 = 200, cfl::Float64 = 0.8, nstage::Int64 = 2)

    dx = pi/ncell

    #Timestep determined dynamically
    dt = 0

    uinf = 1.; c = 1.; re = 1000.

    x = zeros(ncell+2); ue = zeros(ncell+2); uex = zeros(ncell+2); uet = zeros(ncell+2)
    E = zeros(ncell+2); B = zeros(ncell+2); F = zeros(ncell+2); S = zeros(ncell+2); dfde = zeros(ncell+2);
    del = zeros(ncell+2);

    sol = zeros(2,ncell+2); solt = zeros(2,ncell+2)
    lamb1 = zeros(ncell+2); lamb2 = zeros(ncell+2)
    Jsep = zeros(2,ncell+2)
    Csep = zeros(ncell+2)

    flux = zeros(2,2,ncell+2)
    cf = zeros(ncell+2)

    #Set sources (cell centred)
    for i = 2:ncell+1
        x[i] = (i-1.5)*dx
        ue[i] = 2*sin(x[i])
        uex[i] = 2*cos(x[i])
    end
    x[1] = 2*x[2] - x[3]
    x[ncell+2] = 2*x[ncell+1] - x[ncell]
    uet[:] = zeros(ncell+2)
    ue[1] = 0.
    ue[ncell+2] = 0.
    uex[1] = cos(0.)
    uex[ncell+2] = cos(pi)

    #Set initial conditions
    E_init = 0.4142
    for i = 1:ncell+2
        E[i] = E_init
        B[i] = 131.9*E_init^3 - 167.32*E_init^2 + 76.642*E_init - 11.068
        del[i] = sqrt.(B[i]*0.005)
        F[i] = 4.8274*E_init^4 - 5.9816*E_init^3 + 4.0274*E_init^2 + 0.23247*E_init + 0.15174
    end

    t = 0.005

    sol[1,:] = del[:]
    sol[2,:] = del[:].*(E[:] .+ 1)

    rhs = zeros(2, ncell+2)

    istep = 0
    while t < tTot
        istep += 1

        #Update ghost cells
        sol[:,1] = 2*sol[:,2] - sol[:,3]
        sol[:,ncell+2] = 2*sol[:,ncell+1] - sol[:,ncell]

        #Calculate derived quantities
        E, F, B, S, dfde, del = calc_shapes(ncell, sol)

        #Compute eigenvalues
        lamb1, lamb2 = calc_eigen(ncell, E, F, dfde, ue)

        #Compute timestep
        dt = calc_dt(ncell, cfl, dx, lamb1, lamb2)

        flux = calc_fluxes(ncell, E, F, del, ue, lamb1, lamb2)

        rhs[:,:] = zeros(2, ncell+2)

        #The ghost rhs are never used - dont matter
        for i = 2:ncell+1
            rhs[1,i] = B[i]/(2*del[i]) - del[i]*uet[i]/ue[i] - (E[i]+1)*del[i]*uex[i]
            rhs[2,i] = S[i]/del[i] - 2*E[i]*del[i]*uet[i]/ue[i] - 2*F[i]*del[i]*uex[i]
        end

        for i = 2:ncell
            solt[1,i] = sol[1,i] - (dt/dx)*(flux[1,1,i] - flux[1,1,i-1] + flux[2,1,i+1]
                                            -flux[2,1,i]) + dt*rhs[1,i]
            solt[2,i] = sol[2,i] - (dt/dx)*(flux[1,2,i] - flux[1,2,i-1] + flux[2,2,i+1]
                                            -flux[2,2,i]) + dt*rhs[2,i]
        end

        i = ncell+1
        solt[1,i] = sol[1,i] - (dt/dx)*(flux[1,1,i] - flux[1,1,i-1]) + dt*rhs[1,i]
        solt[2,i] = sol[2,i] - (dt/dx)*(flux[1,2,i] - flux[1,2,i-1]) + dt*rhs[2,i]

        #Update ghost cells
        solt[:,1] = 2*solt[:,2] - solt[:,3]
        solt[:,ncell+2] = 2*solt[:,ncell+1] - solt[:,ncell]

        #Calculate derived quantities
        E, F, B, S, dfde, del = calc_shapes(ncell, solt)

        #Compute eigenvalues
        lamb1, lamb2 = calc_eigen(ncell, E, F, dfde, ue)

        fluxt = calc_fluxes(ncell, E, F, del, ue, lamb1, lamb2)

        rhs[:,:] = zeros(2, ncell+2)
        #The ghost rhs are never used - dont matter
        for i = 2:ncell+1
            rhs[1,i] = B[i]/(2*del[i]) - del[i]*uet[i]/ue[i] - (E[i]+1)*del[i]*uex[i]
            rhs[2,i] = S[i]/del[i] - 2*E[i]*del[i]*uet[i]/ue[i] - 2*F[i]*del[i]*uex[i]
        end

        i = 2
        sol[1,i] = 0.5*(sol[1,i] + solt[1,i]) - (0.5*dt/dx)*(flux[1,1,i] - flux[1,1,i-1]
                                                             - flux[2,1,i] + 2*flux[2,1,i+1] - flux[2,1,i+2] + fluxt[1,1,i] - fluxt[1,1,i-1]
                                                             + fluxt[2,1,i+1] - fluxt[2,1,i]) + 0.5*dt*rhs[1,i]
        sol[2,i] = 0.5*(sol[2,i] + solt[2,i]) - (0.5*dt/dx)*(flux[1,2,i] - flux[1,2,i-1]
                                                             - flux[2,2,i] + 2*flux[2,2,i+1]-flux[2,2,i+2]+fluxt[1,2,i] - fluxt[1,2,i-1]
                                                             + fluxt[2,2,i+1] - fluxt[2,2,i])+0.5*dt*rhs[2,i]
        for i = 3:ncell
            sol[1,i] = 0.5*(sol[1,i] + solt[1,i]) - (0.5*dt/dx)*(flux[1,1,i] - 2*flux[1,1,i-1] + flux[1,1,i-2]
                                                                 - flux[2,1,i] + 2*flux[2,1,i+1] - flux[2,1,i+2] + fluxt[1,1,i] - fluxt[1,1,i-1]
                                                                 + fluxt[2,1,i+1] - fluxt[2,1,i]) + 0.5*dt*rhs[1,i]
            sol[2,i] = 0.5*(sol[2,i] + solt[2,i]) - (0.5*dt/dx)*(flux[1,2,i] - 2*flux[1,2,i-1] + flux[1,2,i-2]
                                                                 - flux[2,2,i] + 2*flux[2,2,i+1] - flux[2,2,i+2] + fluxt[1,2,i] - fluxt[1,2,i-1]
                                                                 + fluxt[2,2,i+1] - fluxt[2,2,i]) + 0.5*dt*rhs[2,i]
        end
        i = ncell+1
        sol[1,i] = 0.5*(sol[1,i] + solt[1,i]) - (0.5*dt/dx)*(flux[1,1,i] - 2*flux[1,1,i-1] + flux[1,1,i-2]
                                                             +fluxt[1,1,i] - fluxt[1,1,i-1]) + 0.5*dt*rhs[1,i]
        sol[2,i] = 0.5*(sol[2,i] + solt[2,i]) - (0.5*dt/dx)*(flux[1,2,i] - 2*flux[1,2,i-1] + flux[1,2,i-2]
                                                             +fluxt[1,2,i] - fluxt[1,2,i-1]) + 0.5*dt*rhs[2,i]

        t = t + dt
        if t > tTot
            break
        end

        for i = 2:ncell+1
            cf[i] = B[i]*ue[i]*c/(del[i]*uinf*re)
            Csep[i] = (del[i+1] - del[i])/(del[i] - del[i-1])
            if Csep[i] > 10.
                println("singularity (separation) detected at t=$t, x=$(x[i])")
            end
        end
        #println("$istep \t", "$dt \t", "$t \t", "$(maximum(Csep)) \t", "$(minimum(Csep)) \n")
        println("$t \t")
        display(plot(x/pi,[del, E] , layout=(2,1), legend = false))
        sleep(0.1)

    end

mat = zeros(0, 8)
mat = mat'
for i = 2:ncell+1
    mat = hcat(mat,[x[i], ue[i], uex[i], uet[i], del[i], E[i], cf[i], Csep[i]])
end
mat = mat'

return mat
end

# f = open("t$istep", "w")
# write(f, ["#x \t", "#ue \t", "#uex \t", "#uet \t", "#del \t", "#E \t", "#cf \t", "#Csep \n"])
# writedlm(f, mat)
# close(f)

max_stars_repo_name: desanga/IBMFVM

max_issues_repo_name: desanga/IBMFVM

max_forks_repo_name: desanga/IBMFVM


text:
using DeepQLearning
using POMDPModels
using POMDPSimulators
using POMDPPolicies
using Flux
using Random
using RLInterface
using StaticArrays
using Test
Random.seed!(7)
GLOBAL_RNG = MersenneTwister(1) # for test consistency

include("test_env.jl")

function evaluate(mdp, policy, rng, n_ep=100, max_steps=100)
    avg_r = 0.
    sim = RolloutSimulator(rng=rng, max_steps=max_steps)
    for i=1:n_ep
        DeepQLearning.resetstate!(policy)
        avg_r += simulate(sim, mdp, policy)
    end
    return avg_r/=n_ep
end

@testset "vanilla DQN" begin
    mdp = TestMDP((5,5), 4, 6)
    model = Chain(x->flattenbatch(x), Dense(100, 8, tanh), Dense(8, length(actions(mdp))))
    max_steps = 10000
    exploration = EpsGreedyPolicy(mdp, LinearDecaySchedule(start=1.0, stop=0.01, steps=max_steps/2),
                                  rng=GLOBAL_RNG)
    solver = DeepQLearningSolver(qnetwork = model, max_steps=max_steps, learning_rate=0.005,
                                 exploration_policy = exploration,
                                 eval_freq=2000,num_ep_eval=100,
                                 log_freq = 500,
                                 double_q = false, dueling=false, prioritized_replay=false)

    policy = solve(solver, mdp)
    r_basic = evaluate(mdp, policy, GLOBAL_RNG)
    @test r_basic >= 1.5
    @test size(actionvalues(policy, rand(GLOBAL_RNG, initialstate(mdp)))) == (length(actions(mdp)),)
end

@testset "double Q DQN" begin
    mdp = TestMDP((5,5), 4, 6)
    model = Chain(x->flattenbatch(x), Dense(100, 8, tanh), Dense(8, length(actions(mdp))))
    max_steps = 10000
    exploration = EpsGreedyPolicy(mdp, LinearDecaySchedule(start=1.0, stop=0.01, steps=max_steps/2),
                                  rng=GLOBAL_RNG)
    solver = DeepQLearningSolver(qnetwork=model,max_steps=max_steps, learning_rate=0.005, eval_freq=2000,num_ep_eval=100,
                                 exploration_policy = exploration,
                                 log_freq = 500,
                                 double_q = true, dueling=false, prioritized_replay=false)

    policy = solve(solver, mdp)
    r_double =  evaluate(mdp, policy, GLOBAL_RNG)
    @test r_double >= 1.5
end

@testset "dueling DQN" begin
    mdp = TestMDP((5,5), 4, 6)
    model = Chain(x->flattenbatch(x), Dense(100, 8, tanh), Dense(8, length(actions(mdp))))
    max_steps = 10000
    exploration = EpsGreedyPolicy(mdp, LinearDecaySchedule(start=1.0, stop=0.01, steps=max_steps/2),
                                  rng=GLOBAL_RNG)
    solver = DeepQLearningSolver(qnetwork = model, max_steps=max_steps, learning_rate=0.005,
                                 eval_freq=2000,num_ep_eval=100,
                                 exploration_policy = exploration,
                                 log_freq = 500,
                                 double_q = false, dueling=true, prioritized_replay=false)

    policy = solve(solver, mdp)
    r_duel =  evaluate(mdp, policy, GLOBAL_RNG)
    @test r_duel >= 1.5
end

@testset "Prioritized DDQN" begin
    mdp = TestMDP((5,5), 4, 6)
    model = Chain(x->flattenbatch(x), Dense(100, 8, tanh), Dense(8, length(actions(mdp))))
    max_steps = 10000
    exploration = EpsGreedyPolicy(mdp, LinearDecaySchedule(start=1.0, stop=0.01, steps=max_steps/2),
                                  rng=GLOBAL_RNG)
    solver = DeepQLearningSolver(qnetwork = model, max_steps=max_steps, learning_rate=0.005,
                                 exploration_policy = exploration,
                                 eval_freq=2000,num_ep_eval=100,
                                 log_freq = 500,
                                 double_q = true, dueling=true, prioritized_replay=true)

    policy = solve(solver, mdp)
    r_ddqn =  evaluate(mdp, policy, GLOBAL_RNG)
    @test r_ddqn >= 1.5
end

# DRQN tests

@testset "TestMDP DRQN" begin
    mdp = TestMDP((5,5), 1, 6)
    model = Chain(x->flattenbatch(x), LSTM(25, 8), Dense(8, length(actions(mdp))))
    max_steps = 10000
    exploration = EpsGreedyPolicy(mdp, LinearDecaySchedule(start=1.0, stop=0.01, steps=max_steps/2),
                                  rng=GLOBAL_RNG)
    solver = DeepQLearningSolver(qnetwork = model, max_steps=max_steps, learning_rate=0.005,
                                 exploration_policy = exploration,
                                 eval_freq=2000,num_ep_eval=100,
                                 log_freq = 500,
                                 double_q = true, dueling=false, recurrence=true)
    policy = solve(solver, mdp)
    r_drqn =  evaluate(mdp, policy, GLOBAL_RNG)
    @test r_drqn >= 0.
end

@testset "GridWorld DDRQN" begin
    mdp = SimpleGridWorld();
    model = Chain(x->flattenbatch(x), LSTM(2, 32), Dense(32, length(actions(mdp))))
    max_steps = 10000
    exploration = EpsGreedyPolicy(mdp, LinearDecaySchedule(start=1.0, stop=0.01, steps=max_steps/2),
                                  rng=GLOBAL_RNG)
    solver = DeepQLearningSolver(qnetwork = model, prioritized_replay=false, max_steps=max_steps,
                                 exploration_policy = exploration,
                                 learning_rate=0.001,log_freq=500,
                                 recurrence=true,trace_length=10,
                                 double_q=true, dueling=true)

    policy = solve(solver, mdp)
    sim = RolloutSimulator(rng=GLOBAL_RNG, max_steps=10)
    r_drqn =  evaluate(mdp, policy, GLOBAL_RNG)
    @test r_drqn >= 0.
end

@testset "TigerPOMDP DDRQN" begin
    pomdp = TigerPOMDP(0.01, -1.0, 0.1, 0.8, 0.95);
    input_dims = reduce(*, size(convert_o(Vector{Float64}, first(observations(pomdp)), pomdp)))
    model = Chain(x->flattenbatch(x), LSTM(input_dims, 4), Dense(4, length(actions(pomdp))))
    max_steps = 10000
    exploration = EpsGreedyPolicy(pomdp, LinearDecaySchedule(start=1.0, stop=0.01, steps=max_steps/2),
                                  rng=GLOBAL_RNG)
    solver = DeepQLearningSolver(qnetwork = model, prioritized_replay=false, max_steps=max_steps,
                             learning_rate=0.0001, exploration_policy = exploration,
                             log_freq=500, target_update_freq = 1000,
                             recurrence=true,trace_length=10, double_q=true, dueling=true, max_episode_length=100)

    policy = solve(solver, pomdp)
    @test size(actionvalues(policy, true)) == (length(actions(pomdp)),)
end

mutable struct StaticArrayMDP <: MDP{typeof(SVector(1)), Int64}
    state::typeof(SVector(1))
end
POMDPs.discount(::StaticArrayMDP) = 0.95f0
function POMDPs.initialstate(m::StaticArrayMDP)
    ImplicitDistribution() do rng
        m.state
    end
end

function POMDPs.gen(m::StaticArrayMDP, s, a, rng::AbstractRNG)
    return (sp=s + SVector(a), r=m.state[1]^2)
end

POMDPs.isterminal(::StaticArrayMDP, s) = s[1] >= 3
POMDPs.actions(::StaticArrayMDP) = [0,1]

@testset "Static Array Env" begin
    mdp = StaticArrayMDP(SVector(1))

    model = Chain(Dense(1, 32), Dense(32, length(actions(mdp))))

    max_steps = 10000
    exploration = EpsGreedyPolicy(mdp, LinearDecaySchedule(start=1.0, stop=0.01, steps=5),
                                  rng=GLOBAL_RNG)

    solver = DeepQLearningSolver(qnetwork = model, max_steps=10, exploration_policy=exploration,
                                learning_rate=0.005,log_freq=500,
                                recurrence=false,double_q=true, dueling=true, prioritized_replay=true)
    policy = solve(solver, mdp)

    @test evaluate(mdp, policy, GLOBAL_RNG) > 1.0
end


text:
#= Custom calibration configuration file. =#

using Distributions
using StatsBase
using LinearAlgebra
using CalibrateEDMF
using CalibrateEDMF.ModelTypes
using CalibrateEDMF.LESUtils
using CalibrateEDMF.TurbulenceConvectionUtils
# Import prior constraints
using EnsembleKalmanProcesses.ParameterDistributions
using JLD2

# Cases defined as structs for quick access to default configs
struct GridSearchCases end

function get_config()
    config = Dict()
    # Define regularization of inverse problem
    config["regularization"] = get_regularization_config()
    # Define reference used in the inverse problem
    config["reference"] = get_reference_config(GridSearchCases())
    # Define the SCM static configuration
    config["scm"] = get_scm_config()
    # grid search configuration
    config["grid_search"] = get_grid_search_config()
    return config
end

function get_regularization_config()
    config = Dict()
    # Regularization of observations: mean and covariance
    config["perform_PCA"] = true # Performs PCA on data
    config["variance_loss"] = 1.0e-2 # Variance truncation level in PCA
    config["normalize"] = true  # whether to normalize data by pooled variance
    config["tikhonov_mode"] = "relative" # Tikhonov regularization
    config["tikhonov_noise"] = 1.0e-4 # Tikhonov regularization
    config["dim_scaling"] = true # Dimensional scaling of the loss

    # Parameter regularization: L2 regularization with respect to prior mean.
    #  - Set to `nothing` to use prior covariance as regularizer,
    #  - Set to a float for isotropic parameter regularization.
    #  - Pass a dictionary of lists similar to config["prior_mean"] for
    #       anisotropic regularization. The dictionary must be complete.
    #       If you want to avoid regularizing a certain parameter, set the entry
    #       to [0].
    # To turn off regularization, set config["process"]["augmented"] to false.
    config["l2_reg"] = nothing
    return config
end

function get_reference_config(::GridSearchCases)
    config = Dict()
    config["case_name"] = ["TRMM_LBA", "Rico", "TRMM_LBA"]
    n_cases = length(config["case_name"])
    # Flag to indicate source of data (LES or SCM) for reference data and covariance
    config["y_reference_type"] = SCM()
    # config["Σ_reference_type"] = LES()
    # Fields to learn from during training
    config["y_names"] = repeat([["ql_mean", "qt_mean", "total_flux_qt"]], n_cases)
    # LES data can be stored as an Artifact and downloaded lazily
    # config["y_dir"] = [...]  # only needed for loss map computation
    # provide list of dirs if different from `y_dir`
    # config["Σ_dir"] = [...]
    config["t_start"] = repeat([100.0], n_cases)
    config["t_end"] = repeat([200.0], n_cases)
    config["scm_parent_dir"] = repeat(["scm_init"], n_cases)
    config["scm_suffix"] = repeat(["0"], n_cases)
    return config
end

function get_scm_config()
    config = Dict()
    # List of tuples like [("time_stepping", "dt_min", 1.0)], or nothing
    config["namelist_args"] = nothing
    config["namelist_args"] = [
        ("time_stepping", "t_max", 200.0),
        ("turbulence", "EDMF_PrognosticTKE", "entrainment", "moisture_deficit"),
        ("turbulence", "EDMF_PrognosticTKE", "stochastic_entrainment", "prognostic_noisy_relaxation_process"),
        ("time_stepping", "adapt_dt", false),
    ]
    return config
end

function get_grid_search_config()
    config = Dict()
    # Grid search is performed over each pair of parameters, across all specified values
    config["parameters"] = Dict(
        "general_stochastic_ent_params_{1}" => [0.3, 0.4],
        "general_stochastic_ent_params_{2}" => [0.2, 0.3],
        "entrainment_factor" => [0.1],
    )
    # Number of simulations to run with identical configuration (except random seed)
    config["ensemble_size"] = 1
    # grid search output data stored in `<output_root_dir>/output/YYmmdd_abc`
    # config["output_root_dir"] = pwd()  # provided explicitly in `grid_search_test.jl`
    # Perform grid search and loss map calculation for either the `reference` or `validation` set.
    config["sim_type"] = "reference"
    return config
end


text:
@example :DICT_MERGE begin
	merge(["l"=>"l"],["l"=>6])
end

@example :STR_CONCAT begin
	"l"+"l"
end

@example :REGEX_ITER begin
	l = [6:6]
	for i in eachmatch(r"\d","6")
		l[i]
	end
end

@example :REGEX_ITER begin
	l = {"l"=>"l"}
	for i in eachmatch(r"[a-z]", "l")
		l[i]
	end
end

@example :INVARIANTS begin
	function l(::Array{Number,1});end
	l([6])
end

@example :NOT_PYTHON begin
	eval(parse("if true"))
end

@example :NOT_PYTHON begin
	eval(parse("end"))
end

@example :BASE_ANDOR begin
	eval(parse("Base.&&"))
end

@example :BASE_ANDOR begin
	eval(parse("Base.||"))
end

@example :ARRAY_SIZE begin
	function l(::Array{Number,6});end
	l(Number[6])
end

@example :FUNC_COLON begin
	eval(parse("function f11(): a = 6;end;f11()"))
end

@example :FUNC_COLON begin
	eval(parse("function f12(): print();end;f12()"))
end

@example :FUNC_COLON begin
	eval(parse("function f13(): Base.print();end;f13()"))
end

@example :DICTS_ITER begin
	for l in {"l"=>6}
		l+6
	end
end

@example :IMPORT_LVL begin
	function l();import Base;end
	l()
end

@example :IS_DEFINED begin
	isdefined(isdefined)
end

@example :IS_DEFINED begin
	isdefined("l")
end

@example :IS_DEFINED begin
	isdefined(l)
end

@example :IF_BOOLEAN begin
	if 6?6:6
	end
end

@example :IF_BOOLEAN begin
	if ""?"":""
	end
end

@example :IF_BOOLEAN begin
	if []?[]:[]
	end
end

@example :IF_BOOLEAN begin
	if 1 & 2 | 3
	end
end

@example :IF_BOOLEAN begin
	if length([])
	end
end

@example :PUSHAPPEND begin
	l = Int[]
	push!(l,[6])
end

@example :PUSHAPPEND begin
	l = Int[]
	append!(l,6)
end

@example :INDEX_BASE begin
	[][0]
end

@example :STRCONVERT begin
	"l"*6
end

@example :RANGEARRAY begin
	function l(::Array);end
	l(6:6)
end

@example :KEYWORDARG begin
	function l(q,w=6;e=6);end
	l(6,w=6)
end

@example :TYPEASSERT begin
	6::String
end

@example :ISOPERATOR begin
	6<:String
end

@example :SYMBOLEXPR begin
	function l(::Expr);end
	l(:(6))
end

@example :SYMBOLEXPR begin
	function l(::Symbol);end
	l(:(6+6))
end

@example :DOT_EQUALS begin
	eval(parse("[].===6"))
end

@example :TYPE_PARAM begin
	eval(parse("function f33{::Any}();end"))
end

@example :INDIRECTFN begin
	function l(::Any);end
	function li(a::Any)
		l() = l(a)
		l()
	end
	li(6)
end

@example :COLLECTSUM begin
	sum(values(Dict()))
end

@example :DIMS_MATCH begin
	ones((4,5)) * ones((4,5))
end

@example :BOOLASSIGN begin
	eval(parse("if a=0;end"))
end

@example :MAKESTRING begin
	String()
end

@example :MAKESYMBOL begin
	Symbol()
end


text:
include("agent/DQN.jl")
include("agent/OnlineDQN.jl")
include("agent/AT-DQN.jl")


text:
"""
    ArrayDictionary(indices, values)

Construct an `ArrayDictionary <: AbstractDictionary` from two arbitrary Julia iterables, one
specifying the `indices` and one specifying the `values`. Lookup uses naive iteration.

If `indices` and `values` are `Vector`s, then the result is `isinsertable` and
`istokenizable`, making simple and flexible dictionaries using naive search that may be
optimal for small collections.
"""
struct ArrayDictionary{I, T, Inds <: ArrayIndices{I}, Vals <: AbstractArray{T}} <: AbstractDictionary{I, T}
	indices::Inds
    values::Vals

    @inline function ArrayDictionary{I, T, Indices, Values}(indices::Indices, values::Values) where {I, T, Indices <: ArrayIndices{I}, Values <: AbstractArray{T}}
        @boundscheck if LinearIndices(parent(indices)) != LinearIndices(values)
            error("Dictinary indices and values inputs do not match")
        end
        return new(indices, values)
    end
end

"""
    ArrayDictionary{I,T}()

Construct an `ArrayDictionary` from empty `Vector`s, with `I` and `T` default to `Any`.
"""
@propagate_inbounds ArrayDictionary() = ArrayDictionary{Any, Any}([], [])
@propagate_inbounds ArrayDictionary{I}() where {I} = ArrayDictionary{I, Any}(I[], [])
@propagate_inbounds ArrayDictionary{I, T}() where {I, T} = ArrayDictionary{I, T}(I[], T[])

@propagate_inbounds ArrayDictionary(inds, vals) = ArrayDictionary(ArrayIndices(inds), vals)
@propagate_inbounds ArrayDictionary{I}(inds, vals) where {I} = ArrayDictionary{I}(ArrayIndices{I}(inds), vals)
@propagate_inbounds ArrayDictionary{I, T}(inds, vals) where {I, T} = ArrayDictionary{I, T}(ArrayIndices{I}(inds), vals)

@propagate_inbounds ArrayDictionary(inds::ArrayIndices, vals) = ArrayDictionary(inds, collect(vals))
@propagate_inbounds ArrayDictionary{I}(inds::ArrayIndices{I}, vals) where {I} = ArrayDictionary{I}(inds, collect(vals))
@propagate_inbounds ArrayDictionary{I, T}(inds::ArrayIndices{I}, vals) where {I, T} = ArrayDictionary{I, T}(inds, collect(T, vals))

@propagate_inbounds ArrayDictionary(inds::ArrayIndices, vals::AbstractArray) = ArrayDictionary{eltype(inds), eltype(vals)}(inds, vals)
@propagate_inbounds ArrayDictionary{I}(inds::ArrayIndices{I}, vals::AbstractArray) where {I} = ArrayDictionary{I, eltype(vals)}(inds, vals)

@propagate_inbounds ArrayDictionary{I,T}(inds::ArrayIndices{I}, vals::AbstractArray) where {I,T} = ArrayDictionary{I,T}(inds, convert(AbstractArray{T}, vals))
@propagate_inbounds ArrayDictionary{I,T}(inds::ArrayIndices{I}, vals::AbstractArray{T}) where {I,T} = ArrayDictionary{I,T,typeof(inds),typeof(vals)}(inds, vals)

"""
    ArrayDictionary(indices, undef::UndefInitializer)

Construct a `ArrayDictionary` from an iterable of `indices`, where the values are
undefined/uninitialized.
"""
@propagate_inbounds ArrayDictionary{I, T}(inds, ::UndefInitializer) where {I, T} = ArrayDictionary{I, T}(ArrayIndices{I}(inds), undef)
@propagate_inbounds ArrayDictionary{I, T}(inds::ArrayIndices{I}, ::UndefInitializer) where {I, T} = ArrayDictionary{I, T}(inds, similar(parent(inds), T))

"""
    ArrayDictionary(indexable)

Construct a `ArrayDictionary` from an indexable container `indexable` with the same `keys` and
`values`, equivalent to `ArrayDictionary(keys(indexable), values(indexable))`. Note that
`indexable` may not be copied.
"""
ArrayDictionary(indexable) = ArrayDictionary(keys(indexable), values(indexable))
ArrayDictionary{I}(indexable) where {I} = ArrayDictionary{I}(keys(indexable), values(indexable))
ArrayDictionary{I, T}(indexable) where {I, T} = ArrayDictionary{I, T}(keys(indexable), values(indexable))

Base.parent(d::ArrayDictionary) = getfield(d, :values)

function Base.keys(d::ArrayDictionary{I}) where {I}
    return getfield(d, :indices)
end

# token interface
istokenizable(::ArrayDictionary) = true
istokenassigned(d::ArrayDictionary, t::Int) = isassigned(parent(d), t)
@propagate_inbounds gettokenvalue(d::ArrayDictionary, t::Int) = parent(d)[t]

# settable interface
issettable(::ArrayDictionary) = true # Need an array trait for this...

@propagate_inbounds function settokenvalue!(d::ArrayDictionary{<:Any, T}, t::Int, value::T) where {T}
    parent(d)[t] = value
    return d
end

function Base.similar(inds::ArrayIndices, ::Type{T}) where {T}
    return ArrayDictionary(inds, similar(parent(inds), T))
end

# insertable interface
isinsertable(::ArrayDictionary) = true # Need an array trait for this...

@propagate_inbounds function gettoken!(d::ArrayDictionary{I}, i::I) where {I}
    (hadtoken, token) = gettoken!(keys(d), i)
    if !hadtoken
        resize!(parent(d), length(parent(d)) + 1)
    end
    return (hadtoken, token)
end

@propagate_inbounds function deletetoken!(d::ArrayDictionary, t::Int)
    deletetoken!(keys(d), t)
    deleteat!(parent(d), t)
    return d
end

function Base.empty!(d::ArrayDictionary)
    empty!(keys(d))
    empty!(parent(d))
    return d
end

function Base.empty(d::ArrayDictionary, ::Type{I}, ::Type{T}) where {I, T}
    return ArrayDictionary(empty(keys(d), I), empty(parent(d), T))
end


text:
@defcomp emissions begin
    regions     = Index()                           # The regions index must be specified for each component

    E           = Variable(index=[time, regions])   # Total greenhouse gas emissions
    E_Global    = Variable(index=[time])            # Global emissions (sum of regional emissions)
    sigma       = Parameter(index=[time, regions])  # Emissions output ratio
    YGROSS      = Parameter(index=[time, regions])  # Gross output - Note that YGROSS is now a parameter

    # function init(p, v, d)
    # end

    function run_timestep(p, v, d, t)
        # Define an equation for E
        for r in d.regions
            v.E[t,r] = p.YGROSS[t,r] * p.sigma[t,r]
        end

        # Define an equation for E_Global
        for r in d.regions
            v.E_Global[t] = sum(v.E[t,:])
        end
    end

end


text:
## Get a QMolecule object starting with the pyscf driver,
## all within qiskit_nature

@pyimport qiskit_nature.drivers as drivers

function get_molecule_nature()
    atom = "H .0 .0 .0; H .0 .0 0.7414"
    #    atom = "Li .0 .0 .0; H .0 .0 1.4"
#    atom = "O 0. 0. 0.; H 0.757 0.586 0.; H -0.757 0.586 0."
    driver = drivers.PySCFDriver(atom=atom,
                         unit=drivers.UnitsType.ANGSTROM,
                         charge=0,
                         spin=0,
                         basis="sto3g")
    molecule = driver.run()
    return molecule
end

nat_mol = get_molecule_nature()


text:
test = """
    mxmxvkd kfcds sqjhc nhms (contains dairy, fish)
    trh fvjkl sbzzf mxmxvkd (contains dairy)
    sqjhc fvjkl (contains soy)
    sqjhc mxmxvkd sbzzf (contains fish)
    """

function parse_input(io)
    foods = []
    while !eof(io)
        line = readline(io)
        (ingredients, allergens) = split(line, " (contains ")
        ingredients = Symbol.(split(ingredients, r"\s+"; keepempty = false))
        allergens = isempty(allergens) ? String[] :
            Symbol.(split(rstrip(==(')'), allergens), ", "))
        push!(foods, (ingredients, allergens))
    end
    return foods
end

function q1(foods)::Int
    assmt = assignment(foods)
    allergic = values(assmt)
    ingredients = vcat(first.(foods)...)
    return deleteat!(ingredients, findall(∈(allergic), ingredients)) |> length
end

function q2(foods)::String
    assmt = sort(collect(assignment(foods)))
    return join(last.(assmt), ',')
end

function assignment(foods::AbstractVector)::Dict
    assmt = Dict{Symbol, Vector{Symbol}}()
    for (ingredients, allergens) in foods
        for allergen in allergens
            if !haskey(assmt, allergen)
                assmt[allergen] = copy(ingredients)
            else
                intersect!(assmt[allergen], ingredients)
            end
        end
    end

    n_unambig = 0
    unambig = findall(==(1) ∘ length, assmt)
    while n_unambig != length(unambig)
        covered = [only(v) for (k, v) in assmt if k ∈ unambig]
        for k in setdiff(keys(assmt), unambig)
            setdiff!(assmt[k], covered)
        end
        n_unambig = length(unambig)
        unambig = findall(==(1) ∘ length, assmt)
    end

    return Dict(k => only(v) for (k, v) in assmt)
end

let foods = parse_input(IOBuffer(test))
    @assert q1(foods) == 5
    @assert q2(foods) == "mxmxvkd,sqjhc,fvjkl"
end

open("day21.in", "r") do io
    foods = parse_input(io)
    println("Q1: ", q1(foods))
    println("Q2: ", q2(foods))
end

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
############################
## DataFrame construction ##
############################

function composeDataFrame(vals, nams)
    ## compose DataFrame from data and column names
    ##
    ## Input:
    ## 	vals 		Array{Float64,2} of values
    ## 	nams		column names as Array{Symbol,1} as returned by
    ## 				names(df)
    ##
    ## Output:
    ## 	DataFrame with values given by vals and column names given by
    ## 	nams

    errMsg = "to compose a DataFrame, number of columns and number of
names must match"

    if size(vals, 2) != length(nams)
        println("size of values: $(size(vals))")
        println("number of names: $(length(nams))")
        error(errMsg)
    end

    df = convert(DataFrame, vals)
    names!(df, nams)
    return df
end

import Base.convert
function convert(::Type{DataFrame}, td::AbstractTimedata)
    df = [DataFrame(idx = td.idx) td.vals]
end

############################
## round dataframe values ##
############################

function tryRound(x)
    try
        round(x, 2)
    catch
        x
    end
end

function tryRound(x, nDgts::Int)
    try
        round(x, nDgts)
    catch
        x
    end
end

import Base.round
function round(df::DataFrame, nDgts::Int)
    vals = Any[tryRound(df[ii, jj], nDgts) for ii=1:size(df, 1),
        jj=1:size(df, 2)]
    return composeDataFrame(vals, names(df))
end

function round(df::DataFrame)
    vals = Any[tryRound(df[ii, jj]) for ii=1:size(df, 1),
        jj=1:size(df, 2)]
    return composeDataFrame(vals, names(df))
end

macro roundDf(expr::Expr)
    quote
        res = $(esc(expr))
        resRnd = Any[tryRound(res[ii, jj]) for ii=1:size(res, 1), jj=1:size(res, 2)] |>
                 x -> composeDataFrame(x, names(res))

        display(resRnd)
    end
end

##################################
## basic mathematical operators ##
##################################

## importall Base
## for op = (:.+, :.-, :.*, :./, :.^)
##   eval(quote
##       function $op(df1::DataFrame, df2::DataFrame)
##           try
##               (array(df1), array(df2))
##               catch
##               error("mathematical operators only work on DataFrames
## with numeric values only")
##           end

##           res = $op(array(df1), array(df2))
##           ## if size(res, 1) == 1
##               ## df = composeDataFrame(res', names(df1))
##           ## else
##           df = composeDataFrame(res, names(df1))
##           ## end
##       end
##   end)
## end

##############
## setDfRow ##
##############

function setDfRow!(df::DataFrame, arr::Array, rowInd::Int)
    nElem = length(arr)
    if nElem != size(df, 2)
        error("wrong number of columns during assignment")
    end

    for ii=1:nElem
        df[rowInd, ii] = arr[ii]
    end
    return df
end

function setDfRow!(df::DataFrame, dfRow::DataFrame, rowInd::Int)
    if size(df, 2) != size(dfRow, 2)
        error("wrong number of columns during assignment")
    end
    if size(dfRow, 1) != 1
        error("only 1xn DataFrame can be used for row assignment")
    end

    for ii=1:size(df, 2)
        df[rowInd, ii] = dfRow[1, ii]
    end
    return df
end

function setDfRow!(df::DataFrame, dfRow::DataFrameRow, rowInd::Int)
    if size(df, 2) != length(dfRow)
        error("wrong number of columns during assignment")
    end

    for ii=1:size(df, 2)
        df[rowInd, ii] = dfRow[ii]
    end
    return df
end

######################
## org-babel output ##
######################

function Base.writedlm(io::IO, df::DataFrame, dlm; opts...)
    nams = names(df)
    vals = Any[df[ii, jj] for ii=1:size(df, 1), jj=1:size(df, 2)]
    arrAny = [nams', vals]
    writedlm(io, arrAny, dlm; opts...)
end

###################
## preallocation ##
###################

function preallocateDf(td)
    ## DataFrame of same size with same column types as TimeData
    return DataFrame(eltypes(td.vals), names(td), size(td, 1))
end

function preallocateBoolDf(nams::Array{Symbol,1}, nObs::Int)
    nVars = length(nams)
    typs = [Bool for ii=1:nVars]
    return DataFrame(typs, nams, nObs)
end


text:
using Test
using Krippendorff
using Tables

# write some tests
# should do some unit tests and then some tests against known reliability values

# krippendorffs_alpha

# examples from: "Computing Krippendorff’s Alpha-Reliability, Klaus Krippendorff, kkrippendorff at asc.upenn.edu, 2011.1.25"
binarytest = [ 0 1 0 0 0 0 0 0 1 0   # α ≈ 0.095 in his example
               1 1 1 0 0 1 0 0 0 0 ]
@test round(Krippendorff.alpha(binarytest, units = :cols), sigdigits = 2) == 0.095
@test krippendorffs_alpha(binarytest', units = :rows) == Krippendorff.alpha(binarytest, units = :columns)
# the next test passes only with approx, indicating some minor difference between the approaches
@test Krippendorff.alpha(binarytest, R = 0:1) ≈ Krippendorff.alpha(binarytest, R = :continuous)

# _dispatch_units_iterator
# _units_eltype_and_size
# _units_eltype_and_size(Matrix)
# _dispatch_responses(precomputed R)
# _dispatch_responses(Symbol/String)
# _optimize_possible_reponses
# _map_responses_to_indices(mapping)
# _dispatch_metric(Symbol/String)
@test Krippendorff._dispatch_metric("nominal", nothing) == (!=)
kripp_intervall = Krippendorff._dispatch_metric(:interval, Int)
test_intervall = (x,y)->abs2(x-y)
intervaltest_rands = rand(Int,(100,2))
@test all([kripp_intervall(rands[1],rands[2]) for rands in eachrow(intervaltest_rands)] .== [test_intervall(rands[1],rands[2]) for rands in eachrow(intervaltest_rands)])
@test_throws ArgumentError Krippendorff._dispatch_metric(:blob, Int)
testmetric = (x,y)->(x-y)
@test Krippendorff._dispatch_metric(testmetric, Int) == testmetric
# _compute_distance_matrix
# _compute_alpha
# _fill_coincidence_matrix
# _evaluate_coincidence_matrix
# compute_alpha_generical
# compute_alpha_with_coincidences
# prepare_iterator
@test all(iterate(Krippendorff.prepare_iterator(reshape(1:20,(4,5))))[1] .== [1,5,9,13,17])
@test collect(Iterators.flatten(Krippendorff.prepare_iterator([[1,2,missing],[missing,missing,missing],[3,missing]]))) == [1,2,3]
# istable

max_stars_repo_name: FPGro/Krippendorff

max_issues_repo_name: FPGro/Krippendorff

max_forks_repo_name: FPGro/Krippendorff


text:
using Genie.Configuration

const config = Settings(
  server_port                     = 8000,
  server_host                     = "127.0.0.1",
  log_level                       = Logging.Debug,
  log_to_file                     = true,
  server_handle_static_files      = true,
  websocket_server                = false
)

ENV["JULIA_REVISE"] = "off"


text:
## A set of functions for exploring convergence
## in an interactive manner:

### Benchmarking tests
## We have
##
## * visualize_diagonostics(which): to see summaries of the methods
## over the different functions with which in (:summary, :counts,
## :residuals)
##
## * identify_regression(): to identify regressions in counts,
## residuals, failures, or exact answers, as compared with a past
## diagnostic run
##
## * write_out(): to save results of a diagnostic run. Must be run in the Roots/test directory
##
## * compare_convergence(Methods): to generate a table of convergence orders

using Roots
using Test
using JSON
using Printf

## This uses functions defined here:
include(joinpath(@__DIR__, "test_derivative_free.jl"))

mutable struct Wrapper
f
n
end
Wrapper(f) = Wrapper(f, 0)
(F::Wrapper)(x) = (F.n+=1; F.f(x))

struct TallyCard
cnts::Vector{Int}     # fn evals
res::Vector{Float64} # rts
rs::Vector{Float64}   # residual
exact::Vector{Int}    # exact 1, 0, -1 (for fail)
ds::Vector{Float64}   # difference from alpha(s)
rng::Vector{Float64}  # range sampled
end
TallyCard() = TallyCard(Int[], Float64[], Float64[], Int[],Float64[], Float64[])

# test a single function, create tally card
function test_function(method, F)
    tc = TallyCard()
    m,M = Inf, -Inf
    for p in F.params
        f = Wrapper(x->F.val(p,x))
        result, residual = nothing, nothing
        x0 = F.x0(p)
        m = x0 < m ? x0 : m
        M = x0 > M ? x0 : M
        try
            result = method(f, x0)
            cnt = f.n
            residual = f(result)
            d = minimum(abs.(result .- F.alpha))
            exact = (iszero(residual) || f(result)*f(nextfloat(result)) < 0 || f(result) * f(prevfloat(result)) < 0) ? 1 : 0

            push!(tc.cnts, cnt)
            push!(tc.res, result)
            push!(tc.ds, d)
            push!(tc.rs, residual)
            push!(tc.exact, exact)
        catch err
            push!(tc.cnts, -1)
            push!(tc.res, NaN)
            push!(tc.ds, NaN)
            push!(tc.rs, NaN)
            push!(tc.exact, -1)
        end
    end
    append!(tc.rng, (m, M))
    tc
end

# vector of vectors to an array
function vvta(vs, T=eltype(vs[1]))
    n = length(vs)
    m = length(vs[1])
    A = zeros(T, n, m)
    for j in 1:n
        A[j,:] = vs[j]
    end
    A
end

function vvta1(vs, T)
    n = length(vs)
    m = length(vs[1])
    A = zeros(T, m,n)
    for j in 1:n
        A[:,j] = vs[j]
    end
    A
end

## Return Dict of arrays
function create_diagonostics(Ms, Fs, nms)
    @assert length(nms) == length(Ms)

    out = Array{Any}(undef, length(Ms), length(Fs))

    for (i,M) in enumerate(Ms)
        for (j,F) in enumerate(Fs)
            out[i,j] = test_function(M, F)
        end
    end

    m,n = length(Ms), length(Fs)

    D = Dict("counts" =>Dict(),
             "rts" => Dict(),
             "residuals" =>Dict(),
             "exact" => Dict(),
             "delta" => Dict(),
             "failed"=>Dict(),
             "nms" => nms,
             "size" => (m,n))

    # Counts
    for j in eachindex(Fs)
        jj = string(j) # for JSON serialization
        cnts = [out[i,j].cnts for i in eachindex(Ms)]
        rts = [out[i,j].res for i in eachindex(Ms)]
        rs = [out[i,j].rs for i in eachindex(Ms)]
        exact = [out[i,j].exact for i in eachindex(Ms)]
        ds = [out[i,j].ds for i in eachindex(Ms)]

        M = length(cnts)
        N = div(M,m)

        D["counts"][jj] = vvta(cnts)
        D["rts"][jj] = vvta(rts)
        D["residuals"][jj] = vvta(rs)
        D["exact"][jj] =  vvta(exact)
        D["delta"][jj] =  vvta(ds)

        # failed
        fc = D["counts"][jj] .< 0
        fd = D["delta"][jj] .> 1.0
        fr = D["residuals"][jj] .> 1.0
        D["failed"][jj] = [fc[i,j] || fd[i,j] || fr[i,j] for i in 1:size(fc)[1], j in 1:size(fc)[2]]
    end

    D
end

function visualize_diagnostics(::Val{:summary}, D)

    ## Some diagnostics
    n = length(D["counts"]) # number of functions
    m = size(D["counts"]["1"])[1]

    for i in 1:m
        fs, cnt, exact, maxresidual = 0, 0, 0, 0.0
        for j in 1:n
            jj = string(j)

            fail = D["failed"][jj][i,:]
            cnts =  D["counts"][jj][i,:]
            resids = D["residuals"][jj][i,:]
            exacts = D["exact"][jj][i,:]

            fs += sum(fail)
            cnt += sum(cnts[.!fail])
            exact += sum(exacts .== 1)
            rs = resids[.!fail]
            if !isempty(rs)
                maxresidual = max(maxresidual, maximum(rs))
            end
        end
        nm = D["nms"][i]
        print(rpad(nm, 15))
        print(rpad(cnt, 6))
        print(rpad(exact, 6))
        print(rpad(fs, 5))
        println(maxresidual)
    end

end

function visualize_diagnostics(::Val{:counts}, D, j)

    ## Some diagnostics
    n = length(D["counts"]) # number of functions
    m = size(D["counts"]["1"])[1]

#    for j in 1:n
        jj = string(j)
        counts = D["counts"][jj]
        fails = D["failed"][jj]
        println(" --- Function $jj ---")
        for i in 1:m
            nm = D["nms"][i]
            print(rpad(nm, 15))
            tot,ntot = 0,0
            for (fail, cnt) in zip(fails[i,:], counts[i,:])
                if fail
                    print(lpad(".", 5))
                else
                    tot += cnt
                    ntot += 1
                    print(lpad(cnt, 5))
                end
            end
            avg = round(tot/ntot, digits=1)
            println("\t$avg")
        end
#    end
end

function visualize_diagnostics(::Val{:residuals}, D, j)

    ## Some diagnostics
    n = length(D["counts"]) # number of functions
    m = size(D["counts"]["1"])[1]

    jj = string(j)
    resids = D["residuals"][jj]
    exacts = D["exact"][jj]
    fails = D["failed"][jj]
    println(" --- Function $jj ---")
    for i in 1:m
        nm = D["nms"][i]
        print(rpad(nm, 15))
        tot,ntot = 0,0
        for (fail, res, exact) in zip(fails[i,:], resids[i,:], exacts)
            if fail
                ch = "X"
            elseif res > 1e-3
                ch = "x"
            elseif res > 1e-12
                    ch = "~"
            elseif exact == 1
                ch = "∘"
            else
                ch = "."
            end
            print(lpad(ch, 2))
        end
        println("")
    end
end

# which in (:all, :summary, :counts, :residuals)
function visualize_diagnostics(D, which=:summary) # :cnts, ...

    ## Some diagnostics
    n = length(D["counts"]) # number of functions
    m = size(D["counts"]["1"])[1]

    println("Visualize diagnostics: :$which ∈ (:all, :summary, :counts, :residuals)\n")

    # summary
    if which in (:all, :summary)
        println("Method         fncnt exact fail maxresidual")
        # one row summarizing each method
        visualize_diagnostics(Val(:summary), D)
    end

    # Counts
    if which in (:all, :counts)
        for j in 1:n
            visualize_diagnostics(Val(:counts), D,j)
        end
    end

    # residuals
    if which in (:all, :residuals)

        println("Residual summary")
        println("Key:\t`∘` - exact\n\t`.` - res < 1e-12\n\t`~` - 1e-12 < res < 1e-3`\n\t`x` - 1e-3 < res\n\t`X` - failed")

        for j in 1:n
            visualize_diagnostics(Val(:residuals), D,j)
        end
    end

end

## write out/ read in summaries to file
function write_out(fname, D)
    io = open(fname, "w")
    JSON.Writer.print(io, D)
    close(io)
end

# annoyingly need to convert values to proper type after
# JSON serializations
function read_in(fname)
    D = JSON.parsefile(fname)
    E = Dict()
    E["nms"] = string.(D["nms"])
    E["size"] = Tuple(D["size"])

    for nm in ("counts", "exact", "failed")
        E[nm] = Dict()
        for (k,v) in D[nm]
            E[nm][k] = vvta1(v, Int)
        end
    end

    E["residuals"] = Dict()
    for (k,v) in D["residuals"]
        for vi in v
            vi[vi .=== nothing] .= NaN
        end
        E["residuals"][k] = vvta1(v, Float64)
    end

    E
end

# compare D to E
function identify_regressions(Dnew,Dold)
    out = String[]

    Dnew["nms"] == Dold["nms"] || return "Names are different"
    Dnew["size"] == Dold["size"] || return "sizes are different"

    for (k, v) in Dnew["counts"]
        A, A1 = v, Dold["counts"][k]
        sum(A[A .> 0]) <= sum(A1[A1 .> 0]) || push!(out, "counts increased for function $k")
        sum(A[A .> 0]) < sum(A1[A1 .> 0]) && push!(out, "✓ counts decreased for function $k")
    end

     for (k, v) in Dnew["exact"]
         A, A1 = v, Dold["exact"][k]
         sum(A .== 1) >= sum(A1 .== 1)  || push!(out, "exact decreased for function $k")
         sum(A .== 1) > sum(A1 .== 1)  && push!(out, "✓ exact increased for function $k")
    end

    for (k, v) in Dnew["failed"]
        A, A1 = v, Dold["failed"][k]
        sum(A) <= sum(A1)  || push!(out, "failed increased for function $k")
        sum(A) < sum(A1) && push!(out, "✓ failed decreased for function $k")
    end

    for (k, v) in Dnew["residuals"]
        A, A1 = v, Dold["residuals"][k]
        for i in eachindex(A)
            newi, oldi = A[i], A1[i]
            if abs(newi) > 1.1 * abs(oldi)
                push!(out, "residuals increased for function $k")
                break
            end
            if abs(newi) < 0.9 * abs(oldi)
                push!(out, "✓ residuals decreased for function $k")
                break
            end
        end
    end

    return out

end

## Main interface for interactive use
fname = joinpath(@__DIR__, "derivative_free_diagnostics.json")
elide_ascii(x, n=12) = length(x) > n ? x[1:n-3] * "..." * x[end-1:end] : x
function create_diagonostics()
    meths = [Order0(),
             Order1(), Roots.Order1B(), Roots.King(),
             Order2(), Roots.Steffensen(), Roots.Order2B(), Roots.Esser(),
             Order5(), Roots.KumarSinghAkanksha(),
             Order8(), Roots.Thukral8(),
             Order16(), Roots.Thukral16()]

    Ms = [(f, b) -> find_zero(f, b, M) for M in meths] # F(f,b), name
    nms = elide_ascii.([replace(string(M), r"^Roots."=>"") for M in meths])
    Fs = known_functions
    create_diagonostics(Ms, Fs, nms)
end
## write out current diagnostic test
function write_out()
    println("Creating diagonostics to save")
    write_out(fname, create_diagonostics())
end

## visualize state
"""
   visualize_diagnostics(which=:summary)

Show diagnostics summary

`which` is one of  `(:all, :summary, :counts, :residuals)`
"""
visualize_diagnostics(which=:summary) = visualize_diagnostics(create_diagonostics(), which)

## identify regressions from currently saved state
"""
    identify_regressions()

Compare current state to saved state.

Use `write_out` to save a state.
"""
function identify_regressions()

    if !isfile(fname)
        println("No previous diagnostic to compare with")
        return String[] # empty
    end

    Dnew = create_diagonostics()
    Dold = read_in(fname)
    out = identify_regressions(Dnew, Dold)

    out
end

## tests for newton, halley
import ForwardDiff: derivative
D(f,n=1) = n > 1 ? D(D(f), n-1) : x -> derivative(f, float(x))
function derivative_based_diagonostics()
    Ms = ((f,b) -> Roots.find_zero((f, D(f)), b, Roots.Newton()),
          (f,b) -> Roots.find_zero((f, D(f),D(f,2)), b, Roots.Halley()),
          (f,b) -> Roots.find_zero((f, D(f),D(f,2)), b, Roots.Schroder()),
          (f,b) -> Roots.find_zero((f, D(f)), b, Roots.Order5())
          )
    nms = ("Newton", "Halley","Schroder","Order5")
    Fs = known_functions
    create_diagonostics(Ms, Fs, nms)
end

## Order of convergence
## assuming e_n = x_n - alpha and
## e_{n+1} = C e_n^q this tries to find q by:
## alpha known: e_{n+1}/e_n = (e_n/e_{n-1})^q,
## alpha unknown: f(x_n) - f(alpha) ~ f'(alpha)*(en), so en ~ Cf(x_n)
function COC(M, f, x0, alpha=missing)
    op = precision(BigFloat)
    setprecision(BigFloat, 8 * 256)
    tracks = Roots.Tracks(BigFloat[], BigFloat[])
    try
        find_zero(f, big(x0), M, tracks=tracks)
    catch err
        setprecision(BigFloat, op)
        rethrow()
    end
    setprecision(BigFloat, op)

    if ismissing(alpha)
        fs = tracks.fs
        [Float64(log(abs(fs[k+2]/fs[k+1])) /   log(abs(fs[k+1]/fs[k]))) for k in 1:length(fs)-4]
    else
        xs = tracks.xs
        es = xs .- alpha
        [Float64(log(abs(es[k+1]))/log(abs(es[k]))) for k in 1:length(xs)-3]
    end
end

### Traditional tests: start nearby and compare convergence
"""

Compare convergences. For each method and for several test functions, computes
* abs(x1-alpha), abs(x2-alpha), abs(x3-alpha);
* the computational order of convergence rc = log|f(xk+1)/f(xk)|/log|f(xk)/f(xk-1)|
* the computed zero
Example
```
compare_convergence((Order1(), Order2(), Order8()))
##
## Example (errors from Sidi, Unified treatment...)
## Let error be En = xn - alpha
## Secant: E_{n+1} = f[x_n, x_{n-1}, alpha] / f[x_n, x_{n-1}] E_n E_{n-1}
## Newton: E_{n+1} = f[xn, xn, alpha]/f'(xn) E_n^2
## Steffensen: E_{n+1} = f[xn, xn+fxn, alpha] / f[xn, xn + fxn] ⋅ (1 f[xn,alpha] ⋅ En^2
using ForwardDiff; D(f) = x -> ForwardDiff.derivative(f, float(x))
struct MNewton <: Roots.AbstractSecant end
function Roots.update_state(M::MNewton, f, o, opts)
    o.xn0, o.fxn0 = o.xn1, o.fxn1
    o.xn1 -= f(o.xn1) / D(f)(o.xn1)
    o.fxn1 = f(o.xn1)
    nothing
end
compare_convergence((Roots.Secant(), MNewton(), Roots.Steffensen()))

```
"""
function compare_convergence(Ms; F=identity)
    fns = ((x -> log(x^2 + x + 2) - x + 1, big"4.4",
            4.152590736757158274996989004767139785813809448259893154635015805935085336704627),
           (x -> sin(x)^2-x^2 + 1, big"1.4",
            1.404491648215341226035086817786868077176602575918625035145218238569654850906246),
           (x -> exp(-x^2 + x+2)-cos(x+1)+x^3 + 1, -big"0.5", -big"1.0"),
           (x -> x^11 + x +1, -big"1",
            -8.4439752879202298306029802319765064181966469706324549911391920468474414245176e-01),
    (x ->  (x-2)*(x^10 + x + 1)*exp(-x-1), big"1.9", big"2"),

    (x->sin(3x) + x*cos(x), big"1",
     1.197769535216271165938579472950989827411047786536025790115116815210444571657156),
    (x -> exp(-x) - 1 + x/5, big"4.5",
     4.96511423174427630369875913132289394405558498679725097281444614478046398795746),
    (x -> exp(sin(x)) - x + 1, big"2.3",
     2.630664147927903633975327052350598568584731954733163386430717083451519883744738),
    (x -> (exp(x^2 +x - 6)-1)*(x^2+1), big"2.2", big"2"),
    (x -> x*exp(x^2) - sin(x)^2 + 3cos(x) + 5, -big"1.2",
     -1.207647827130918927009416758356084097760235818949538815205924601763336168539887),

    (x -> x^8 - x^7 + x^4 - x^3 + x -1, big"1.2", big"1")
    )

    # formatting fns
    Log10 = x -> begin # formatting
        if iszero(x)
            return (0,0)
        end
        n = trunc(Int,log10(x))-1
        x/(10.0^n), -n
    end
    elide_ascii(x, n=12) = length(x) > n ? x[1:n-2] * ".." * x[end-1:end] : x

    for (i,u) in enumerate(fns)
        fn_,x0,xstar = u
        fn = F(fn_)
        for M in Ms
            tracks = Roots.Tracks(BigFloat, BigFloat)
            a = try
                find_zero(fn, x0, M, tracks=tracks)
            catch err
                NaN * x0
            end
            x1,x2,x3 = tracks.xs[1:3]

            coc = log(abs(fn_(x3)/fn_(x2))) / log(abs(fn_(x2)/fn_(x1)))
            X1,X2,X3 = abs.((x1,x2,x3) .- xstar)

            @printf "%s: %s\t%.3f(-%02i)\t%.3f(-%02i)\t%.3f(-%02i)\t%.2f\t%.8f\n"  "f$i"  elide_ascii(replace(string(M), "Roots." => ""),8)  Log10(X1)... Log10(X2)... Log10(X3)... coc a

        end
        println("")
    end
end


text:
module BAInvariants

using StaticArrays

# ------------------------------------------------------------------------
#             2-BODY Invariants
#             fully equivalent to BL invariants
# ------------------------------------------------------------------------

# x = (r12,)

invariants(x::SVector{1, T}) where {T} =
      copy(x),
      SVector{1, T}(1.0)

invariants_d(x::SVector{1, T}) where {T} =
      (@SVector [ SVector(one(T))  ]),
      (@SVector [ SVector(zero(T)) ])

invariants_ed(x::SVector{1,T}) where {T} =
      copy(x),
      SVector{1, T}(1.0),
      (@SVector [ SVector(one(T))  ]),
      (@SVector [ SVector(zero(T)) ])

tdegrees(::Val{2}) = (1,), (0,)

# this is a terrible hack that mostly likely will not work as we want.
# see Slack discussion for more details.
rtdegrees(::Val{2}) =
   ( (1,0), ), ( (0,0), )

# ------------------------------------------------------------------------
#             3-BODY Invariants
# ------------------------------------------------------------------------

# x = (r1, r2, θ12)

# the 1.0 is a "secondary invariant"
invariants(x::SVector{3, T}) where {T} =
      (@SVector T[ x[1] + x[2], x[1] * x[2], x[3] ]),
      (@SVector T[ 1.0 ])

invariants_d(x::SVector{3, T}) where {T} =
      (@SVector [ (@SVector T[1.0, 1.0, 0.0]),
                  (@SVector T[x[2], x[1], 0.0]),
                  (@SVector T[0.0, 0.0, 1.0]) ]),
      (@SVector [ (@SVector T[0.0, 0.0, 0.0]) ])

invariants_ed(x::SVector{3, T}) where {T} =
      (@SVector T[ x[1] + x[2], x[1] * x[2], x[3] ]),
      (@SVector T[ 1.0 ]),
      (@SVector [ (@SVector T[1.0, 1.0, 0.0]),
                  (@SVector T[x[2], x[1], 0.0]),
                  (@SVector T[0.0, 0.0, 1.0]) ]),
      (@SVector [ (@SVector T[0.0, 0.0, 0.0]) ])

tdegrees(::Val{3}) = (1, 2, 1), (0,)

# this is a terrible hack that mostly likely will not work as we want.
# see Slack discussion for more details.
rtdegrees(::Val{3}) =
   ( (1,0), (2,0), (0,1) ), ( (0,0), )

# ------------------------------------------------------------------------
#             4-BODY Invariants
# ------------------------------------------------------------------------

include("generated/BA_4B/BA_4B_invariants_hack.jl")

# ------------------------------------------------------------------------
#             5-BODY Invariants
# ------------------------------------------------------------------------

include("generated/BA_5B/BA_5B_invariants_hack.jl")

end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
function wordcount(sentence::AbstractString)
    sentence = lowercase(sentence)
    words = Dict()
    for rx = eachmatch(r"[a-z]+'[a-z]+|[0-9a-z]+", sentence)
        word = rx.match
        words[word] = get(words, word, 0) + 1
    end
    words
end


text:
"""
    getpages(s)

Given an arbitray nested structure of Arrays of pairs of Arrays etc
This will return an array with all the markdown filenames
"""
getpages(s::String) = s
getpages(p::Pair)   = getpages(p.second)
getpages(a::Array)  = reduce(vcat,getpages.(a))


text:
using MetaStrategist
using Documenter

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

makedocs(;
    modules=[MetaStrategist],
    authors="Khalil CHRIT",
    repo="https://github.com/JuliaConstraints/MetaStrategist.jl/blob/{commit}{path}#{line}",
    sitename="MetaStrategist.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaConstraints.github.io/MetaStrategist.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaConstraints/MetaStrategist.jl",
    devbranch="main",
)


text:
# La ventana principal y una rejilla para acomodar los botones

g = @Grid()
win = @Window(g, "Gato", 300, 300, false, true)

# Los botones...

b1 = @Button("1")
b2 = @Button("2")
b3 = @Button("3")
b4 = @Button("4")
b5 = @Button("5")
b6 = @Button("6")
b7 = @Button("7")
b8 = @Button("8")
b9 = @Button("9")

# ... y sus posiciones

g[1, 1] = b1
g[1, 2] = b2
g[1, 3] = b3
g[2, 1] = b4
g[2, 2] = b5
g[2, 3] = b6
g[3, 1] = b7
g[3, 2] = b8
g[3, 3] = b9

# Propiedades para que los botones llenen la ventana, con una separación adecuada

setproperty!(g, :column_homogeneous, true)
setproperty!(g, :column_spacing, 15)
setproperty!(g, :row_homogeneous, true)
setproperty!(g, :row_spacing, 15)
showall(win)

# Señales cuando los botones son apretados

signal_connect(button_clicked_cb, b1, :clicked)
signal_connect(button_clicked_cb, b2, :clicked)
signal_connect(button_clicked_cb, b3, :clicked)
signal_connect(button_clicked_cb, b4, :clicked)
signal_connect(button_clicked_cb, b5, :clicked)
signal_connect(button_clicked_cb, b6, :clicked)
signal_connect(button_clicked_cb, b7, :clicked)
signal_connect(button_clicked_cb, b8, :clicked)
signal_connect(button_clicked_cb, b9, :clicked)

# Para correr el código como script

signal_connect(win, :destroy) do widget
	Gtk.gtk_quit()
end
Gtk.gtk_main()

max_stars_repo_name: rodrigolece/juego-gato

max_issues_repo_name: rodrigolece/juego-gato

max_forks_repo_name: rodrigolece/juego-gato


text:
function accuracy(y_pred, y)
    return mean(onecold(cpu(softmax(y_pred))) .== onecold(cpu(y)))
end


text:
using Images

function to_array(history::AutomataBinaryHistory)
    start_width = length(history.states[1])
    width = length(history.states[end])
    height = length(history)
    img = zeros(Float32, height, width)

    step = 1
    for (state, uns) = history
        offset = (width - start_width) ÷ 2 - step + 1
        img[step, :] .= uns
        for i = 1:length(state)
            img[step, i + offset] = state[i]
        end
        step += 1
    end

    return img
end
function to_array(history::AutomataColorHistory)
    start_width = length(history.states[1])
    width = length(history.states[end])
    height = length(history)
    img = zeros(Float32, height, width)

    to_color(s) = Float32(s) / (Float32(history.colors) - 1)

    step = 1
    for (state, uns) = history
        offset = (width - start_width) ÷ 2 - step + 1
        img[step, :] .= to_color(uns)
        for i = 1:length(state)
            img[step, i + offset] = to_color(state[i])
        end
        step += 1
    end

    return img
end
function to_array(history::AutomataContinuousHistory)
    start_width = length(history.states[1])
    width = length(history.states[end])
    height = length(history)
    img = zeros(Float32, height, width)

    step = 1
    for (state, uns) = history
        offset = (width - start_width) ÷ 2 - step + 1
        img[step, :] .= uns
        for i = 1:length(state)
            img[step, i + offset] = state[i]
        end
        step += 1
    end

    return img
end

export render
function render(filename::String, arr::Array{Float32, 2}, filters=[])
    for arr_filter = filters
        arr = arr_filter(arr)
    end
    save(filename, Gray.(1 .- arr))
end
function render(filename::String, history::T, filters=[]) where {T <: AutomataHistory}
    render(filename, to_array(history), filters)
end


text:
#loads the value weighted crsp index
#also does some light preprocessing
function loadidxcrsp(idxcrspname::String, retainedidxcrspcols::Vector{Symbol};
  idxcrsppath::String = IDXCRSP_PATH,
  idxcrspdateformat::DateFormat = IDXCRSP_DATE_FORMAT,
  idxdatecol = IDX_DATE_COL)

  idxcrsp::DataFrame = CSV.read("$idxcrsppath\\$idxcrspname.csv") |> DataFrame
  idxcrsp.date = (i::Int->Date("$i", idxcrspdateformat)).(idxcrsp[!,idxdatecol])

  TEST_OUTPUT && (idxcrsp |> CSV.write("$OUT_PATH\\$idxcrspname.csv"))
  select!(idxcrsp, retainedidxcrspcols)

  return idxcrsp
end

#load the ff data
function loadffdailyports(ffportname::String = FF_PORT_NAME,
  ffpath::String = FF_PATH,
  ffportrename::Dict = FF_PORT_RENAME,
  ffportvals::Vector{Symbol} = FF_PORT_VALS,
  ffdateformat::DateFormat = FF_DATE_FORMAT)

  ff::DataFrame = CSV.read("$ffpath\\$ffportname.csv") |> DataFrame
  ff.date = (i::Int->Date("$i", ffdateformat)).(ff.crspdate)

  names!(ff, (s::Symbol->Symbol(lowercase(replace(string(s)," "=>"")))).(names(ff)))
  rename!(ff, ffportrename)
  TEST_OUTPUT && (ff |> CSV.write("$OUT_PATH\\$ffportname.csv"))
  #now convert the returns to indices
  for (i,r) ∈ enumerate(eachrow(ff))
    if i == 1 #start the index
      for s ∈ ffportvals
        r[s] = 1.0
      end
      continue
    end

    #=Threads.@threads =#
    for s ∈ ffportvals
      v = r[s]

      #handle missing values- assume no change in the index if missing
      ((v == -99.99) || (v < -100.)) && (v = 0.0)
      r[s] = (1. + v/100.) * ff[i-1,s]
    end
  end

  #TEST_OUTPUT && (ff |> CSV.write("$OUT_PATH\\$ffportname.csv"))
  select!(ff, [:date; ffportvals;])

  #rename index fields to avoid confusion with actual return fields
  rename!(ff, Dict(s=>Symbol(s,:idx) for s ∈ ffportvals))

  return ff
end

#interpolates the index to every day in an interval using the appropriate compounding
function interpolateindex(df::DataFrame, idxcol::Symbol,
  datecol::Symbol = :date, tol::Float64 = 10^-10)

  local idf::DataFrame

  df = df[completecases(df[:, [idxcol, datecol]]), :]

  begindate::Date = minimum(df[!,datecol])
  enddate::Date = maximum(df[!,datecol])
  alldates::Vector{Date} = collect(begindate:Day(1):enddate)

  #now create the interpolated dataframe
  Nidf::Int = length(alldates)
  idf = DataFrame(rid = collect(1:Nidf))
  idf[!,datecol] = alldates
  idf = join(idf, df, on=datecol, kind=:left)
  sort!(idf, datecol)

  #set the values to be interpolated between
  T::Type = eltype(idf[!,idxcol])
  idf.nextvalid = Vector{T}(undef, Nidf)
  idf.prevvalid = Vector{T}(undef, Nidf)
  idf.nextvaliddate = Vector{MDate}(undef, Nidf)
  idf.prevvaliddate = Vector{MDate}(undef, Nidf)

  #previous values
  for (i,r) ∈ enumerate(eachrow(idf))
    (i==1) && continue

    if  ismissing(idf[i-1, idxcol])
      r.prevvalid = idf[i-1, :prevvalid]
      r.prevvaliddate = idf[i-1, :prevvaliddate]
    else
      r.prevvalid = idf[i-1, idxcol]
      r.prevvaliddate = idf[i-1, datecol]
    end
  end

  #nextvalues
  for i ∈ Nidf:-1:1
    if ismissing(idf[i,idxcol])
      idf[i, :nextvalid] = idf[i+1, :nextvalid]
      idf[i, :nextvaliddate] = idf[i+1, :nextvaliddate]
    else
      idf[i, :nextvalid] = idf[i, idxcol]
      idf[i, :nextvaliddate] = idf[i, datecol]
    end
  end

  idf.period = ((prev,next)->(!ismissing(prev)) && (!ismissing(next)) ?
    (next - prev).value : missing).(idf.prevvaliddate, idf.nextvaliddate)
  idf.ccret = (log).(idf.nextvalid ./ idf.prevvalid) ./ idf.period

  #finally do the interpolation
  for (i,r) ∈ enumerate(eachrow(idf))
    if ismissing(r[idxcol])
      r[idxcol] = exp(r.ccret)*idf[i-1, idxcol]
    elseif i > 1 #do an error check
      Δ::Float64 = abs(exp(r.ccret)*idf[i-1, idxcol] - r[idxcol])
      (Δ > tol) && @warn("Tolerance $tol exceeded with $Δ on row $i. Row: $r")
    end
  end

  idf |> CSV.write("$OUT_PATH\\idf-$idxcol.csv")

  #drop extraneous values
  select!(idf, Not([:nextvalid, :nextvaliddate, :prevvalid,
    :prevvaliddate, :ccret, :period]))

  return idf
end

function indexreturns!(df::DataFrame, idxcol::Symbol, retcol::Symbol;
  logret::Bool = true, sortdf::Bool = true,
  datecol::Symbol = :date, uselogs::Bool = USE_LOGS)

  #pre-allocate and sort if needed
  N = size(df,1)
  sortdf && (sort!(df, datecol))
  df[!,retcol] = Vector{MFloat64}(undef, N)

  for (i,r) ∈ enumerate(eachrow(df))
    (i==1) && continue #special case since we are working with lags

    #compute the return
    r[retcol] = r[idxcol]/df[i-1, idxcol] - 1.
    (uselogs) && (r[retcol] = log(1. + r[retcol]))
  end

  return df
end

function mergeoutside!(fact::DataFrame; ffportvals::Vector{Symbol} = FF_PORT_VALS)
  local vwcrsp::DataFrame
  local rfcrsp::DataFrame
  local ff::DataFrame

  #first, load the index portfolios
  vwcrsp = loadidxcrsp(VWCRSP_NAME, RETAINED_VWCRSP_COLS)
  rfcrsp = loadidxcrsp(RFCRSP_NAME, RETAINED_RFCRSP_COLS)
  ff = loadffdailyports()

  #interpolate monthly indices where appropriate
  rfcrsp = interpolateindex(rfcrsp, :t30ind)

  #join the indices
  N = size(fact,1)
  fact = join(fact, vwcrsp, on=:date, kind=:inner)
  fact = join(fact, rfcrsp, on=:date, kind=:inner)
  fact = join(fact, ff, on=:date, kind=:inner)

  (N ≠ size(fact,1)) && @warn("Index dates do not match factor dates.
    Factor df size reduced from $N to $(size(fact,1))")

  #compute the index returns
  indexreturns!(fact, :vwindd, :mktgross)
  indexreturns!(fact, :t30ind, :rfr)
  fact.mkt = fact.mktgross .- fact.rfr
  for s ∈
    (s->indexreturns!(fact, Symbol(s,:idx), s)).(ffportvals)
  end

  return fact
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
@with_kw struct TimeDomainMonitior{TQuantity,GeometryType} <: AbstractMonitor
	geometry::GeometryType
	quantity::TQuantity
	start_time::Float = 0.0 #  The time to start recording
	stop_time::Float = Inf # The time to stop recording
end


text:
export
    RGBDDataset,
    NYUDataset,
    get_image,
    CDN,
    CDNImage,
    CDNColorant,
    LabelledCDNImage,
    color_image,
    z_image,
    to_array,
    to_rgbdnimage,
    to_image,
    rgb2lab!,
    preprocess,
    colors,
    distances,
    normals,
    z,
    name
import Base.download
using FITSIO
using MAT
using ..RGBDSegmentation
using ColorTypes
using StaticArrays
import Base: show, +,-,*,/,getindex, setindex!,zeros, zero
import ColorTypes: comp1, comp2, comp3
import Images: accum

abstract type RGBDDataset end
import Base.getindex
function getindex(d::RGBDDataset, i::Integer)
    return get_image(d, i)
end

const CDNImage{T} = Array{T, 3} where T<:Number
@enum ImageType Color=1 Depth=2 Normal=3
function CDNImage(RGB::Array{T1, 3}, D::Array{T2, 3}, N::Array{T3, 3}) where T1 where T2 where T3
    if size(RGB,3) == 3
        RGB = permutedims(RGB./255.0, [3 1 2])
    end
    if size(D, 3) == 3
        D = permutedims(D, [3 1 2])
    end
    if size(N, 3) == 3
        N = permutedims(N, [3 1 2])
    end
    return cat(RGB, D,  N, dims=1)
end
function number_of_pixels(img::CDNImage)
    return size(img,2)*size(img,3)
end
function colors(img::CDNImage)
    return colorview(RGB, img[1:3, :, :])
end
function distances(img::CDNImage)
    return img[4:6, :, :]
end
function z(img::CDNImage)
    return img[6, :, :]
end
function normals(img::CDNImage)
    return img[7:9, :, :]
end
struct LabelledCDNImage
    image::CDNImage
    ground_truth::SegmentedImage
end
function LabelledCDNImage(image::CDNImage, labels::AbstractArray{T}) where T<:Integer
    return LabelledCDNImage(image,
                            SegmentedImage(image, labels))
end


text:
using Revise

using Statistics, Distributions, StatsBase, Random, LinearAlgebra, DelimitedFiles, Distributed
# @everywhere using Pkg
# @everywher Pkg.activate(homedir()*"/GIT/SuSiEGLMM.jl")
 using SuSiEGLMM

Seed(124)

n =100;p=10;q=5;
 L=1; B=100;
#  τ2=0.4;
#GLM
b_true=zeros(p);
b_1s=zeros(B);

res=[];Tm=zeros(B);

## random covariates

for j = 1:B
    b_true[1]= randn(1)[1]
    b_1s[j] = b_true[1]
    # b_true[1]=b_1s[j]
    X=randn(n,p)
    X₀ = randn(n,q)
    δ=randn(q)
    # g=rand(MvNormal(τ2*K))
    Y= logistic.(X*b_true+X₀*δ) .>rand(n) #generating binary outcome
    Y=convert(Vector{Float64},Y)
    # writedlm("./testdata/dataY-julia.csv",Y)
    t0=@elapsed res0= susieGLM(L, ones(p)/p,Y,X,X₀;tol=1e-4)
  # t0=@elapsed  res0= fineQTL_glm(G,Y,X;L=L,tol=1e-4)
    res=[res;res0]; Tm[j]=t0
end

println("min, median, max times for susie-glm including covariates are $(minimum(Tm)), $(median(Tm)),$(maximum(Tm)).")

b̂=[res[j].α[1]*res[j].ν[1] for j=1:B]
α̂ = [res[j].α[1] for j=1:B]
using Plots
ll=@layout[a;b];
# l2=@layout[a b;c d]
p1=scatter(b_1s,b̂,xlabel= "True effects", ylabel="Posterior estimate",label=false,title="SuSiE-GLM+random covariates")
p2=scatter(b_1s,α̂, xlabel="True effects",ylabel="PIP",label=false)
plot(p1,p2,layout=ll)

# correlated covariates with predictors
n =100;p=10;
 L=1; B=100;
#  τ2=0.4;
#GLM
b_true=zeros(p);
b_1s=zeros(B);

res1=[];Tm0=zeros(B);

for j = 1:B
    b_true[1]= randn(1)[1]
    b_1s[j] = b_true[1]
    # b_true[1]=b_1s[j]
    X=randn(n,p)
    #correlated covariates
    X₁ = rand(MvNormal([1 0.9;0.9 1]),n)
    X[:,1]=X₁[1,:]
    X₀=convert(Matrix{Float64},X₁[[2],:]')
    δ=randn(1)

    # g=rand(MvNormal(τ2*K))
    Y= logistic.(X*b_true+X₀*δ) .>rand(n) #generating binary outcome
    Y=convert(Vector{Float64},Y)
    # writedlm("./testdata/dataY-julia.csv",Y)
    t0=@elapsed res0= susieGLM(L, ones(p)/p,Y,X,X₀;tol=1e-4)
  # t0=@elapsed  res0= fineQTL_glm(G,Y,X;L=L,tol=1e-4)
    res1=[res1;res0]; Tm0[j]=t0
end

println("min, median, max times for susie-glm including covariates are $(minimum(Tm0)), $(median(Tm0)),$(maximum(Tm0)).")

b̂=[res1[j].α[1]*res1[j].ν[1] for j=1:B]
α̂ = [res1[j].α[1] for j=1:B]

p3=scatter(b_1s,b̂,xlabel= "True effects", ylabel="Posterior estimate",label=false,title="SuSiE-GLM+correlated covariates")
p4=scatter(b_1s,α̂, xlabel="True effects",ylabel="PIP",label=false)
plot(p3,p4,layout=ll)


text:
module Neighborhood
using Distances
export Euclidean, Chebyshev, Cityblock, Minkowski

include("util.jl")
include("api.jl")
include("theiler.jl")
include("bruteforce.jl")
include("kdtree.jl")
include("Testing.jl")

"Currently supported search structures"
const SSS = [BruteForce, KDTree]

end  # module Neighborhood