text:
if !any([dirname(path) == dirname(@__FILE__) for path in LOAD_PATH])
    push!(LOAD_PATH, dirname(@__FILE__))
end

ENV["LINES"] = 10

using RBT: Entity, ToyDatabase, Schema, Class, Arrow, Instance, Iso, Opt, Seq, Monetary

RBT.formatdf()

CSV_DIR = joinpath(dirname(@__FILE__), "../data/tpch")

schema = Schema(
    Class(
        :region,
        Arrow(:name, ASCIIString),
        Arrow(:comment, ASCIIString),
        Arrow(
            :nation,
            lunique=false, ltotal=false, runique=true, rtotal=true,
            inverse=:region),
        select=(:name, :comment)),
    Class(
        :nation,
        Arrow(:name, ASCIIString),
        Arrow(:region, select=:name),
        Arrow(:comment, ASCIIString),
        Arrow(
            :customer,
            lunique=false, ltotal=false, runique=true, rtotal=true,
            inverse=:nation),
        Arrow(
            :supplier,
            lunique=false, ltotal=false, runique=true, rtotal=true,
            inverse=:nation),
        select=(:name, :region, :comment)),
    Class(
        :customer,
        Arrow(:name, ASCIIString),
        Arrow(:address, ASCIIString),
        Arrow(:nation, select=:name),
        Arrow(:phone, ASCIIString),
        Arrow(:acctbal, Monetary{:USD}),
        Arrow(:mktsegment, ASCIIString),
        Arrow(:comment, ASCIIString),
        Arrow(
            :order,
            lunique=false, ltotal=false, runique=true, rtotal=true,
            inverse=:customer),
        select=(:name, :address, :nation, :phone, :acctbal, :mktsegment, :comment)),
    Class(
        :supplier,
        Arrow(:name, ASCIIString),
        Arrow(:address, ASCIIString),
        Arrow(:nation, select=:name),
        Arrow(:phone, ASCIIString),
        Arrow(:acctbal, Monetary{:USD}),
        Arrow(:comment, ASCIIString),
        Arrow(
            :partsupp,
            lunique=false, ltotal=false, runique=true, rtotal=true,
            inverse=:supplier),
        Arrow(
            :lineitem,
            lunique=false, ltotal=false, runique=true, rtotal=true,
            inverse=:supplier),
        select=(:name, :address, :nation, :phone, :acctbal, :comment)),
    Class(
        :part,
        Arrow(:name, ASCIIString),
        Arrow(:mfgr, ASCIIString),
        Arrow(:brand, ASCIIString),
        Arrow(:type_, ASCIIString),
        Arrow(:size, Int),
        Arrow(:container, ASCIIString),
        Arrow(:retailprice, Monetary{:USD}),
        Arrow(:comment, ASCIIString),
        Arrow(
            :partsupp,
            lunique=false, ltotal=false, runique=true, rtotal=true,
            inverse=:part),
        Arrow(
            :lineitem,
            lunique=false, ltotal=false, runique=true, rtotal=true,
            inverse=:part),
        select=(:name, :mfgr, :brand, :type_, :size, :container, :retailprice, :comment)),
    Class(
        :partsupp,
        Arrow(:part, select=:name),
        Arrow(:supplier, select=:name),
        Arrow(:availqty, Int),
        Arrow(:supplycost, Monetary{:USD}),
        Arrow(:comment, ASCIIString),
        Arrow(
            :lineitem,
            lunique=false, ltotal=false, runique=true, rtotal=true,
            inverse=:partsupp),
        select=(:part, :supplier, :availqty, :supplycost, :comment)),
    Class(
        :order,
        Arrow(:customer, select=:name),
        Arrow(:orderstatus, ASCIIString),
        Arrow(:totalprice, Monetary{:USD}),
        Arrow(:orderdate, Date),
        Arrow(:orderpriority, ASCIIString),
        Arrow(:clerk, ASCIIString),
        Arrow(:shippriority, Int),
        Arrow(:comment, ASCIIString),
        Arrow(
            :lineitem,
            lunique=false, ltotal=false, runique=true, rtotal=true,
            inverse=:order),
        select=(:customer, :orderstatus, :totalprice, :orderdate, :orderpriority, :clerk, :shippriority, :comment)),
    Class(
        :lineitem,
        Arrow(:order, select=:customer),
        Arrow(:part, select=:name),
        Arrow(:supplier, select=:name),
        Arrow(:partsupp, select=(:part, :supplier)),
        Arrow(:linenumber, Int),
        Arrow(:quantity, Int),
        Arrow(:extendedprice, Monetary{:USD}),
        Arrow(:discount, Float64),
        Arrow(:tax, Float64),
        Arrow(:returnflag, ASCIIString),
        Arrow(:linestatus, ASCIIString),
        Arrow(:shipdate, Date),
        Arrow(:commitdate, Date),
        Arrow(:receiptdate, Date),
        Arrow(:shipinstruct, ASCIIString),
        Arrow(:shipmode, ASCIIString),
        Arrow(:comment, ASCIIString),
        select=(
            :order, :part, :supplier, :linenumber, :quantity, :extendedprice, :discount, :tax,
            :returnflag, :linestatus, :shipdate, :commitdate, :receiptdate, :shipinstruct, :shipmode, :comment)))

R = Entity{:region}
N = Entity{:nation}
C = Entity{:customer}
S = Entity{:supplier}
P = Entity{:part}
PS = Entity{:partsupp}
O = Entity{:order}
L = Entity{:lineitem}

rkey = Dict{Int, R}()
nkey = Dict{Int, N}()
ckey = Dict{Int, C}()
skey = Dict{Int, S}()
pkey = Dict{Int, P}()
pskey = Dict{Tuple{Int,Int}, PS}()
okey = Dict{Int, O}()

rs = Vector{R}()
ns = Vector{N}()
cs = Vector{C}()
ss = Vector{S}()
ps = Vector{P}()
pss = Vector{PS}()
os = Vector{O}()
ls = Vector{L}()

r_name = Vector{Iso{ASCIIString}}()
r_comment = Vector{Iso{ASCIIString}}()
r_nation = Vector{Seq{N}}()

n_name = Vector{Iso{ASCIIString}}()
n_region = Vector{Iso{R}}()
n_comment = Vector{Iso{ASCIIString}}()
n_customer = Vector{Seq{C}}()
n_supplier = Vector{Seq{S}}()

c_name = Vector{Iso{ASCIIString}}()
c_address = Vector{Iso{ASCIIString}}()
c_nation = Vector{Iso{N}}()
c_phone = Vector{Iso{ASCIIString}}()
c_acctbal = Vector{Iso{Monetary{:USD}}}()
c_mktsegment = Vector{Iso{ASCIIString}}()
c_comment = Vector{Iso{ASCIIString}}()
c_order = Vector{Seq{O}}()

s_name = Vector{Iso{ASCIIString}}()
s_address = Vector{Iso{ASCIIString}}()
s_nation = Vector{Iso{N}}()
s_phone = Vector{Iso{ASCIIString}}()
s_acctbal = Vector{Iso{Monetary{:USD}}}()
s_comment = Vector{Iso{ASCIIString}}()
s_partsupp = Vector{Seq{PS}}()
s_lineitem = Vector{Seq{L}}()

p_name = Vector{Iso{ASCIIString}}()
p_mfgr = Vector{Iso{ASCIIString}}()
p_brand = Vector{Iso{ASCIIString}}()
p_type = Vector{Iso{ASCIIString}}()
p_size = Vector{Iso{Int}}()
p_container = Vector{Iso{ASCIIString}}()
p_retailprice = Vector{Iso{Monetary{:USD}}}()
p_comment = Vector{Iso{ASCIIString}}()
p_partsupp = Vector{Seq{PS}}()
p_lineitem = Vector{Seq{L}}()

ps_part = Vector{Iso{P}}()
ps_supplier = Vector{Iso{S}}()
ps_availqty = Vector{Iso{Int}}()
ps_supplycost = Vector{Iso{Monetary{:USD}}}()
ps_comment = Vector{Iso{ASCIIString}}()
ps_lineitem = Vector{Seq{L}}()

o_customer = Vector{Iso{C}}()
o_orderstatus = Vector{Iso{ASCIIString}}()
o_totalprice = Vector{Iso{Monetary{:USD}}}()
o_orderdate = Vector{Iso{Date}}()
o_orderpriority = Vector{Iso{ASCIIString}}()
o_clerk = Vector{Iso{ASCIIString}}()
o_shippriority = Vector{Iso{Int}}()
o_comment = Vector{Iso{ASCIIString}}()
o_lineitem = Vector{Seq{L}}()

l_order = Vector{Iso{O}}()
l_part = Vector{Iso{P}}()
l_supplier = Vector{Iso{S}}()
l_partsupp = Vector{Iso{PS}}()
l_linenumber = Vector{Iso{Int}}()
l_quantity = Vector{Iso{Int}}()
l_extendedprice = Vector{Iso{Monetary{:USD}}}()
l_discount = Vector{Iso{Float64}}()
l_tax = Vector{Iso{Float64}}()
l_returnflag = Vector{Iso{ASCIIString}}()
l_linestatus = Vector{Iso{ASCIIString}}()
l_shipdate = Vector{Iso{Date}}()
l_commitdate = Vector{Iso{Date}}()
l_receiptdate = Vector{Iso{Date}}()
l_shipinstruct = Vector{Iso{ASCIIString}}()
l_shipmode = Vector{Iso{ASCIIString}}()
l_comment = Vector{Iso{ASCIIString}}()

csv = readcsv("$CSV_DIR/region.csv", comments=false)
for k = 1:size(csv, 1)
    key, name, comment = csv[k, :]
    r = R(k)
    rkey[key] = r
    push!(rs, r)
    push!(r_name, Iso{ASCIIString}(name))
    push!(r_comment, Iso{ASCIIString}(comment))
    push!(r_nation, Seq{N}(N[]))
end

csv = readcsv("$CSV_DIR/nation.csv", comments=false)
for k = 1:size(csv, 1)
    key, name, regionkey, comment = csv[k, :]
    region = rkey[regionkey]
    n = N(k)
    nkey[key] = n
    push!(ns, n)
    push!(n_name, Iso{ASCIIString}(name))
    push!(n_region, Iso{R}(region))
    push!(n_comment, Iso{ASCIIString}(comment))
    push!(n_customer, Seq{C}(C[]))
    push!(n_supplier, Seq{S}(S[]))
    push!(r_nation[region.id].data, n)
end

csv = readcsv("$CSV_DIR/customer.csv", comments=false)
for k = 1:size(csv, 1)
    key, name, address, nationkey, phone, acctbal, mktsegment, comment = csv[k, :]
    nation = nkey[nationkey]
    c = C(k)
    ckey[key] = c
    push!(cs, c)
    push!(c_name, Iso{ASCIIString}(name))
    push!(c_address, Iso{ASCIIString}(address))
    push!(c_nation, Iso{N}(nation))
    push!(c_phone, Iso{ASCIIString}(phone))
    push!(c_acctbal, Iso{Monetary{:USD}}(acctbal))
    push!(c_mktsegment, Iso{ASCIIString}(mktsegment))
    push!(c_comment, Iso{ASCIIString}(comment))
    push!(c_order, Seq{O}(O[]))
    push!(n_customer[nation.id].data, c)
end

csv = readcsv("$CSV_DIR/supplier.csv", comments=false)
for k = 1:size(csv, 1)
    key, name, address, nationkey, phone, acctbal, comment = csv[k, :]
    nation = nkey[nationkey]
    s = S(k)
    skey[key] = s
    push!(ss, s)
    push!(s_name, Iso{ASCIIString}(name))
    push!(s_address, Iso{ASCIIString}(address))
    push!(s_nation, Iso{N}(nation))
    push!(s_phone, Iso{ASCIIString}(phone))
    push!(s_acctbal, Iso{Monetary{:USD}}(acctbal))
    push!(s_comment, Iso{ASCIIString}(comment))
    push!(s_partsupp, Seq{PS}(PS[]))
    push!(s_lineitem, Seq{L}(L[]))
    push!(n_supplier[nation.id].data, s)
end

csv = readcsv("$CSV_DIR/part.csv", comments=false)
for k = 1:size(csv, 1)
    key, name, mfgr, brand, type_, size, container, retailprice, comment = csv[k, :]
    p = P(k)
    pkey[key] = p
    push!(ps, p)
    push!(p_name, Iso{ASCIIString}(name))
    push!(p_mfgr, Iso{ASCIIString}(mfgr))
    push!(p_brand, Iso{ASCIIString}(brand))
    push!(p_type, Iso{ASCIIString}(type_))
    push!(p_size, Iso{Int}(size))
    push!(p_container, Iso{ASCIIString}(container))
    push!(p_retailprice, Iso{Monetary{:USD}}(retailprice))
    push!(p_comment, Iso{ASCIIString}(comment))
    push!(p_partsupp, Seq{PS}(PS[]))
    push!(p_lineitem, Seq{L}(L[]))
end

csv = readcsv("$CSV_DIR/partsupp.csv", comments=false)
for k = 1:size(csv, 1)
    partkey, suppkey, availqty, supplycost, comment = csv[k, :]
    ps_ = PS(k)
    part = pkey[partkey]
    supplier = skey[suppkey]
    pskey[partkey, suppkey] = ps_
    push!(pss, ps_)
    push!(ps_part, Iso{P}(part))
    push!(ps_supplier, Iso{S}(supplier))
    push!(ps_availqty, Iso{Int}(availqty))
    push!(ps_supplycost, Iso{Monetary{:USD}}(supplycost))
    push!(ps_comment, Iso{ASCIIString}(comment))
    push!(ps_lineitem, Seq{L}(L[]))
    push!(s_partsupp[supplier.id].data, ps_)
    push!(p_partsupp[part.id].data, ps_)
end

csv = readcsv("$CSV_DIR/orders.csv", comments=false)
for k = 1:size(csv, 1)
    orderkey, custkey, orderstatus, totalprice, orderdate, orderpriority, clerk, shippriority, comment = csv[k, :]
    o = O(k)
    customer = ckey[custkey]
    okey[orderkey] = o
    push!(os, o)
    push!(o_customer, Iso{C}(customer))
    push!(o_orderstatus, Iso{ASCIIString}(orderstatus))
    push!(o_totalprice, Iso{Monetary{:USD}}(totalprice))
    push!(o_orderdate, Iso{Date}(Date(orderdate)))
    push!(o_orderpriority, Iso{ASCIIString}(orderpriority))
    push!(o_clerk, Iso{ASCIIString}(clerk))
    push!(o_shippriority, Iso{Int}(shippriority))
    push!(o_comment, Iso{ASCIIString}(comment))
    push!(o_lineitem, Seq{L}(L[]))
    push!(c_order[customer.id].data, o)
end

csv = readcsv("$CSV_DIR/lineitem.csv", comments=false)
for k = 1:size(csv, 1)
    orderkey, partkey, suppkey, linenumber, quantity, extendedprice, discount, tax,
    returnflag, linestatus, shipdate, commitdate, receiptdate,
    shipinstruct, shipmode, comment = csv[k, :]
    l = L(k)
    order = okey[orderkey]
    part = pkey[partkey]
    supplier = skey[suppkey]
    partsupp = pskey[partkey, suppkey]
    push!(ls, l)
    push!(l_order, Iso{O}(order))
    push!(l_part, Iso{P}(part))
    push!(l_supplier, Iso{S}(supplier))
    push!(l_partsupp, Iso{PS}(partsupp))
    push!(l_linenumber, Iso{Int}(linenumber))
    push!(l_quantity, Iso{Int}(quantity))
    push!(l_extendedprice, Iso{Monetary{:USD}}(extendedprice))
    push!(l_discount, Iso{Float64}(discount))
    push!(l_tax, Iso{Float64}(tax))
    push!(l_returnflag, Iso{ASCIIString}(returnflag))
    push!(l_linestatus, Iso{ASCIIString}(linestatus))
    push!(l_shipdate, Iso{Date}(Date(shipdate)))
    push!(l_commitdate, Iso{Date}(Date(commitdate)))
    push!(l_receiptdate, Iso{Date}(Date(receiptdate)))
    push!(l_shipinstruct, Iso{ASCIIString}(shipinstruct))
    push!(l_shipmode, Iso{ASCIIString}(shipmode))
    push!(l_comment, Iso{ASCIIString}(comment))
    push!(o_lineitem[order.id].data, l)
    push!(p_lineitem[part.id].data, l)
    push!(s_lineitem[supplier.id].data, l)
    push!(ps_lineitem[partsupp.id].data, l)
end

instance = Instance(
    Dict(
        :region => rs,
        :nation => ns,
        :customer => cs,
        :supplier => ss,
        :part => ps,
        :partsupp => pss,
        :order => os,
        :lineitem => ls),
    Dict(
        (:region, :name) => r_name,
        (:region, :comment) => r_comment,
        (:region, :nation) => r_nation,
        (:nation, :name) => n_name,
        (:nation, :region) => n_region,
        (:nation, :comment) => n_comment,
        (:nation, :customer) => n_customer,
        (:nation, :supplier) => n_supplier,
        (:customer, :name) => c_name,
        (:customer, :address) => c_address,
        (:customer, :nation) => c_nation,
        (:customer, :phone) => c_phone,
        (:customer, :acctbal) => c_acctbal,
        (:customer, :mktsegment) => c_mktsegment,
        (:customer, :comment) => c_comment,
        (:customer, :order) => c_order,
        (:supplier, :name) => s_name,
        (:supplier, :address) => s_address,
        (:supplier, :nation) => s_nation,
        (:supplier, :phone) => s_phone,
        (:supplier, :acctbal) => s_acctbal,
        (:supplier, :comment) => s_comment,
        (:supplier, :partsupp) => s_partsupp,
        (:supplier, :lineitem) => s_lineitem,
        (:part, :name) => p_name,
        (:part, :mfgr) => p_mfgr,
        (:part, :brand) => p_brand,
        (:part, :type_) => p_type,
        (:part, :size) => p_size,
        (:part, :container) => p_container,
        (:part, :retailprice) => p_retailprice,
        (:part, :comment) => p_comment,
        (:part, :partsupp) => p_partsupp,
        (:part, :lineitem) => p_lineitem,
        (:partsupp, :part) => ps_part,
        (:partsupp, :supplier) => ps_supplier,
        (:partsupp, :availqty) => ps_availqty,
        (:partsupp, :supplycost) => ps_supplycost,
        (:partsupp, :comment) => ps_comment,
        (:partsupp, :lineitem) => ps_lineitem,
        (:order, :customer) => o_customer,
        (:order, :orderstatus) => o_orderstatus,
        (:order, :totalprice) => o_totalprice,
        (:order, :orderdate) => o_orderdate,
        (:order, :orderpriority) => o_orderpriority,
        (:order, :clerk) => o_clerk,
        (:order, :shippriority) => o_shippriority,
        (:order, :comment) => o_comment,
        (:order, :lineitem) => o_lineitem,
        (:lineitem, :order) => l_order,
        (:lineitem, :part) => l_part,
        (:lineitem, :supplier) => l_supplier,
        (:lineitem, :partsupp) => l_partsupp,
        (:lineitem, :linenumber) => l_linenumber,
        (:lineitem, :quantity) => l_quantity,
        (:lineitem, :extendedprice) => l_extendedprice,
        (:lineitem, :discount) => l_discount,
        (:lineitem, :tax) => l_tax,
        (:lineitem, :returnflag) => l_returnflag,
        (:lineitem, :linestatus) => l_linestatus,
        (:lineitem, :shipdate) => l_shipdate,
        (:lineitem, :commitdate) => l_commitdate,
        (:lineitem, :receiptdate) => l_receiptdate,
        (:lineitem, :shipinstruct) => l_shipinstruct,
        (:lineitem, :shipmode) => l_shipmode,
        (:lineitem, :comment) => l_comment))

tpch = ToyDatabase(schema, instance)

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
using BinaryProvider # requires BinaryProvider 0.5.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = LibraryProduct[]

dependencies = [
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.4/build_Zlib.v1.2.11.jl",
    "https://github.com/JuliaWeb/MbedTLSBuilder/releases/download/v0.20.0/build_MbedTLS.v2.6.1.jl",
    "https://github.com/JuliaWeb/LibCURLBuilder/releases/download/v0.5.1/build_LibCURL.v7.64.1.jl",
]

for dependency in dependencies
    file = joinpath(@__DIR__, basename(dependency))

    # Use the BinaryProvider download function as it supports proxies.
    # See: https://github.com/JuliaWeb/LibCURL.jl/issues/69
    isfile(file) || BinaryProvider.download(dependency, file)

    # Build the dependencies
    # Note: It is a bit faster to run the build in an anonymous module instead of starting a new
    # julia process
    m = @eval module Anon end
    m.include(file)
    append!(products, m.products)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)

ext: jl


text:
module TOMLX

using Dates

module Internals
    const HAVE_DATES = true
    include("parser.jl")
end # module

"""
    Parser()

Constructor for a TOML `Parser`. After creation the
function [`TOML.reinit!`](@ref) is used to initialize
the parser and then `TOML.parse` is called to parse
the data.
Note that in most cases one does not need to explicitly create
a `Parser` but instead one directly use
use [`parsefile`](@ref) or [`parsestring`](@ref).
"""
const Parser = Internals.Parser

"""
    parsefile(f::AbstractString)
    parsefile(p::Parser, f::AbstractString)

Parses a file `f` and returns the resulting
table (dictionary). Throws a [`ParserError`](@ref)
upon failure.

See also [`tryparsefile`](@ref)
"""
parsefile(f::AbstractString) =
    Internals.parse(Parser(read(f, String); filepath=abspath(f)); raise=true)
parsefile(p::Parser, f::AbstractString) =
    Internals.parse(Internals.reinit!(p, read(f, String); filepath=abspath(f)); raise=true)

"""
    tryparsefile(f::AbstractString)
    tryparsefile(p::Parser, f::AbstractString)

Parses a file `f` and returns the resulting
table (dictionary). Returns a [`ParserError`](@ref)
upon failure.

See also [`parsefile`](@ref)
"""
tryparsefile(f::AbstractString) =
    Internals.parse(Parser(read(f, String); filepath=abspath(f)); raise=false)
tryparsefile(p::Parser, f::AbstractString) =
    Internals.parse(Internals.reinit!(p, read(f, String); filepath=abspath(f)); raise=false)

"""
    parsestring(str::AbstractString)
    parsestring(p::Parser, str::AbstractString)

Parses a string `str` and returns the resulting
table (dictionary). Returns a [`ParserError`](@ref)
upon failure.

See also [`tryparsestring`](@ref)
"""
parsestring(str::AbstractString) =
    Internals.parse(Parser(String(str)); raise=true)
parsestring(p::Parser, str::AbstractString) =
    Internals.parse(Internals.reinit!(p, String(str)); raise=true)

"""
    tryparsestring(str::AbstractString)
    tryparsestring(p::Parser, str::AbstractString)

Parses a string `str` and returns the resulting
table (dictionary). Returns a [`ParserError`](@ref)
upon failure.
1
See also [`parsestring`](@ref)
"""
tryparsestring(str::AbstractString) =
    Internals.parse(Parser(String(str)); raise=false)
tryparsestring(p::Parser, str::AbstractString) =
    Internals.parse(Internals.reinit!(p, String(str)); raise=false)

"""
    ParserError
Type that is returned from [`tryparsestring`](@ref) and
[`tryparsefile`](@ref) when parsing fails. It contains the following
fields:
- `pos`, the position in the string when the error happened
- `table`, the result that so far was successfully parsed
- `type`, an error type, different for different type of errors
"""
const ParserError = Internals.ParserError

end


text:
module TestScheduler

using TapirSchedulers:
    TapirSchedulers, WORK_STEALING_SCHEDULER, DEPTH_FIRST_SCHEDULER, rollback_priority
using Test

test_simple_spawn_work_stealing() = check_simple_spawn(WORK_STEALING_SCHEDULER[])

function test_simple_spawn_depth_first()
    rollback_priority() do
        check_simple_spawn(DEPTH_FIRST_SCHEDULER[])
    end
end

function check_simple_spawn(sch)
    w = TapirSchedulers.spawn!(sch) do
        1 + 1
    end
    @test fetch(w) == 2
end

function fib(N, sch)
    if N <= 1
        return N
    end
    t = TapirSchedulers.spawn!(sch) do
        fib(N - 2, sch)
    end
    y = fib(N - 1, sch)
    return fetch(t)::Int + y
end

fib_ws(N) = fib(N, WORK_STEALING_SCHEDULER[])
function fib_df(N)
    rollback_priority() do
        fib(N, DEPTH_FIRST_SCHEDULER[])
    end
end

test_fib_work_stealing() = check_fib(fib_ws)
test_fib_depth_first() = check_fib(fib_df)

function check_fib(fib)
    @test fib(1) == 1
    @test fib(2) == 1
    @test fib(3) == 2
    @test fib(4) == 3
    @test fib(5) == 5
    @test fib(6) == 8
    @test fib(10) == 55
end

end  # module


text:
"""
    Shuffle

Support for a number of deterministic and random shuffling algorithms. Provides functions
[`shuffle`](@ref), [`shuffle!`](@ref), [`nshuffle`](@ref) and [`nshuffle!`](@ref)
as well as the following shuffling algorithms:
- [faro (or weave) shuffle](https://luapulu.github.io/Shuffle.jl/stable/reference/#Shuffle.Faro),
- a [cut](https://luapulu.github.io/Shuffle.jl/stable/reference/#Shuffle.Cut),
- [random shuffle](https://luapulu.github.io/Shuffle.jl/stable/reference/#Shuffle.RandomShuffle)
    (uses [`Random.shuffle`](https://docs.julialang.org/en/v1/stdlib/Random/#Random.shuffle)) and
- [Gilbert-Shannon-Reeds model](https://luapulu.github.io/Shuffle.jl/stable/reference/#Shuffle.GilbertShannonReeds).
"""
module Shuffle

import Base: show
import Random # used so Random.shuffle! can be seperated from shuffle!

using Random: AbstractRNG, default_rng
using Base: copymutable

export AbstractShuffle, AbstractDeterministicShuffle, AbstractRandomShuffle
export Faro, Weave, Cut, RandomShuffle, GilbertShannonReeds
export infaro, outfaro, inweave, outweave, randshuffle, gsrshuffle
export shuffle, shuffle!, nshuffle, nshuffle!

## Abstract Types ##

"""
    AbstractShuffle

Supertype for [`AbstractDeterministicShuffle`](@ref) and [`AbstractRandomShuffle`](@ref).

# Implementation

New shuffling algorithm types should be sub types of either
[`AbstractDeterministicShuffle`](@ref) or [`AbstractRandomShuffle`](@ref). If an algorithm
can be defined in-place, only [`shuffle!`](@ref) needs to be extended. [`shuffle`](@ref),
[`nshuffle`](@ref) and [`nshuffle!`](@ref) will make a copy / repeatedly call
[`shuffle!`](@ref) as needed. If the algorithm can not be defined in-place, define
[`shuffle`](@ref) and [`nshuffle!`](@ref). [`nshuffle`](@ref) will make a copy and call
[`nshuffle!`](@ref).
"""
abstract type AbstractShuffle end

"""
    AbstractDeterministicShuffle <: AbstractShuffle

Supertype for fully deterministic shuffle algorithms such as [`Faro`](@ref) and
[`Cut`](@ref).

# Implementation

New deterministic shuffling algorithms should implement either [`shuffle`](@ref) or
[`shuffle!`](@ref), as appropriate. The new shuffle method should take a collection to be
shuffled as the first argument and an instance of the new shuffling algorithm type as the
second argument.

For example:
```julia
struct MyShuffle <: AbstractDeterministicShuffle
    parameter
end

function shuffle(c::AbstractArray, s::MyShuffle, prealloc_out=similar(c))
    # Do the shuffling

    return prealloc_out
end
```

Pre-allocating the output makes it easy to also write an efficient [`nshuffle!`](@ref)
function for algorithms that cannot be written in-place.

See also: [`AbstractShuffle`](@ref)
"""
abstract type AbstractDeterministicShuffle <: AbstractShuffle end

"""
    AbstractRandomShuffle <: AbstractShuffle

Supertype for partly or fully random shuffle algorithms such as [`RandomShuffle`](@ref) and
[`GilbertShannonReeds`](@ref).

# Implementation

New random shuffling algorithms should implement either [`shuffle`](@ref) or
[`shuffle!`](@ref), as appropriate. The new shuffle method should take a
[`Random.AbstractRNG`](https://docs.julialang.org/en/v1/stdlib/Random/#Random.AbstractRNG)
as the first argument, a collection to be shuffled as the second argument and an instance of
the new shuffling algorithm type as the third argument.

For example:
```julia
struct MyShuffle <: AbstractRandomShuffle end

function shuffle!(r::AbstractRNG, c, s::MyShuffle)
    # Do the shuffling, passing r to any random number generating functions

    return c
end
```

See also: [`AbstractShuffle`](@ref)
"""
abstract type AbstractRandomShuffle <: AbstractShuffle end

## Include files ##

include("DeterministicShuffles.jl")
include("RandomShuffles.jl")

## Defaults ##

mutable struct Defaults
    shuffle::AbstractShuffle
end

"""
    DEFAULTS

binding to a mutable struct containing the default shuffling algorithm.
[`RandomShuffle`](@ref) is set initially.

# Examples
```jldoctest
julia> Shuffle.DEFAULTS.shuffle
RandomShuffle()

julia> mt1 = MersenneTwister(1234); mt2 = MersenneTwister(1234);

julia> shuffle(mt1, collect(1:100)) == shuffle(mt2, collect(1:100), RandomShuffle())
true

julia> Shuffle.DEFAULTS.shuffle = infaro; Shuffle.DEFAULTS.shuffle
Faro{:in}()

julia> shuffle([1, 2, 3, 4, 5, 6, 7, 8])
8-element Vector{Int64}:
 5
 1
 6
 2
 7
 3
 8
 4

julia> Shuffle.DEFAULTS.shuffle = RandomShuffle();
```
"""
const DEFAULTS = Defaults(RandomShuffle())

## Shuffle! and shuffle ##

"""
    shuffle!([rng=GLOBAL_RNG,], c, s::AbstractShuffle=RandomShuffle())

shuffle collection `c` in-place using shuffling algorithm `s`. A random-number generator `rng`
may be supplied for random shuffling algorithms.

# Examples
```jldoctest
julia> mt = MersenneTwister(1234);

julia> a = collect(1:6);

julia> shuffle!(mt, a); a
6-element Vector{Int64}:
 2
 1
 3
 6
 4
 5

julia> shuffle!(a, Faro{:out}()); a
6-element Vector{Int64}:
 2
 6
 1
 4
 3
 5
```
"""
shuffle!(c) = shuffle!(c, DEFAULTS.shuffle)
shuffle!(r::AbstractRNG, c) = shuffle!(r, c, DEFAULTS.shuffle)
shuffle!(c, s::AbstractRandomShuffle) = shuffle!(default_rng(), c, s)

"""
    shuffle([rng=GLOBAL_RNG,], c, s::AbstractShuffle=RandomShuffle())

shuffle collection `c` using shuffling algorithm `s`. To shuffle `c` in-place,
see [`shuffle!`](@ref). A random-number generator `rng` may be supplied for random shuffling
algorithms.

# Examples
```jldoctest
julia> mt = MersenneTwister(1234);

julia> shuffle(mt, [1, 2, 3, 4, 5, 6, 7, 8], GilbertShannonReeds())
8-element Vector{Int64}:
 6
 7
 1
 2
 8
 3
 4
 5

julia> shuffle([6, 5, 4, 3, 2, 1], Faro{:in}())
6-element Vector{Int64}:
 3
 6
 2
 5
 1
 4
```
"""
shuffle(c, s::AbstractShuffle=DEFAULTS.shuffle) = shuffle!(copymutable(c), s)
shuffle(r::AbstractRNG, c, s::AbstractRandomShuffle=DEFAULTS.shuffle) = shuffle!(r, copymutable(c), s)
shuffle(c, s::AbstractRandomShuffle) = shuffle(default_rng(), c, s)

## nshuffle! and nshuffle ##

"""
    nshuffle!([rng=GLOBAL_RNG,], c, n::Integer, s::AbstractShuffle=RandomShuffle())

shuffle collection `c` in-place `n` times using shuffling algorithm `s`. A random-number
generator `rng` may be supplied for random shuffling algorithms.

# Examples
```jldoctest
julia> mt = MersenneTwister(1234);

julia> a = collect(1:7);

julia> nshuffle!(mt, a, 3, GilbertShannonReeds()); a
7-element Vector{Int64}:
 5
 6
 1
 7
 2
 3
 4
```
"""
nshuffle!(c, n::Integer) = nshuffle!(c, n, DEFAULTS.shuffle)

function nshuffle!(r::AbstractRNG, c, n::Integer, s::AbstractRandomShuffle=DEFAULTS.shuffle)
    for i in 1:n
        shuffle!(r, c, s)
    end
    return c
end

nshuffle!(c, n::Integer, s::AbstractRandomShuffle) = nshuffle!(default_rng(), c, n, s)

function nshuffle!(c, n::Integer, s::AbstractDeterministicShuffle)
    for i in 1:n
        shuffle!(c, s)
    end
    return c
end

"""
    nshuffle([rng=GLOBAL_RNG,], c, n::Integer, s::AbstractShuffle=RandomShuffle())

shuffle collection `c` `n` times using shuffling algorithm `s`. A random-number generator
`rng` may be supplied for random shuffling algorithms.
To shuffle `c` in-place see [`nshuffle!`](@ref)

# Examples
```jldoctest
julia> nshuffle(collect(1:8), 3, Faro{:in}())
8-element Vector{Int64}:
 8
 7
 6
 5
 4
 3
 2
 1
```
"""
nshuffle(c, n::Integer, s::AbstractShuffle=DEFAULTS.shuffle) = nshuffle!(copymutable(c), n, s)
nshuffle(c, n::Integer, s::AbstractRandomShuffle) = nshuffle(default_rng(), c, n, s)
nshuffle(r::AbstractRNG, c, n::Integer, s::AbstractRandomShuffle=DEFAULTS.shuffle) =
    nshuffle!(r, copymutable(c), n, s)

end # module


text:
# include the strategy for this language -
include("./strategy/JuliaStrategy.jl")

"""
    generate(julia_model_object::VLJuliaModelObject;
    intermediate_representation_dictionary::Union{Nothing,Dict{String,Any}} = nothing,
    logger::Union{Nothing,SimpleLogger} = nothing)

Generate model code from a Julia data object and write generated code to the user-defined output directory.

Input arguments:
`julia_model_object::VLJuliaModelObject` - abstract data object for generating model code in Julia language.
`intermediate_representation_dictionary::Union{Nothing,Dict{String,Any}}` - an intermediate data dictionary for building model code (optional).
`logger::Union{Nothing,SimpleLogger}` - logger to log information during the code generation process (optional).

Output arguments:
None
"""
function generate(julia_model_object::VLJuliaModelObject;
    intermediate_representation_dictionary::Union{Nothing,Dict{String,Any}} = nothing,
    logger::Union{Nothing,SimpleLogger} = nothing)

    # initialize -
    src_component_set = Set{NamedTuple}()
    config_component_set = Set{NamedTuple}()
    root_component_set = Set{NamedTuple}()

    try

        # check: do we have an intermediate representation?
        ir_dictionary = intermediate_representation_dictionary
        if (isnothing(intermediate_representation_dictionary) == true)

            # parse the vff document -
            result = parse_vff_model_document(julia_model_object)
            if (isa(result.value,Exception) == true)
                rethrow(result.value) # TODO: what is the diff between re-throw -vs- throw?
            end
            ir_dictionary = result.value
        end

        # Step 1: Copy the "distribution" files to thier location -
        path_to_output_dir = julia_model_object.path_to_output_dir

        # Transfer distrubtion jl files to the output -> these files are shared between model types
        transfer_distribution_files("$(_PATH_TO_SRC)/distribution/julia/src", "$(path_to_output_dir)/src",".jl")

        # Transfer root jl files -> these files are moved to the root dir of the file
        transfer_distribution_files("$(_PATH_TO_SRC)/distribution/julia/root", "$(path_to_output_dir)",".jl")

        # Transfer root toml files -> these files are moved to the root dir of the model
        transfer_distribution_files("$(_PATH_TO_SRC)/distribution/julia/root", "$(path_to_output_dir)",".toml")

        # start building the custom program components -
        # Build the Data.jl data dictionary -
        program_component = build_data_dictionary_program_component(ir_dictionary)
        if (isa(program_component.value, Exception) == true)
            throw(program_component.value)
        end
        data_dictionary_component = program_component.value
        push!(src_component_set, data_dictionary_component)

        # Build the Kinetics.jl file which holds all the rates -
        program_component = build_kinetics_program_component(ir_dictionary)
        if (isa(program_component.value, Exception) == true)
            throw(program_component.value)
        end
        kinetics_file_component = program_component.value
        push!(src_component_set, kinetics_file_component)

        # Build the Control.jl file which holds the control values for each rate -
        program_component = build_control_program_component(ir_dictionary)
        if (isa(program_component.value, Exception) == true)
            throw(program_component.value)
        end
        control_file_component = program_component.value
        push!(src_component_set, control_file_component)

        # Generate the stoichiometric_matrix -
        stm_generation_result = generate_stoichiometric_matrix(ir_dictionary)
        if (isa(stm_generation_result.value,Exception) == true)
            throw(stm_generation_result.value)
        end
        stoichiometric_matrix = stm_generation_result.value
        stm_program_component = (matrix=stoichiometric_matrix, filename="Network.dat", component_type=:matrix)
        push!(src_component_set, stm_program_component)

        # dump src and config components to disk -
        _output_path_to_src_distribution_files = joinpath(path_to_output_dir,"src")
        write_program_components_to_disk(_output_path_to_src_distribution_files, src_component_set)

        # let the user know the model has been generated ...
        @info "Model code generation has completed. Please check: $(path_to_output_dir) for your model files."

    catch error

        # let the user know that something went wrong, and prompt them to the log file -
        @info "Model code generation failed before completion. If logging is enabled, please check the log file."

        # ok, so if we get an error, log it (if we have a logger), then explode -
        if (isnothing(logger) == false)
            with_logger(logger) do
                # generate a back trace -
                @error "Ooops! we encountered an error" exception=(error, catch_backtrace())
            end
        end
    end
end

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
using StanOptimize, Test

@testset "Bernoulli optimize" begin
  include(joinpath(@__DIR__, "../examples/Bernoulli/bernoulli.jl"))
  @test optim1["theta"][end] ≈ 0.3 atol=0.1
  @test optim2["theta"][end] ≈ 0.3 atol=0.1
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CSDP_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CSDP")
JLLWrappers.@generate_main_file("CSDP", UUID("9ce75daa-2788-5e2c-ba1d-cf8c48367b27"))
end  # module CSDP_jll


text:
# Test data
using irtfun
using RCall, DataFrames

Science = R"""
library(mirt)
Science
"""
Science = convert(Matrix, Science)

# Graded Response model
using Plots
pars = graded(a = [1.0], d = [2, 0, -1.5])
[irf.(Ref(pars), -4:0.1:4, k, 1) for k in 0:3] |> plot

# 3PLM (guessing) model
pars = guessing(a = [1.0], d = [0.0], c = [0.2])
[irf.(Ref(pars), -4:0.1:4, k, 1) for k in 0:1] |> plot

# Gauss Hermite
using Distributions, FastGaussQuadrature, LinearAlgebra
n, w = gausshermite(20)
# Rescaling
dot(pdf.(Normal(0, 1), n./√2), w.*exp.(n.^2)) / √2

# Cholesky decomposition
using LinearAlgebra
A = [1.0 0.6; 0.6 1.0]
B = cholesky(A)
# B.L'B.L
B.U'B.U

D = combine(groupby([Science DataFrame(__group = fill(1, nrow(Science)))], All()), nrow => :__count)
U = D[:, Not([:__count, :__group])]
F = D[:, :__count]
G = D[:, :__group]
U = mapcols(c -> rescore(c), U) # rescoreは，オプションで選択可能にする。エラー検知もできるとよい。
pind, G = takepersonindex(G)
nind = Dict(zip(axes(U, 1), F))
iind = [takeitemindex(@view U[i,:]) for i in axes(U, 1)]
giind = [findall(map(c -> !all(ismissing(c)), eachcol(U[pind[g], :]))) for g in axes(pind, 1)]
P, I = size(U)
nnode = 5
lnL = initializematrix(P, nnode)
post = initializematrix(P, nnode)
Q = explanatoryMIRT(1, 4)
pars = [init(graded(), Q[i,:], Science[:, i]) for i in 1:I]
[setgroupindex!(pars[j], U[:,j], G) for j in 1:I]
gh = [mcov(size(Q, 2)) for _ in 1:P]
gh = [initializenode(size(Q, 2), nnode) for _ in 1:length(unique(G))]
bgh = initializenode(size(Q, 2), nnode; dentype = :GH)
bgh = initializenode(size(Q, 2), nnode; dentype = :DN)
pmoments = [Moments(size(Q, 2)) for _ in unique(G)]

# EM cycle
BAEMcycle(lnL, post, pind, nind, iind, gh, bgh, U, pars, Q; Mstep = Mstep())
@profview BAEMcycle(lnL, post, pind, nind, iind, gh, bgh, U, pars, Q; Mstep = Mstep())

# Estep
# DN
BAEstep!(lnL, post, pind, iind, gh, U, pars, Q)
@code_warntype BAEstep!(lnL, post, pind, iind, gh, U, pars, Q)
# GH
BAEstep!(lnL, post, pind, iind, gh, bgh, U, pars, Q)
@code_warntype BAEstep!(lnL, post, pind, iind, gh, bgh, U, pars, Q)

# using Plots; plot(bgh.n, post')

N, r = expectedcount(post, pind, nind, iind, pars, U)
@code_warntype expectedcount(post, pind, nind, iind, pars, U)

# Adapt the quadrature points
moments = [_eap(post[p,:], bgh) for p in axes(post, 1)]
[adapt!(gh[p], moments[p]) for p in axes(gh, 1)]

# Estimate the moments of the population
@code_warntype updateDN!(pind, nind, post, gh, pmoments)
updateAGH!(pind, nind, moments, pmoments, 1)

# Mstep

# ForwardDiff version
using ForwardDiff
# DN
@time test = ell([1.2, 0, -1.0, -2.0], pars[1], bgh, r[1], Q[1,:]) # -787くらい
@code_warntype ell([1.2, 0, -1.0, -2.0], pars[1], bgh, r[1], Q[1,:])
g = (j::Int64) -> ForwardDiff.gradient(x′ -> ell(x′, pars[j], bgh, r[1], Q[1,:]), take(pars[j], Q[j,:]))
H = (j::Int64) -> ForwardDiff.hessian(x′ -> ell(x′, pars[j], bgh, r[1], Q[1,:]), take(pars[j], Q[j,:]))
@code_warntype g(1)
@code_warntype Newton(pars, g, H, 1)
@code_warntype H(1)
# GH
@time test = ell([1.2, 0, -1.0, -2.0], pars[1], gh, bgh, post, pind, nind, Q[1,:], U[:, 1]) # -787くらい
@code_warntype ell([1.2, 0, -1.0, -2.0], pars[1], gh, bgh, post, pind, nind, Q[1,:], U[:, 1])
g = (j::Int64) -> ForwardDiff.gradient(x′ -> ell(x′, pars[j], gh, bgh, post, pind, nind, Q[j,:], @view U[:,j]), take(pars[j], Q[j,:]))
H = (j::Int64) -> ForwardDiff.hessian(x′ -> ell(x′, pars[j], gh, bgh, post, pind, nind, Q[j,:], @view U[:,j]), take(pars[j], Q[j,:]))

for j in 1:4
    @time @fastmath Newton(pars, g, H, j, Q[j,:])
end

oldpars = irtfun.copy.(pars)
@code_warntype BAEMcycle(lnL, post, pind, nind, iind, gh, bgh, U, pars, Q)
checkconvergence(pars, oldpars)

# Adaptive quadrature
moments = [_eap(post[p,:], gh[p]) for p in axes(post, 1)]
[adapt!(gh[p], bgh, moments[p]) for p in axes(gh, 1)]
update!(pind, nind, moments, pmoments, size(Q, 2))
[rescale!(pars[i], pmoments[1].μ, pmoments[1].Σ, Q[i,:]) for i in axes(pars, 1)]

@time estimate(Science, 1, graded(), E = Estep(quadpts = 11, N = 10))
@time estimate(Science, 1, graded(), E = Estep(quadpts = 61))
@profview estimate(Science, 1, graded(), E = Estep(N = 10, quadpts = 121))
@code_warntype estimate(Science, 1, graded(), E = Estep(N = 2))

@rlibrary mirt
@time mirt(Science, 1)

# Multiple group
# Test data
using irtfun
using RCall, DataFrames

mgdata = R"""
library(mirt)
set.seed(12345)
a <- matrix(abs(rnorm(15,1,.3)), ncol=1)
d <- matrix(rnorm(15,0,.7),ncol=1)
itemtype <- rep('2PL', nrow(a))
N <- 1000
dataset1 <- simdata(a, d, N, itemtype)
dataset2 <- simdata(a, d, N, itemtype, mu = .1, sigma = matrix(1.5))
dat <- rbind(dataset1, dataset2)
group <- c(rep('D1', N), rep('D2', N))
list(dat, group)
"""
R"""
models <- 'F1 = 1-15'
mod_configural <- multipleGroup(dat, models, group = group)
"""
data = DataFrame(convert(Matrix, mgdata[1]))
group = convert(Vector, mgdata[2])

D = combine(groupby([data DataFrame(__group = group)], All()), nrow => :__count)
U = D[:, Not([:__count, :__group])]
F = D[:, :__count]
G = D[:, :__group]
U = mapcols(c -> rescore(c), U) # rescoreは，オプションで選択可能にする。エラー検知もできるとよい。
pind, G = takepersonindex(G)
nind = Dict(zip(axes(U, 1), F))
iind = [takeitemindex(@view U[i,:]) for i in axes(U, 1)]
P, I = size(U)
lnL = initializematrix(P, 61; use_static = true)
[initializematrix(1, 61; use_static = true) for _ in 1:2000]
post = initializematrix(P, 61; use_static = true)
Q = fill(true, I)
pars = [init(graded(), Q[i,:], data[:, i]) for i in 1:I]
[setgroupindex!(pars[j], U[:,j], G) for j in 1:I]
# gh = [initializenode(size(Q, 2), 21) for _ in 1:P]
gh = [initializenode(size(Q, 2), 61; use_static = true) for _ in 1:length(unique(G))]
bgh = initializenode(size(Q, 2), 61; use_static = true)
pmoments = [Moments(size(Q, 2)) for _ in unique(G)]


text:
# want:
# - truew -- true width of Rchannel in m
#
# have
# - ppm: pix per m at surface
#

"""
    pixel2meter(pixels, exi::ExpImgs)

This gives the pixel to meter conversion for the plane in which the
R-channel lies.
"""
function pixel2meter(pixels, exi::ExpImgs)
    r1 = pixels/pixels_per_m(exi)
    r1_to_R(r1, exi)
end

pixels_per_m(exi::ExpImgs) = exi.p2m.ppm/exi.thin_num

"""
    r1_to_R(r1, p2m::Pixel2Meter)
    r1_to_R(r1, exi::ExpImgs)

Transforms a distance as measured at the surface of the block into a
distance inside the ice block.

Notes:
- see Yuri's thesis
- green is about 500nm wavelength
- water refractive index for green at 0C is about 1.339
  http://optics.sgu.ru/_media/optics/staff/bashkatov/bashkatov_spie_03_5068_393.pdf
- perspex is about 1.48
"""
function r1_to_R(r1, exi::ExpImgs)
    p2m = exi.p2m
    r1_to_R(r1, p2m.h1, p2m.h2, p2m.h3, p2m.n1, p2m.n2, p2m.n3)
end
function r1_to_R(r1, h1, h2, h3, n1=1.0, n2=[1.339,1.48][1], n3=1.31)
    α = atan(r1/h1)
    β = asin(sin(α)*n1/n2)
    γ = asin(n2/n3*sin(β))
    Δr2 = tan(β)*h2
    r3 = (tan(γ)*h3 + r1 + Δr2) / (1+(tan(γ))^2)
    return r3/cos(γ)
end


text:
# This file was generated, do not modify it. # hide
ms = models() do m
    Vector{Continuous} <: m.target_scitype
end
foreach(m -> println(m.name), ms)


text:
module Models

include("Layers.jl")
include("Optimizer.jl")
include("DataSet.jl")
include("Func.jl")

export MLP,SimpleRNN,ManyToOneRNN,predict,learn,optimizer,shaping_rnn,save

using Random,JSON

#==============================
共通メンバ変数
・params : 各レイヤに重みをリスト化
・grads : 各レイヤに勾配をリスト化 → Layers.forward()でLayerのメンバに代入したら合わせて変更される
・layers : 各レイヤをリスト化
・loss_layer : 損失関数レイヤ
・optimizer : 最適化アルゴリズム

共通メソッド
・モデル名(;neuron,loss_layer) :コンストラクタ
    ・neuron : 入力～中間～出力　までのニューロン数を指定
    ・loss_layer : 損失関数（フルネーム指定）
        ・重み，バイアス，レイヤ，損失関数レイヤの初期化
        ・生成したレイヤからparams,gradsを取り出しメンバ変数へ保管

・predict(model_params;data) : 推論　→　スコア算出

・learn(model_params;batch_size,max_epoch,data,t_data=[])
    ※以下をmax_epoch回ループ,max_ite回ループ
    ・predict()実行　→　損失算出(loss_layerにt_dataが指定してあれば損失を返す)
    ・Model.layersを逆からbackward() → 勾配算出
    ・params更新

===============================#

#========
Model_params
========#
abstract type Model_params end
mutable struct MLP_params <:Model_params
    params
    grads
    layers::Array
    optimizer
    debug::Bool
end
mutable struct SimpleRNN_params <:Model_params
    params
    grads
    layers::Array
    rnn_layer
    optimizer
end
mutable struct ManyToOneRNN_params <: Model_params
    params
    grads
    layers
    optimizer
    padding
end

#==========
コンストラクタ
===========#
function MLP(;activation="Sigmoid",loss="Sigmoid_with_CrossEntropy",neuron,debug=false)
    #neuronには入力，中間1，中間2，...中間n，出力までの各ニューロン数を格納

    if debug
        Random.seed!(123456)
    end

    optimizer = []
    params = []
    grads = []

    #レイヤ初期化
    layers = []
    for i in 1:length(neuron)-1
        layers_tmp = []
        #パラメータ生成
        W = randn(neuron[i],neuron[i+1])
        b = randn(neuron[i+1],1)

        if i == length(neuron)-1
            #活性化関数&損失関数レイヤ
            activation_a = string("Layers.",loss,"()")
        else
            #活性化関数レイヤ
            activation_a = string("Layers.",activation,"()")
        end
        activation_layer = Meta.parse(activation_a)

        #レイヤ生成
        layers_tmp = [
            Layers.Affine(W,b),
            eval(activation_layer),
       ]

       #params,gradsをリストへ
       append!(params,layers_tmp[1].params) #Affineのparamsを取り出し
       append!(grads,layers_tmp[1].grads) #Affineのgradsを取り出し

       append!(layers,layers_tmp)
    end

    return MLP_params(params,grads,layers,optimizer,debug)
end
function SimpleRNN(;data_size,hidden_size,output_size)
    #hidden_size RNNレイヤの隠れ状態の次元数
    D, H, O = data_size, hidden_size, output_size

    #重みの初期化
    rnn_Wx = randn(D,H)/sqrt(D)
    rnn_Wh = randn(H,H)/sqrt(H)
    rnn_b = zeros(1,H)
    affine_W = randn(H,O)/sqrt(H)

    #その他パラメータ
    stateful = true
    padding = nothing

    affine_b = zeros(1,O)

    #勾配の初期化は各レイヤで行う

    #レイヤの初期化
    layers = [
        Layers.TimeRNN(rnn_Wx,rnn_Wh,rnn_b,stateful,padding),
        Layers.TimeAffine(affine_W,affine_b),
        Layers.TimeMse()
    ]

    #params,gradsをリストへ
    params = vcat(layers[1].params,layers[2].params)
    grads = vcat(layers[1].grads,layers[2].grads)

    optimizer = []

    return SimpleRNN_params(params,grads,layers,layers[1],[])
end
function ManyToOneRNN(input_size,hidden_size,output_size)
    #p : SimpleRNNとの差異

    #hidden_size RNNレイヤの隠れ状態の次元数
    D, H, O = input_size, hidden_size, output_size

    #その他パラメータ
    stateful = true
    padding = nothing

    #重みの初期化
    # rnn_Wx = randn(D,H)/sqrt(D)
    # rnn_Wh = randn(H,H)/sqrt(H)
    # rnn_b = zeros(1,H)
    rnn_Wx = randn(D,4*H)/sqrt(D)
    rnn_Wh = randn(H,4*H)/sqrt(H)
    rnn_b = zeros(1,4*H)

    affine_W = randn(H,O)/sqrt(H)
    affine_b = zeros(O,1)

    #レイヤの初期化
    layers = [
        Layers.TimeLSTM(rnn_Wx,rnn_Wh,rnn_b,true,padding),
        Layers.Affine(affine_W,affine_b), #p

        Layers.Mean_Squared_Error() #p
    ]

    #params,gradsをリストへ
    params = vcat(layers[1].params,layers[2].params)
    grads = vcat(layers[1].grads,layers[2].grads)

    optimizer = []
    padding = nothing

    return ManyToOneRNN_params(params,grads,layers,optimizer,padding)
end

function hoge()
    return 0
end

#================
推論
================#
function predict(params::MLP_params,data)
    out = data
    for i in 1:length(params.layers)
        out = Layers.forward(params.layers[i],out)
    end

    #推論と学習に応じて，スコアか損失か変わる(学習データの有無で判別)
    return out
end
function predict(params::SimpleRNN_params,data)
    hs = Layers.forward(params.layers[1],data)
    ys = Layers.forward(params.layers[2],hs)
    return ys
end
function predict(this::ManyToOneRNN_params,data)
    hs = Layers.forward(this.layers[1],data)
    #使うのはTの出力のみ
    h = hs[:,end,:]
    y = Layers.forward(this.layers[2],h)

    return y
end

#================
学習
================#
function learn(this::MLP_params;batch_size,max_epoch,data,t_data=[])
    ite = size(data)[1] ÷ batch_size
    epoch_loss = 0
    min_loss = 0
    min_loss_ep = 0
    loss_list = [] #エポック毎の損失の平均を格納
    #1エポック=データを1巡
    for i in 1:max_epoch
        ite_avg_loss_sum = 0 #1イテレーション毎の平均損失の総和
        #データのインデックスシャッフル
        idx = shuffle(1:size(data)[1])
        data = data[idx,:]
        t_data = t_data[idx,:]
        epoch_avg_loss = 0

        for k in 1:ite
            batch_data = data[1+(k-1)*batch_size:k*batch_size,1:end]
            this.layers[end].t = t_data[1+(k-1)*batch_size:k*batch_size,1:end]
            #順伝播
            loss = predict(this,batch_data)
            #逆伝播
            dout=1
            for j in length(this.layers):-1:1
                dout = Layers.backward(this.layers[j],dout)
            end
            ite_avg_loss_sum += sum(loss)/length(loss) #1イテレーションの平均損失を足していく

            #更新
            Optimizer.update(this.optimizer,this)
        end

        epoch_avg_loss = ite_avg_loss_sum/ite #1エポック毎の平均損失

        #10回に一回出力する
        if i%(max_epoch/10) == 0
            println("ep.$i : Loss :　",epoch_avg_loss)
            println("grads:",this.grads[1][1])
        end

        if i==1 || min_loss > epoch_avg_loss
            min_loss = epoch_avg_loss
            min_loss_ep = i
        end

        append!(loss_list,epoch_avg_loss)
    end

    println("min_loss : $min_loss")
    println("min_loss_ep : $min_loss_ep")

    #predict単体で実行するために初期化しておく
    this.layers[end].t = []

    return loss_list
end
function learn(this::SimpleRNN_params;batch_size,max_epoch,window_size,data,t_data=[])
    X = size(data,1) #総データ数
    D = size(data,2) #データ次元数
    T = window_size #RNNレイヤ数
    N = batch_size #バッチ数

    max_ite = (X/N)*(1/T) #イテレーション数
    loss_list = [] #avg_lossのリスト

    for epoch in 1:max_epoch
        total_loss = 0 #損失合計
        avg_loss = 0 #1エポックの平均損失
        for ite in 1:max_ite
            #ミニバッチ作成
            xs = zeros(Float64,(N,T,D))
            ts = zeros(Float64,(N,T,D))
            for n in 1:N
                st = 1+(T*(ite-1))+(X/N)*(n-1)
                ed = T+(T*(ite-1))+(X/N)*(n-1)
                xs[Int(n),:,:] = data[Int(st):Int(ed),:]
                ts[Int(n),:,:] = t_data[Int(st):Int(ed),:]
            end
            #順伝播
            ys = predict(this,xs)
            loss = Layers.forward(this.layers[end],ys,ts)
            total_loss += sum(loss)
            #逆伝播
            dl = ones(Float64,size(loss))
            dys = Layers.backward(this.layers[end],dl)
            dhs = Layers.backward(this.layers[2],dys)
            dxs = Layers.backward(this.layers[1],dhs)
            #更新
            Optimizer.update(this.optimizer,this)
        end

        avg_loss = total_loss/(T*D*N)
        append!(loss_list,avg_loss)

        #10回に一回出力する
        if epoch%(max_epoch/10) == 0
            println("ep.$epoch : Loss :　",avg_loss)
            println("grads:",this.grads[1][1])
        end
    end

    return loss_list
end
function learn(this::ManyToOneRNN_params;max_epoch,window_size,data,t_data)
    #先にshapingでデータ加工
    #Tはsize(data,2)を割れる値にする

    D = size(data,3) #データ次元数
    T = window_size #RNNレイヤ数
    N = size(data,1) #バッチ数

    max_ite = size(data,2)÷T #イテレーション数
    loss_list = [] #avg_lossのリスト

    grads_list = []

    for epoch in 1:max_epoch
        ite_total_loss = 0 #損失合計
        avg_loss = 0 #1エポックの平均損失
        st = 0 #data切り取り位置
        ed = 0
        for ite in 1:max_ite
            #ミニバッチ作成
            st = Int(1+(ite-1)*T)
            ed = Int(T*ite)
            xs = data[:,st:ed,:]
            t = t_data[:,ite,:]
            #順伝播
            y = predict(this,xs)
            this.layers[end].t = t #教師データ挿入
            loss = Layers.forward(this.layers[end],y)
            #ite毎の平均損失を加算
            ite_total_loss += sum(loss)/length(loss)
            #逆伝播
            dy = Layers.backward(this.layers[end],0)
            dh = Layers.backward(this.layers[2],dy)
            dxs = Layers.backward(this.layers[1],dh)
            #勾配クリッピング
            # append!(grads_list,this.grads[1][1])
            #更新
            Optimizer.update(this.optimizer,this)
        end

        avg_loss = ite_total_loss/max_ite
        append!(loss_list,avg_loss)

        if epoch == 1 || epoch%(max_epoch÷10) == 0
            println("ep.$epoch : Loss :　",avg_loss)
            # println("grads:",this.grads[1][1])
            # println("grads_l:",this.layers[1].grads[2][1])
        end
    end

    #学習の最後に隠れベクトルをリセット->predictする際，ミニバッチサイズが異なる為
    Layers.reset_state(this.layers[1])
    this.layers[end].t = [] #正解データリセット

    return loss_list, grads_list
end

#================
最適化 optimizer
=================#
function optimizer(this::Model_params;name="Adam",learning_rate=0.001,p1=0.95,p2=0.99,e=10^(-12))
    vs = []
    ss = []
    for i in 1:length(this.params)
        v = zeros(Float64,size(this.params[i]))
        s = zeros(Float64,size(this.params[i]))
        append!(vs,[v])
        append!(ss,[s])
    end
    hyper_params = [
        learning_rate,
        p1,
        p2,
        e,
        vs,
        ss
    ]
    #optimizer生成
    optimizer_a = string("Optimizer.$name($hyper_params)")
    this.optimizer = eval(Meta.parse(optimizer_a))
end

#================
データ整形
=================#
function shaping_rnn(model::ManyToOneRNN_params,data,N,padding)
    model.padding = padding
    #re_data (N, X/N, D)
    #ite毎に切り出して使用する
    #Tは2軸目のサイズに応じて後から決める

    X,D = size(data,1),size(data,2)
    re_data = []

    #データの補充
    if size(data,1)%N != 0
        #(n,D)分，paddingで埋める
        n = N - size(data,1)%N
        padding_arr = fill(padding,n,D)
        data = vcat(data,padding_arr) #補充
        X,D = size(data)
    end

    re_data = reshape(data,Int(X/N),N,D)
    re_data = permutedims(re_data,(2,1,3)) #軸の入れ替え

    return re_data
end
function get_rate_of_change(array)
    rate_arr = zeros(length(array)-1)
    for i in 1:length(array)-1
        rate = (array[i+1]-array[i])/array[i]
        rate_arr[i] = rate
    end
    #行列化
    rate_arr = reshape(rate_arr,(length(rate_arr),1))
    return rate_arr
end

#================
サンプルデータ
=================#
function dataset(type::String,n::Int64)
    type_name = string("DataSet.",type,"($n)")
    type_function = Meta.parse(type_name)
    return eval(type_function)
end

#=================
その他
==================#
function save(model::Model_params,file)
    JSON.json(model)
    f = open("$file","w") #ファイル生成
    JSON.print(f, model)
    close(f)
end

function load(model::Model_params,file)
    if !isfile(file)
        println("no file.")
        return nothing
    end

    f = open("$file","r")
    model_a = JSON.parse(f)

    model.params = model_a["params"]
    model.grads = model_a["grads"]
end

end


text:
type TwContributor
  id::Int64 #ID of user who contributed to this Tweet
  id_str::AbstractString #String representation of ID
  screen_name::AbstractString #Screen Name

  TwContributor(d::Dict{AbstractString, Any}) = new(
    get(d, "id", 0),
    get(d, "id_str", "0"),
    get(d, "screen_name", "")
  )
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
using BSON
using BSONMmap
using Test

dict = Dict(
    "str" => "abcdef",
    "nt" => (α = 1, β = 2, γ = "β"),
    "arr[uint8]" => zeros(UInt8, 7),
    "arr[str]" => ["α" "β"; "γ" "δ"],
    "arr[nt]" => repeat([(a = 1, b = "α", c = 1f0)], 2),
    "group" => Dict(
        "arr[float32]" => zeros(Float32, 11),
        "group" => Dict("arr[float64]" => ones(10)),
        )
    )
bssave("dict.bson", dict)
@test bsload("dict.bson", mmaparrays = false) == dict
@test bsload("dict.bson", mmaparrays = true) == dict

mutable struct Data{TX, TY}
    x::TX
    y::TY
    z::Dict{String, Int}
    w::String
end
data = Data(rand(Float32, 7, 3), rand(2, 2, 2), Dict("a" => 1), "abcd")
bssave("data.bson", data, force = true)
for mmap in (true, false)
    data′ = bsload("data.bson", Data, mmaparrays = mmap)
    for s in fieldnames(Data)
        @test getfield(data, s) == getfield(data′, s)
    end
end

if !Sys.iswindows()
    rm("data.bson", force = true)
    rm("dict.bson", force = true)
end


text:
include("conjugategradient.jl")
export HagerZhang


text:
using ArnoldiMethod
using ArnoldiMethodTransformations
using LinearAlgebra
using SparseArrays
using Statistics
using Test

N = 20
# construct fixed eval matrix in random basis
D = diagm(0=>0:(N-1))

@testset "ordinary eigenvalue problem: " begin
    for i ∈ 1:500
        S = I + .1*randn(N,N) # random basis

        A = S\D*S
        # find eigenpairs closest to 5.001 (cannot be 5 as algorithm is unstable if σ is exactly an eval)
        decomp, hist = partialschur(A,5.001)
        # get evecs
        λ, v = partialeigen(decomp,5.001)
        @test norm(A*v-v*diagm(0=>λ))/prod(size(v)) ≤ sqrt(eps(float(N)))N

        As = sparse(A)
        # find eigenpairs closest to 5.001 (cannot be 5 as algorithm is unstable if σ is exactly an eval)
        decomp, hist = partialschur(As,5.001)
        # get evecs
        λ, v = partialeigen(decomp,5.001)
        @test norm(A*v-v*diagm(0=>λ))/prod(size(v)) ≤ sqrt(eps(float(N)))N
    end
end


text:
using GeometryTypes, GLVisualize, GLAbstraction, Reactive, GLWindow, GLFW, Colors

const description = """
You can drag the edges of this graph,
or remove them with a right click.
Camera is fixed in this example.
"""

if !isdefined(:runtests)
    window = glscreen()
end
window.color = RGBA(0.04f0, 0.02f0, 0.05f0, 1f0)

const record_interactive = true

n = 200
n_connections = 200
w, h = widths(window)
a = map(x->x.*Point2f0(w, h), rand(Point2f0, n))
indices = Signal(rand(1:n, n_connections))
# for the points, we only need unique indices!
unique_indices = map(unique, indices)
scales = map(Vec2f0, rand(5f0:15f0, n))
colors = zeros(RGBA{Float32}, length(a))

# when scale and position is given, its enough to pass the Circle type instead of
# a concrete instance
points = visualize(
    (Circle, a),
    indices=unique_indices,
    scale=scales, offset=map(x-> -x./2f0, scales),
    color=RGBA(1f0,1f0,1f0,1f0),
    stroke_width=3f0, stroke_color=RGBA(1f0,1f0,1f0,0.4f0),
    glow_width=4f0, glow_color=colors
)

const point_robj = points.children[] # temporary way of getting the render object. Shouldn't stay like this
 # best way to get the gpu object. One could also start by creating a gpu array oneself.
 # this is a bit tricky, since not there are three different types.
 # for points and lines you need a GLBuffer. e.g gpu_position = GLBuffer(rand(Point2f0, 50)*1000f0)
const gpu_position = point_robj[:position]
const gpu_scales = point_robj[:scale]
const gpu_colors = point_robj[:glow_color]
const gpu_offset = point_robj[:offset]

# GPUVector adds functionality like push! and splice!
# now the lines and points share the same gpu object
# for linesegments, you can pass indices, which needs to be of some 32bit int type
cmap = map(RGBA{Float32}, colormap("Blues", 50))
lines = visualize(
    gpu_position, :linesegment,
    indices=indices, color=cmap[rand(1:50, n)], # sample random colors from gradient
    thickness=0.5f0
)

# current tuple of renderobject id and index into the gpu array
const m2id = GLWindow.mouse2id(window)
isoverpoint = const_lift(is_same_id, m2id, point_robj)
v0 = (0, scales, colors)
time = bounce(linspace(0f0, 1f0, 20)) # bounces between 0 - 1, for the animation
preserve(foldp(v0, isoverpoint, time) do v0, overpoint, t
    index0, scale0, color0 = v0
    if overpoint
        id, index = value(m2id)
        if id==point_robj.id && length(a) >= index
            new_colors = zeros(RGBA{Float32}, length(a))
            gpu_colors[index] = RGBA{Float32}(0.95, 0.56, t, 0.7)
            newscale = scale0[index] + t*7f0
            gpu_scales[index] = newscale
            gpu_offset[index] = -newscale./2f0
            return index, scale0, color0
        end
    else
        if index0 > 0
            gpu_colors[index0] = color0[index0]
            gpu_scales[index0] = scale0[index0]
            gpu_offset[index0] = -scale0[index0]./2f0
        end
    end
    return v0
end)

# inputs are a dict, materialize gets the keys out of it (equivalent to mouseposition = window.inputs[:mouseposition])
@materialize mouse_buttons_pressed, mouseposition = window.inputs

# single left mousekey pressed (while no other mouse key is pressed)
key_pressed = const_lift(GLAbstraction.singlepressed, mouse_buttons_pressed, GLFW.MOUSE_BUTTON_LEFT)
# righ
right_pressed = const_lift(GLAbstraction.singlepressed, mouse_buttons_pressed, GLFW.MOUSE_BUTTON_RIGHT)
# dragg while key_pressed. Drag only starts if isoverpoint is true
mousedragg = GLAbstraction.dragged(mouseposition, key_pressed, isoverpoint)

# use mousedrag and mouseid + index to actually change the gpu array with the positions
preserve(foldp((value(m2id)..., Point2f0(0)), mousedragg) do v0, dragg
    if dragg == Vec2f0(0) # if drag just started. Not the best way, maybe dragged should return a tuple of (draggvalue, started)
        id, index = value(m2id)
        if id==point_robj.id && length(gpu_position) >= index
            p0 = gpu_position[index]
        else
            p0 = v0[3]
        end
    else
        id, index, p0 = v0
        if id==point_robj.id && length(gpu_position) >= index
            gpu_position[index] = Point2f0(p0) + Point2f0(dragg)
        end
    end
    return id, index, p0
end)
# On right click remove nodes!
s = map(right_pressed) do rp
    id, index = value(m2id)
    if rp && id==point_robj.id && length(gpu_position) >= index
        new_indices = Int[]
        indice_val = value(indices)
        for i=1:2:length(indice_val) #filter out indices
            a, b = indice_val[i], indice_val[i+1]
            if a != index && b != index
                push!(new_indices, a, b)
            end
        end
        push!(indices, new_indices) # update indices!
    end
end
preserve(s)
# _view it!
_view(lines, window, camera=:fixed_pixel)
_view(points, window, camera=:fixed_pixel)

if !isdefined(:runtests)
    renderloop(window)
end


text:
# Julia wrapper for header: rocrand_version.h
# Automatically generated using Clang.jl

# Julia wrapper for header: rocrand_discrete_types.h
# Automatically generated using Clang.jl

# Julia wrapper for header: rocrand.h
# Automatically generated using Clang.jl

function rocrand_create_generator(generator, rng_type)
    @check ccall((:rocrand_create_generator, librocrand), rocrand_status,
                 (Ptr{rocrand_generator}, rocrand_rng_type),
                 generator, rng_type)
end

function rocrand_destroy_generator(generator)
    @check ccall((:rocrand_destroy_generator, librocrand), rocrand_status,
                 (rocrand_generator,),
                 generator)
end

function rocrand_generate(generator, output_data, n)
    @check ccall((:rocrand_generate, librocrand), rocrand_status,
                 (rocrand_generator, Ptr{UInt32}, Cint),
                 generator, output_data, n)
end

function rocrand_generate_char(generator, output_data, n)
    @check ccall((:rocrand_generate_char, librocrand), rocrand_status,
                 (rocrand_generator, Ptr{Cuchar}, Cint),
                 generator, output_data, n)
end

function rocrand_generate_short(generator, output_data, n)
    @check ccall((:rocrand_generate_short, librocrand), rocrand_status,
                 (rocrand_generator, Ptr{UInt16}, Cint),
                 generator, output_data, n)
end

function rocrand_generate_uniform(generator, output_data, n)
    @check ccall((:rocrand_generate_uniform, librocrand), rocrand_status,
                 (rocrand_generator, Ptr{Cfloat}, Cint),
                 generator, output_data, n)
end

function rocrand_generate_uniform_double(generator, output_data, n)
    @check ccall((:rocrand_generate_uniform_double, librocrand), rocrand_status,
                 (rocrand_generator, Ptr{Cdouble}, Cint),
                 generator, output_data, n)
end

function rocrand_generate_uniform_half(generator, output_data, n)
    @check ccall((:rocrand_generate_uniform_half, librocrand), rocrand_status,
                 (rocrand_generator, Ptr{half}, Cint),
                 generator, output_data, n)
end

function rocrand_generate_normal(generator, output_data, n, mean, stddev)
    @check ccall((:rocrand_generate_normal, librocrand), rocrand_status,
                 (rocrand_generator, Ptr{Cfloat}, Cint, Cfloat, Cfloat),
                 generator, output_data, n, mean, stddev)
end

function rocrand_generate_normal_double(generator, output_data, n, mean, stddev)
    @check ccall((:rocrand_generate_normal_double, librocrand), rocrand_status,
                 (rocrand_generator, Ptr{Cdouble}, Cint, Cdouble, Cdouble),
                 generator, output_data, n, mean, stddev)
end

function rocrand_generate_normal_half(generator, output_data, n, mean, stddev)
    @check ccall((:rocrand_generate_normal_half, librocrand), rocrand_status,
                 (rocrand_generator, Ptr{half}, Cint, half, half),
                 generator, output_data, n, mean, stddev)
end

function rocrand_generate_log_normal(generator, output_data, n, mean, stddev)
    @check ccall((:rocrand_generate_log_normal, librocrand), rocrand_status,
                 (rocrand_generator, Ptr{Cfloat}, Cint, Cfloat, Cfloat),
                 generator, output_data, n, mean, stddev)
end

function rocrand_generate_log_normal_double(generator, output_data, n, mean, stddev)
    @check ccall((:rocrand_generate_log_normal_double, librocrand), rocrand_status,
                 (rocrand_generator, Ptr{Cdouble}, Cint, Cdouble, Cdouble),
                 generator, output_data, n, mean, stddev)
end

function rocrand_generate_log_normal_half(generator, output_data, n, mean, stddev)
    @check ccall((:rocrand_generate_log_normal_half, librocrand), rocrand_status,
                 (rocrand_generator, Ptr{half}, Cint, half, half),
                 generator, output_data, n, mean, stddev)
end

function rocrand_generate_poisson(generator, output_data, n, lambda)
    @check ccall((:rocrand_generate_poisson, librocrand), rocrand_status,
                 (rocrand_generator, Ptr{UInt32}, Cint, Cdouble),
                 generator, output_data, n, lambda)
end

function rocrand_initialize_generator(generator)
    @check ccall((:rocrand_initialize_generator, librocrand), rocrand_status,
                 (rocrand_generator,),
                 generator)
end

function rocrand_set_stream(generator, stream)
    @check ccall((:rocrand_set_stream, librocrand), rocrand_status,
                 (rocrand_generator, Cint),
                 generator, stream)
end

function rocrand_set_seed(generator, seed)
    @check ccall((:rocrand_set_seed, librocrand), rocrand_status,
                 (rocrand_generator, Culonglong),
                 generator, seed)
end

function rocrand_set_offset(generator, offset)
    @check ccall((:rocrand_set_offset, librocrand), rocrand_status,
                 (rocrand_generator, Culonglong),
                 generator, offset)
end

function rocrand_set_quasi_random_generator_dimensions(generator, dimensions)
    @check ccall((:rocrand_set_quasi_random_generator_dimensions, librocrand), rocrand_status,
                 (rocrand_generator, UInt32),
                 generator, dimensions)
end

function rocrand_get_version(version)
    @check ccall((:rocrand_get_version, librocrand), rocrand_status,
                 (Ptr{Cint},),
                 version)
end

function rocrand_create_poisson_distribution(lambda, discrete_distribution)
    @check ccall((:rocrand_create_poisson_distribution, librocrand), rocrand_status,
                 (Cdouble, Ptr{rocrand_discrete_distribution}),
                 lambda, discrete_distribution)
end

function rocrand_create_discrete_distribution(probabilities, size, offset,
                                              discrete_distribution)
    @check ccall((:rocrand_create_discrete_distribution, librocrand), rocrand_status,
                 (Ptr{Cdouble}, UInt32, UInt32, Ptr{rocrand_discrete_distribution}),
                 probabilities, size, offset, discrete_distribution)
end

function rocrand_destroy_discrete_distribution(discrete_distribution)
    @check ccall((:rocrand_destroy_discrete_distribution, librocrand), rocrand_status,
                 (rocrand_discrete_distribution,),
                 discrete_distribution)
end


text:
module MessyOutput
using OnlineStats, Base.Test, Distributions

x = randn(500)
x1 = randn(500)
x2 = randn(501)
xs = hcat(x1, x)

@testset "show methods" begin
    display(Mean(x))
    display(Means(xs))
    display(Variance(x))
    display(Variances(xs))
    display(CovMatrix(xs))
    display(Extrema(x))
    display(QuantileSGD(x))
    display(QuantileMM(x))
    display(Moments(x))
    display(QuantRegMM(5))
    display(KMeans(5, 4))
    display(NormalMix(4))
    display(FitCategorical())
    display(FitBeta())
    display(BernoulliBootstrap(Mean(), mean, 1000))
    display(StatLearn(4))

    @testset "Full Factorial of Combinations" begin
        n, p = 500, 5
        x = randn(n, p)
        β = collect(linspace(-1, 1, p))
        β_with_intercept = vcat(0.0, β)
        xβ = x*β
        alg = [SGD(), AdaGrad(), AdaGrad2(), AdaDelta(), ADAM()] #, RDA(), MMGrad()]
        pen = [NoPenalty(), RidgePenalty(.1), LassoPenalty(.1), ElasticNetPenalty(.1, .5)]
        mod = [
            LinearRegression(), L1Regression(), LogisticRegression(),
            PoissonRegression(), QuantileRegression(.5), SVMLike(), HuberRegression(2.)
        ]

        generate(::LinearRegression, xβ) = xβ + randn(size(xβ, 1))
        generate(::L1Regression, xβ) = xβ + randn(size(xβ, 1))
        generate(::LogisticRegression, xβ) = [rand(Bernoulli(1 / (1 + exp(-η)))) for η in xβ]
        generate(::PoissonRegression, xβ) = [rand(Poisson(exp(η))) for η in xβ]
        generate(::QuantileRegression, xβ) = xβ + randn(size(xβ, 1))
        generate(::SVMLike, xβ) = [rand(Bernoulli(1 / (1 + exp(-η)))) for η in xβ]
        generate(::HuberRegression, xβ) = xβ + randn(size(xβ, 1))

        for a in alg, p in pen, m in mod
            y = generate(m, xβ)
            println(replace("    > $a, $p, $m", "OnlineStats.", ""))
            StatLearn(x, y, m, a, p)
            StatLearn(x, y, 10, m, a, p)
        end
    end
end

end#module


text:
@inline function biasedEstimate(model::Model, rpio::RPIO{Particle}) where Particle
  potentials::Vector{Float64} = rpio.potentials
  v::Float64 = 0.0
  _sampleParticles(model, rpio)
  for p in 1:rpio.n
    v += log(_evaluatePotentials(model, rpio, p) / rpio.N)
  end
  return v
end

@inline function biasedEstimateDiscrete(model::Model, rpio::RPIO{Int64Particle})
  _sampleParticles(model, rpio)
  ζcounts::Vector{Int64} = _getCounts(rpio)
  potentials::Vector{Float64} = Vector{Float64}(undef, length(ζcounts))
  v::Float64 = 0.0
  for p in 1:rpio.n
    s::Float64 = _countPotentials(potentials, ζcounts, model, rpio, p)
    v += log(s / rpio.N)
  end
  return v
end


text:
export PoolDims

"""
    PoolDims(x_size::NTuple{M}, k::Union{NTuple{L, Int}, Int};
            stride=k, padding=0, dilation=1)  where {M, L}

Dimensions for a "pooling" operation that can have an arbitrary input size, kernel size,
stride, dilation, and channel count.  Used to dispatch onto efficient implementations at
compile-time.
"""
struct PoolDims{N, K, S, P, D} <: ConvDims{N}
    input_size::NTuple{N, Int}

    kernel_size::NTuple{K, Int}
    channels_in::Int

    stride::NTuple{S, Int}
    padding::NTuple{P, Int}
    dilation::NTuple{D, Int}
end

function PoolDims(
    x_size::NTuple{M}, k::Union{NTuple{L, Int}, Int};
    stride = k, padding = 0, dilation = 1,
) where {M, L}
    _check_kernel(k::Number, N::Int) = ntuple(_ -> Int(k), N)
    _check_kernel(k::NTuple, ::Int) = k

    kernel = _check_kernel(k, M - 2)
    spdf_kernel = NTuple{M, Int}([kernel..., 1, 1])

    sstride, ppadding, ddilation = check_spdf(
        x_size, spdf_kernel, stride, padding, dilation)
    PoolDims(
        x_size[1:(end - 2)], kernel, x_size[end - 1],
        sstride, ppadding, ddilation)
end

PoolDims(x::AbstractArray, k; kwargs...) = PoolDims(size(x), k; kwargs...)

# Useful for constructing a new PoolDims that has only a few elements different
# from the original progenitor object that it inherits shapes from.
PoolDims(
    c::C; I=input_size(c), K=kernel_size(c),
    C_in=channels_in(c), S=stride(c), P=padding(c), D=dilation(c),
) where C <: ConvDims = PoolDims(I, K, C_in, S, P, D)

@inline channels_in(c::PoolDims) = c.channels_in
@inline channels_out(c::PoolDims) = c.channels_in

@inline input_size(c::PoolDims) = c.input_size
@inline kernel_size(c::PoolDims) = c.kernel_size

@inline stride(c::PoolDims) = c.stride
@inline padding(c::PoolDims) = c.padding
@inline dilation(c::PoolDims) = c.dilation
@inline flipkernel(c::PoolDims) = false

function check_dims(x::NTuple{M}, y::NTuple{M}, pdims::PoolDims) where {M}
    # First, check that channel counts are all correct:
    @assert x[end-1] == channels_in(pdims) DimensionMismatch("Data input channel count ($(x[end-1]) vs. $(channels_in(pdims)))")
    @assert y[end-1] == channels_out(pdims) DimensionMismatch("Data output channel count ($(y[end-1]) vs. $(channels_out(pdims)))")

    # Next, check that the spatial dimensions match up
    @assert x[1:end-2] == input_size(pdims) DimensionMismatch("Data input spatial size ($(x[1:end-2]) vs. $(input_size(pdims)))")
    @assert y[1:end-2] == output_size(pdims) DimensionMismatch("Data output spatial size ($(y[1:end-2]) vs. $(output_size(pdims)))")

    # Finally, check that the batch size matches
    @assert x[end] == y[end] DimensionMismatch("Batch size ($(x[end]) vs. $(y[end]))")
end


text:
using DrWatson
@quickactivate
using GenerativeAD
using BSON, FileIO
using ImageTransformations

outpath = datadir("mvtec_ad/downscaled_data")
mkpath(outpath)

# start with wood
function downscale_data(category, size)
	normal, anomalous = GenerativeAD.Datasets.load_mvtec_ad_data(category=category)
	dnormal, danomalous = map(x->imresize(x,(size,size)), (normal, anomalous))
	fname = joinpath(outpath, "$(category)_$(size).bson")
	save(fname, Dict(:normal => dnormal, :anomalous => danomalous))
	@info "saved to $fname"
end

map(f->map(c->downscale_data(c,f), ["wood", "transistor", "grid"]), [256, 128, 64, 32])


text:
# Single RefTx entry
struct RefTx
   info::TxInfo
   don::CoordTuple
   acc::CoordTuple
   length::Float64
end

# Single RefGene entry
struct RefGene
   info::GeneInfo
   don::CoordTuple
   acc::CoordTuple
   txst::CoordTuple
   txen::CoordTuple
   exons::CoordTree
   noveldon::CoordTuple
   novelacc::CoordTuple
   exonexpr::Float64
   length::Float64
   reftx::Vector{RefTx}
end

const RefSet = Dict{GeneName,RefGene}

# this guy maps an array of strings or substrings to Int parsing and then into a tuple
# which is then sliced based on the l (left) and r (right) adjusters, c is the mathmatical adjuster to the data
parse_splice( subarray; l=0, r=0, c=0 ) = tuple(map( x->convert(CoordInt,parse(Int,x)+c), subarray )...)[(1+l):(end-r)]
parse_coordint( i ) = convert(CoordInt, parse(Int, i))
parse_coordtup( i; c=0 ) = tuple(parse_coordint(i)+CoordInt(c))

unique_tuple( tup1::Tuple{}, tup2::Tuple{} ) = ()
unique_tuple( tup1::Tuple{Vararg{T}}, tup2::Tuple{} ) where T = unique_tuple(tup1, tup1)
unique_tuple( tup1::Tuple{}, tup2::Tuple{Vararg{T}} ) where T = unique_tuple(tup2, tup2)

function unique_tuple( tup1::Tuple{Vararg{T}}, tup2::Tuple{Vararg{T}}) where T
   uniq = SortedSet{T,Base.Order.ForwardOrdering}()
   for i in tup1, j in tup2
     push!(uniq, i)
     push!(uniq, j)
   end
   tuple(uniq...)
end

function minimum_threshold!( dict::Dict{K,V}, limit::V, range::UnitRange ) where {K, V <: Number}
   for k in collect(keys(dict))
      if dict[k] < limit || !(k in range)
         delete!( dict, k )
      end
   end
end

function load_refflat( fh; txbool=true )
   error("ERROR: Annotation files in refflat format have been deprecated as of Whippet v0.11, please use --gtf!")
end

function fetch_meta( var::String, meta; off=1 )
   val = "_" # empty
   for i in off:2:length(meta)
      if meta[i] == var && i <= length(meta)
         return string(meta[i+1]),i+1
      end
   end
   val,0
end

function load_gtf( fh; txbool=true, suppress=false, usebam=false, bamreader=Nullable{BAM.Reader}(), bamreads=2, bamoneknown=false )

   # Temporary variables
   gninfo   = Dict{GeneName,GeneInfo}()
   gndon    = Dict{GeneName,CoordTuple}()
   gnacc    = Dict{GeneName,CoordTuple}()
   gntxst   = Dict{GeneName,CoordTuple}()
   gntxen   = Dict{GeneName,CoordTuple}()
   gnlens   = Dict{GeneName,CoordTuple}()
   gnexons  = Dict{GeneName,CoordTree}()
   gnreftx  = Dict{GeneName,Vector{RefTx}}()
   gnlen    = Dict{GeneName,Float64}()
   gncnt    = Dict{GeneName,Int}()

   # RefSet variables
   geneset  = Dict{GeneName,RefGene}()

   used_txn = Dict{GeneName,Bool}()
   warning  = false
   warning_num = 0

   txnum = 75000
   sizehint!(geneset, txnum >> 1)

   curtran = ""
   curgene = ""
   curchrom = ""
   curstran = '+'

   trandon = Vector{CoordInt}()
   tranacc = Vector{CoordInt}()
   txlen   = 0

   function private_add_transcript!( curtran, curgene, curchrom, curstran, trandon, tranacc, txlen )
      sort!(trandon)
      sort!(tranacc)

      txS = minimum(tranacc)
      txE = maximum(trandon)

      txinfo = TxInfo(curtran, txS, txE, UInt16(length(tranacc)))

      don = tuple(trandon[1:(end-1)]...)
      acc = tuple(tranacc[2:end]...)

      if haskey(gninfo, curgene)
         curchrom == gninfo[curgene].name || return # can have only one chrom
         gndon[curgene]  = unique_tuple(gndon[curgene], don)
         gnacc[curgene]  = unique_tuple(gnacc[curgene], acc)
         gntxst[curgene] = unique_tuple(gntxst[curgene], tuple(txS))
         gntxen[curgene] = unique_tuple(gntxen[curgene], tuple(txE))
         gnlen[curgene] += txlen
         gncnt[curgene] += 1
         push!( gnreftx[curgene], RefTx( txinfo, tuple(trandon...), tuple(tranacc...), txlen ) )
      else
         gndon[curgene]  = don
         gnacc[curgene]  = acc
         gninfo[curgene] = GeneInfo(curgene, string(curchrom), curstran)
         gntxst[curgene] = tuple(txS)
         gntxen[curgene] = tuple(txE)
         gnlen[curgene]  = txlen
         gncnt[curgene]  = 1
         gnreftx[curgene]  = RefTx[ RefTx( txinfo, tuple(trandon...), tuple(tranacc...), txlen ) ]
      end
   end

   for l in eachline(fh)

      (l[1] == '#') && continue # ignore comment lines

      (chrom, src, entrytype,
       st, en, _,
       strand, _,
       meta) = split(chomp(l), '\t')

      (entrytype == "exon") || continue

      metaspl = split(meta, [' ',';','"'], keep=false)

      geneid,pos  = fetch_meta( "gene_id", metaspl )
      tranid,pos  = fetch_meta( "transcript_id", metaspl )
      genesym     = fetch_meta( "gene_name", metaspl )
      support,val = fetch_meta( "transcript_support_level", metaspl )

      # if we observe low transcript support levels, we will warn the user
      if (support != "_" && support != "1" && support != "2" &&  support != "NA")

         if suppress
            continue
         elseif !warning
            println(STDERR, "")
            warn("Using low quality Transcript Support Levels (TSL 3+) in your GTF file is not recommended!\nFor more information on TSL, see: http://www.ensembl.org/Help/Glossary?id=492\n\nIf you would like Whippet to ignore these when building its index, use `--suppress-low-tsl` option!\n\n")

            warning=true
         end

      elseif haskey(used_txn, tranid)
         error("ERROR: GTF file is not in valid GTF2.2 format!\n\nERROR: Annotation entries for 'transcript_id' $tranid has already been fully processed and closed.\nHINT: All GTF lines with the same 'transcript_id' must be adjacent in the GTF file and referring to the same transcript and gene!")
      elseif tranid == geneid && tranid != curtran
         if warning_num < 25
            warn("Generally 'transcript_id' should not equal 'gene_id' but does at $tranid == $geneid;")
         elseif warning_num == 25
            warn("... NOTE: 'transcript_id' == 'gene_id' will work OK for single isoform genes, but will not produce expected behavior for multi-isoform genes!\n... similar warnings will be suppressed; disregard if 'transcript_id' == 'gene_id' is intentional\n")
         end
         warning_num += 1
      end

      if tranid != curtran
         # add transcript and clean up
         curtran != "" && private_add_transcript!( curtran, curgene, curchrom, curstran, trandon, tranacc, txlen )

         used_txn[curtran] = true

         curtran = tranid
         curgene = geneid
         curchrom = chrom
         curstran = strand[1]

         empty!(trandon)
         empty!(tranacc)
         txlen = 0

      end

      sval = parse_coordint(st)
      eval = parse_coordint(en)

      push!(tranacc, sval)
      push!(trandon, eval)
      txlen += eval - sval + 1

      insval = Interval{CoordInt}(sval, eval)
      if haskey(gnexons, curgene)
         # make sure we are adding a unique value
         if !haskey(gnexons[curgene], (sval, eval))
            push!(gnexons[curgene], insval)
         end
      else
         gnexons[curgene] = CoordTree()
         push!(gnexons[curgene], insval)
      end
   end

   private_add_transcript!( curtran, curgene, curchrom, curstran, trandon, tranacc, txlen )

   if usebam
      bamnames = map(x->values(x)[1], find(BAM.header(get(bamreader)), "SQ"))
   end

   novel_ss = 0
   annot_ss = 0
   # iterate through genes, if usebam, add novel splice sites and calculate relative exonexpr
   # set RefGene at the end of each iteration
   for gene in keys(gninfo)

      novelacc = Dict{CoordInt,Int}()
      noveldon = Dict{CoordInt,Int}()
      exonexpr = 0.0
      meanleng = gnlen[gene] / gncnt[gene]
      seqname  = gninfo[gene].name
      strand   = gninfo[gene].strand
      bamwasused = false

      annot_ss += length(gndon[gene]) + length(gnacc[gene])

      if usebam && seqname in bamnames
         bamwasused = true
         range   = Int64(minimum(gntxst[gene])):Int64(maximum(gntxen[gene]))
         exoncount = process_records!( get(bamreader), seqname, range, strand,
                                       gnexons[gene], union(gnacc[gene], gndon[gene]), bamoneknown,
                                       novelacc, noveldon )
         exonexpr  = exoncount / (meanleng - 100)

         # clean up cryptic splice-sites
         minimum_threshold!( novelacc, bamreads, range )
         minimum_threshold!( noveldon, bamreads, range )
         novelacctup = map(CoordInt, Tuple(collect(keys(novelacc))))
         noveldontup = map(CoordInt, Tuple(collect(keys(noveldon))))

         # add splice-sites from bam
         annoacc     = gnacc[gene]
         annodon     = gndon[gene]
         gnacc[gene] = unique_tuple( gnacc[gene], novelacctup )
         gndon[gene] = unique_tuple( gndon[gene], noveldontup )
         novel_ss += length(gnacc[gene]) + length(gndon[gene]) - (length(annoacc) + length(annodon))
      end

      # clean up any txst or txen that match a splice site
      gntxst[gene] = Tuple( setdiff( gntxst[gene], gnacc[gene] ) )
      gntxen[gene] = Tuple( setdiff( gntxen[gene], gndon[gene] ) )

      geneset[gene] = RefGene( gninfo[gene],
                               gndon[gene],
                               gnacc[gene],
                               gntxst[gene],
                               gntxen[gene],
                               gnexons[gene],
                               bamwasused ? Tuple( setdiff(noveldontup, annodon) ) : Tuple{}(),
                               bamwasused ? Tuple( setdiff(novelacctup, annoacc) ) : Tuple{}(),
                               exonexpr,
                               meanleng,
                               gnreftx[gene] )
   end

   if usebam
      println(STDERR, "Found $novel_ss new splice-sites from BAM file..")
   end
   println(STDERR, "Loaded $annot_ss annotated splice-sites from GTF file..")

   return geneset
end


text:
using IntervalMatrices: correction_hull

# ==================================
# Abstract interface
# ==================================

abstract type AbstractApproximationModel end

function _default_approximation_model(ivp::IVP{<:AbstractContinuousSystem})
    return Forward()
end

# some algorithms require a polyhedral computations backend
hasbackend(alg::AbstractApproximationModel) = false

# symmetric inteval hull options
sih(X, ::Val{:lazy}) = SymmetricIntervalHull(X)
sih(X, ::Val{:concrete}) = _symmetric_interval_hull(X)

# interval matrix functions
isinterval(A::AbstractMatrix{N}) where {N<:Number} = false
isinterval(A::IntervalMatrix{N, IT}) where {N, IT<:IA.Interval{N}} = true
isinterval(A::AbstractMatrix{IT}) where {IT<:IA.Interval} = true

# options for a-posteriori transformation of a discretized set
# valid options are:
# AbstractDirections, Val{:lazy}, Val{:concrete}, Val{:vrep}, Val{:zono}, Val{:zonotope}
# _alias(setops) = setops # no-op

_alias(setops::AbstractDirections) = setops
_alias(setops::Val{:lazy}) = setops
_alias(setops::Val{:concrete}) = setops
_alias(setops::Val{:vrep}) = setops
_alias(setops::Val{:box}) = setops
_alias(setops::Val{:zono}) = setops
_alias(setops::Val{:zonotope}) = Val(:zono)

"""
    discretize(ivp::IVP, δ, alg::AbstractApproximationModel)

Set-based conservative discretization of a continuous-time initial value problem
into a discrete-time problem.

### Input

- `ivp`   -- initial value problem for a linear ODE in canonical form (see `Notes` below)
- `δ`     -- step size
- `alg`   -- algorithm used to compute the approximation model

### Output

The initial value problem of a discrete system.

### Notes

Different approximation algorithms and their respective options are described
in the docstring of each method, e.g. [`Forward`](@ref). Here is a list of all
the available approximation models:

```jldoctest
julia> subtypes(ReachabilityAnalysis.AbstractApproximationModel)
5-element Vector{Any}:
 Backward
 CorrectionHull
 Forward
 NoBloating
 StepIntersect
```

Initial-value problems considered in this function are of the form

```math
x' = Ax(t) + u(t),\\qquad x(0) ∈ \\mathcal{X}_0,\\qquad (1)
```
and where ``u(t) ∈ U(k)`` add where ``\\{U(k)\\}_k`` is a sequence of sets of
non-deterministic inputs and ``\\mathcal{X}_0`` is the set of initial
states. Other problems, e.g. ``x' = Ax(t) + Bu(t)`` can be brought
to the canonical form with the function [`normalize`](@ref).

For references to the original papers introducing each algorithm, see the docstrings,
e.g. `?Forward`.
"""
function discretize(ivp::IVP, δ, alg::AbstractApproximationModel) end

# =========================================
# Conservative time discretization methods
# =========================================

# Approximation model in discrete time, i.e. without bloating
include("NoBloating.jl")

# Forward approximation
include("Forward.jl")

# Backward approximation
include("Backward.jl")

# Discretize using the correction hull of the matrix exponential
include("CorrectionHull.jl")

# Intersect one step forward in time with one step backward
include("StepIntersect.jl")

# =========================================================================
# Alternatives to apply the set operation depending on the desired output
# =========================================================================

function _apply_setops(X, alg::Forward)
    if hasbackend(alg)
        _apply_setops(X, alg.setops, alg.backend)
    else
        _apply_setops(X, alg.setops)
    end
end

_apply_setops(X::LazySet, ::Val{:lazy}) = X  # no-op
_apply_setops(X::LazySet, ::Val{:concrete}) = concretize(X) # concrete set
_apply_setops(X, template::AbstractDirections) = overapproximate(X, template) # template oa
_apply_setops(X::LazySet, ::Val{:box}) = box_approximation(X) # box oa

_apply_setops(M::AbstractMatrix, X::LazySet, ::Val{:lazy}) = M * X
_apply_setops(M::AbstractMatrix, X::LazySet, ::Val{:concrete}) = linear_map(M, X)

# evantually we should use concretize, but requires fast fallback operations in 2D
# such as Minkowski sum not yet available
function _apply_setops(X::ConvexHull{N, AT, MS}, ::Val{:vrep}, backend=nothing) where {N,
                            AT<:AbstractPolytope{N}, MT,
                            LM<:LinearMap{N, AT, N, MT},
                            BT<:AbstractPolytope, MS<:MinkowskiSum{N, LM}}
    n = dim(X)
    VT = n == 2 ? VPolygon : VPolytope

    # CH(A, B) := CH(X₀, ΦX₀ ⊕ E₊)
    A = X.X
    B = X.Y
    X₀ = convert(VT, A)

    if n == 2
        ΦX₀ = convert(VT, B.X)
        E₊ = convert(VT, B.Y)
        out = convex_hull(X₀, minkowski_sum(ΦX₀, E₊))
    else
        # generic conversion to VPolytope is missing, see LazySets#2467
        ΦX₀ = VPolytope(vertices_list(B.X, prune=false))
        E₊ = convert(VT, B.Y)
        aux = minkowski_sum(ΦX₀, E₊, apply_convex_hull=false)
        out = convex_hull(X₀, aux, backend=backend)
    end

    return out
end

# give X = CH(X₀, ΦX₀ ⊕ E₊), return a zonotope overapproximation
function _apply_setops(X::ConvexHull{N, AT, MS}, ::Val{:zono}, backend=nothing) where {N,
                            AT<:AbstractZonotope{N}, MT,
                            LM<:LinearMap{N, AT, N, MT},
                            BT<:AbstractZonotope, MS<:MinkowskiSum{N, LM}}
    # CH(A, B) := CH(X₀, ΦX₀ ⊕ E₊)
    A = X.X
    B = X.Y
    return overapproximate(CH(A, concretize(B)), Zonotope)
end


text:
const BayesianLDA_ = skda(:LinearDiscriminantAnalysis)
@sk_clf mutable struct BayesianLDA <: MMI.Probabilistic
    solver::String                   = "svd"::(_ in ("svd", "lsqr", "eigen"))
    shrinkage::Union{Nothing,String,Float64} = nothing::(_ === nothing || _ == "auto" || 0 < _ < 1)
    priors::Option{AbstractVector}   = nothing
    n_components::Option{Int}        = nothing
    store_covariance::Bool           = false
    tol::Float64                     = 1e-4::(_ > 0)
end
MMI.fitted_params(m::BayesianLDA, (f, _, _)) = (
    coef       = f.coef_,
    intercept  = f.intercept_,
    covariance = m.store_covariance ? f.covariance_ : nothing,
    means      = f.means_,
    priors     = f.priors_,
    scalings   = f.scalings_,
    xbar       = f.xbar_,
    classes    = f.classes_,
    explained_variance_ratio = f.explained_variance_ratio_
    )
meta(BayesianLDA,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Bayesian Linear Discriminant Analysis."
    )

# ============================================================================
const BayesianQDA_ = skda(:QuadraticDiscriminantAnalysis)
@sk_clf mutable struct BayesianQDA <: MMI.Probabilistic
    priors::Option{AbstractVector} = nothing
    reg_param::Float64             = 0.0::(_ ≥ 0)
    store_covariance::Bool         = false
    tol::Float64                   = 1e-4::(_ > 0)
end
MMI.fitted_params(m::BayesianQDA, (f, _, _)) = (
    covariance = m.store_covariance ? f.covariance_ : nothing,
    means      = f.means_,
    priors     = f.priors_,
    rotations  = f.rotations_,
    scalings   = f.scalings_
    )
meta(BayesianQDA,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Bayesian Quadratic Discriminant Analysis."
    )


text:
struct Gtilde{Y<:AbstractVector{<:Real}, A<:Dirichlet,
              B<:Union{Normal, <:OrderedNormalMeanPrior}, C<:Real,
              D<:LogNormal, E<:Normal, F<:Gamma} <: MCMC.Model
  yC::Y  # finite log expressions from control group
  yT::Y  # finite log expressions from treatment group
  K::Int  # number of mixture components
  beta::Bool  # model indicator
  eta::A
  mu::B
  a_omega::C
  nu::D
  psi::E
  tau::F
end

has_ordered_mu(m::Gtilde) = m.mu isa OrderedNormalMeanPrior
uses_skew(m::Gtilde) = (m.psi !== Normal(0, 0))
uses_tdist(m::Gtilde) = !(m.nu.μ > 10000 && m.nu.σ == 0)

function make_ordered_prior(yC::AbstractVector{<:Real}, yT::AbstractVector{<:Real})
  # TODO
  error("Not implemented!")
end

function Gtilde(yC, yT, K, beta; eta=Dirichlet(K, 1/K), mu=nothing,
              a_omega=2.5, nu=LogNormal(2, .5), psi=Normal(0, 3),
              tau=Gamma(0.5, 1), skew::Bool=true, tdist::Bool=true)
  skew || (psi = Normal(0, 0))
  tdist || (nu = LogNormal(20000, 0))

  if mu === nothing
    y = [yC; yT]
    mu = Normal(mean(y), std(y))
  end

  return Gtilde(yC, yT, K, Bool(beta), eta, mu, a_omega, nu, psi, tau)
end

function print_model_info(m::Gtilde)
  println("Model info:")
  for key in fieldnames(Gtilde)
    if key == :yC
      println("NC (finite): ", length(m.yC))
    elseif key == :yT
      println("NT (finite): ", length(m.yT))
    else
      println("$(key): ", getfield(m, key))
    end
  end
end

function MCMC.make_init_state(m::Gtilde)
  NC = length(m.yC)
  NT = length(m.yT)
  etaC = rand(m.eta)
  etaT = rand(m.eta)
  lambdaC = rand(Categorical(etaC), NC)
  lambdaT = rand(Categorical(etaT), NT)
  tau = rand(m.tau)
  mu = if has_ordered_mu(m)
    cumsum(rand(m.mu))
  else
    rand(m.mu, m.K)
  end
  omega = rand(InverseGamma(m.a_omega, tau), m.K)
  nu = rand(m.nu, m.K)
  psi = rand(m.psi, m.K)

  if uses_tdist(m)
    vC = rand.(Gamma.(nu[lambdaC] / 2, 2 ./ nu[lambdaC]))
    vT = rand.(Gamma.(nu[lambdaT] / 2, 2 ./ nu[lambdaT]))
  else
    vC = one.(m.yC)
    vT = one.(m.yT)
  end

  if uses_skew(m)
    zetaC = rand.(truncated.(Normal.(0, 1 ./ sqrt.(vC)), 0, Inf))
    zetaT = rand.(truncated.(Normal.(0, 1 ./ sqrt.(vT)), 0, Inf))
  else
    zetaC = zero.(m.yC)
    zetaT = zero.(m.yT)
  end

  sigma = @.scalefromaltskewt(sqrt(omega), psi)
  phi = @.skewfromaltskewt(sqrt(omega), psi)

  return (etaC=etaC, etaT=etaT, lambdaC=lambdaC, lambdaT=lambdaT, tau=tau,
          mu=mu, omega=omega, nu=nu, psi=psi, vC=vC, vT=vT, zetaC=zetaC,
          zetaT=zetaT, sigma=sigma, phi=phi)
end

Nsum(m::Gtilde) = length(m.yC) + length(m.yT)

function make_sampler(m::Gtilde; init=nothing)
  init === nothing && (init = MCMC.make_init_state(m))

  if uses_tdist(m)
    cond_nu = RWM(:nu, logprob_nu, mvarwm(init.nu), bijector=Bijectors.Log{1}())
    _update_vC = update_vC
    _update_vT = update_vT
  else
    cond_nu = Conditional(:nu, (m, s) -> fill(20000.0, m.K))
    _update_vC = (m, s) -> one.(m.yC)
    _update_vT = (m, s) -> one.(m.yT)
  end

  if uses_skew(m)
    _update_zetaC = update_zetaC
    _update_zetaT = update_zetaT
    _update_psi = update_psi
  else
    _update_zetaC = (m, s) -> zero.(m.yC)
    _update_zetaT = (m, s) -> zero.(m.yT)
    _update_psi = (m, s) -> zeros(m.K)
  end

  if has_ordered_mu(m)
    _update_mu = (m, s) -> let
      lambda = [s.lambdaC; s.lambdaT]
      v = [s.vC; s.vT]
      zeta = [s.zetaC; s.zetaT]
      y = [m.yC; m.yT]

      s2 = s.omega[lambda] ./ v
      g = y - s.psi[lambda] .* zeta
      update_mu(m.mu, s.mu, g, s2, lambda)
    end
  else
    _update_mu = update_mu
  end

  return Gibbs(m,
               Conditional(:etaC, update_etaC),
               Conditional(:etaT, update_etaT),
               Conditional(:lambdaC, update_lambdaC),
               Conditional(:lambdaT, update_lambdaT),
               Conditional(:vC, _update_vC),
               Conditional(:vT, _update_vT),
               Conditional(:zetaC, _update_zetaC),
               Conditional(:zetaT, _update_zetaT),
               Conditional(:mu, _update_mu),
               Conditional(:tau, update_tau),
               Conditional(:omega, update_omega),
               cond_nu,
               Conditional(:psi, _update_psi),
               Conditional((:sigma, :phi), update_sigma_phi))
end


text:
include("parallelNNConv.jl")
include("parallelConvolve.jl")
include("parallelImg2colConvolve.jl")

###convolution layers forprop

@doc raw"""
    function layerForProp(
        cLayer::ConvLayer,
        Ai::AbstractArray = Array{Any,1}(undef,0);
        FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
        kwargs...
    )

Perform the layer forward propagation for a `ConvLayer`

# Arguments

- `cLayer::ConvLayer`

- `Ai` := optional activation of the previous layer

- `FCache` := a `Dict` holds the outputs of `layerForProp` of the previous `Layer`(s)

# Returns

- `Dict(:Z => Z, :A => Ao)`
"""
function layerForProp(
    cLayer::CL,
    Ai::AbstractArray = Array{Any,1}(undef,0);
    FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
    kwargs...
) where {CL <: ConvLayer}

    kwargs = Dict{Symbol, Any}(kwargs...)
    NNlib = getindex(kwargs, :NNlib; default=true)
    img2col = getindex(kwargs, :img2col; default=false)
    if length(Ai) == 0
        Ai = FCache[cLayer.prevLayer][:A]
    end
    cLayer.inputS = size(Ai)
    paddedS = paddedSize(cLayer, Ai)[1:end-2]
    s = cLayer.s
    f = cLayer.f
    outputS = outDims(cLayer, Ai)[1:end-2]
    ci, m = cLayer.inputS[end-1:end]
    co = cLayer.channels

    if NNlib
        D = NNConv(cLayer, Ai)
    elseif img2col
        ## in case input different size than previous time
        if (prod((outputS..., co)), prod((paddedS..., ci))) != size(cLayer.K)
            cLayer.K = unroll(cLayer, (paddedS..., ci, m))
        end
        D = img2colConvolve(cLayer, Ai)
    else
        D = convolve(cLayer, Ai)
    end

    cLayer.outputS = size(D[:A])

    # Ai = nothing
    cLayer.forwCount += 1
    # Base.GC.gc()
    return D

end #function layerForProp(cLayer::Conv1D)

### Pooling Layers

include("parallelFastPooling.jl")
include("parallelPooling.jl")

#import only the needed parts not to have conflict
import NNlib.maxpool, NNlib.meanpool, NNlib.maxpool!, NNlib.meanpool!, NNlib.PoolDims

@doc raw"""
    layerForProp(
        cLayer::PoolLayer},
        Ai::AbstractArray = Array{Any,1}(undef,0);
        FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
        kwargs...
    )

Perform the layer forward propagation for a `PoolLayer`

# Arguments

- `cLayer::PoolLayer`

- `Ai` := optional activation of the previous layer

- `FCache` := a `Dict` holds the outputs of `layerForProp` of the previous `Layer`(s)

# Returns

- `Dict(:A => Ao)`
"""
function layerForProp(
    cLayer::PL,
    Ai::AbstractArray = Array{Any,1}(undef,0);
    FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
    kwargs...
) where {PL <: PoolLayer}

    kwargs = Dict{Symbol, Any}(kwargs...)
    fastPool = getindex(kwargs, :fastPool; default=true)
    NNlib = getindex(kwargs, :NNlib; default=true)
    if length(Ai) == 0
        Ai = FCache[cLayer.prevLayer][:A]
    end

    cLayer.inputS = size(Ai)
    paddedS = paddedSize(cLayer, Ai)[1:end-2]
    s = cLayer.s
    f = cLayer.f
    outputS = outDims(cLayer, Ai)[1:end-2]
    ci, m = cLayer.inputS[end-1:end]
    co = cLayer.channels
    cLayer.inputS = size(Ai)
    padS = paddingSize(cLayer, Ai)

    Ao = zeros(eltype(Ai), outputS..., co, m)

    if NNlib
        pooldims = PoolDims(Ai, f, stride = s, padding = padS)
        if cLayer isa MaxPoolLayer
            maxpool!(Ao, Ai, pooldims)
        elseif cLayer isa AveragePoolLayer
            meanpool!(Ao, Ai, pooldims)
        end #if cLayer isa MaxPoolLayer
    elseif s == f && fastPool #to use the built-in reshape and maximum and mean
        fastPooling!(cLayer, Ai, Ao)
    else
        pooling!(cLayer, Ai, Ao)
    end #if NNlibConv

    cLayer.outputS = size(Ao)
    # Ai = nothing
    cLayer.forwCount += 1
    # Base.GC.gc()
    return Dict(:A => Ao)
end #unction layerForProp(cLayer::OneD) where {OneD <: Union{MaxPool1D, AveragePool1D}}


text:
module TestIterate
using CompressHashDisplace
using Test

@testset "FrozenDict iterate" begin
    dict = Dict("test$i" => i for i = 1:10)
    fdict = FrozenDict(dict)

    @test length(fdict) == length(dict)
    i = 0
    for (k, v) in fdict
        @test v == dict[k]
        i += 1
    end
    @test i == length(fdict)
end
@testset "Nonstring iterate" begin
    dict = Dict(i => 2i for i = 1:10)
    fdict = FrozenDict(dict)

    @test length(fdict) == length(dict)
    i = 0
    for (k, v) in fdict
        @test v == dict[k]
        i += 1
    end
    @test i == length(fdict)
end

end # module