text:
@testset "MLJ model interface" begin
    data = CSV.read(joinpath(pkgdir(Econometrics), "data", "insure.csv"), DataFrame, select = [:insure, :age, :male, :nonwhite, :site]) |>
        dropmissing |>
        (data -> transform!(data, [:insure, :site] .=> categorical, renamecols = false))
    @testset "Econometrics.jl + MLJModelInterface.jl" begin
        @testset "specify contrasts" begin
            contrasts = Dict(:insure => DummyCoding(base = "Uninsure"))
            formula = @formula(insure ~ age + male + nonwhite + site)
            label = :insure
            model = EconometricsMLJModel(;
                model = EconometricModel,
                contrasts = contrasts,
                formula = formula,
                label = label,
                )
            X = data[!, [:age, :male, :nonwhite, :site]]
            y = data[!, :insure]
            verbosity = 1
            fitresult, _, _ = MLJModelInterface.fit(model, verbosity, X, y)
            @test fitresult isa EconometricModel
            Xnew = deepcopy(X)
            @test_broken ynew = MLJModelInterface.predict(model, fitresult, Xnew)
        end
        @testset "w/o specifying contrasts" begin
            formula = @formula(insure ~ age + male + nonwhite + site)
            label = :insure
            model = EconometricsMLJModel(;
                model = EconometricModel,
                formula = formula,
                label = label,
                )
            X = data[!, [:age, :male, :nonwhite, :site]]
            y = data[!, :insure]
            verbosity = 1
            fitresult, _, _ = MLJModelInterface.fit(model, verbosity, X, y)
            @test fitresult isa EconometricModel
            Xnew = deepcopy(X)
            @test_broken ynew = MLJModelInterface.predict(model, fitresult, Xnew)
        end
    end
    @testset "Econometrics.jl + MLJBase.jl" begin
        contrasts = Dict(:insure => DummyCoding(base = "Uninsure"))
        formula = @formula(insure ~ age + male + nonwhite + site)
        label = :insure
        model = EconometricsMLJModel(;
            model = EconometricModel,
            contrasts = contrasts,
            formula = formula,
            label = label,
            )
        X = data[!, [:age, :male, :nonwhite, :site]]
        y = data[!, :insure]
        mach = MLJBase.machine(model, X, y)
        MLJBase.fit!(mach)
        @test mach.fitresult isa EconometricModel
        Xnew = deepcopy(X)
        @test_broken ynew = MLJBase.predict(mach, Xnew)
        mach = MLJBase.machine(model, X, y)
        @test_broken cv_result = MLJBase.evaluate!(mach, resampling=MLJBase.CV(; nfolds = 6, shuffle = true), measure=MLJBase.brier_score, operation=MLJBase.predict)
    end
end


text:
#=
Given a dictionary of words and a string made up of those words (no spaces), return the original sentence in a list. If there is more than one possible reconstruction, return any of them. If there is no possible reconstruction, then return null.

For example, given the set of words 'quick', 'brown', 'the', 'fox', and the string "thequickbrownfox", you should return ['the', 'quick', 'brown', 'fox'].

Given the set of words 'bed', 'bath', 'bedbath', 'and', 'beyond', and the string "bedbathandbeyond", return either ['bed', 'bath', 'and', 'beyond] or ['bedbath', 'and', 'beyond'].
=#

function find_sentence(s::String, dictionary::Array{String, 1})
    sentence = []
    while length(s) != 0
        word = ""
        char_index = 0
        while !(word in dictionary)
            char_index = char_index + 1
            word = word * s[char_index]
        end
        push!(sentence, word)
        s = s[char_index+1:end]
    end
    return sentence
end

ext: jl

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
module SPDot

export spdot

function spdot(v1::SparseMatrixCSC{Float64, Int64},
               v2::SparseMatrixCSC{Float64, Int64})
    sum = 0.0
    if length(v1.rowval) < length(v2.rowval)
        for iv1 in 1:length(v1.rowval)
            sum += v1.nzval[iv1] * v2[v1.rowval[iv1]]
        end
    else
        for iv2 in 1:length(v2.rowval)
            sum += v2.nzval[iv2] * v1[v2.rowval[iv2]]
        end
    end
    return sum
end

function spdot(v1::Vector{Float64},
               v2::SparseMatrixCSC{Float64, Int64})
    sum = 0.0
    for iv2 in 1:length(v2.rowval)
        sum += v2.nzval[iv2] * v1[v2.rowval[iv2]]
    end
    return sum
end

function spdot(m1::Matrix{Float64}, icol::Int64,
               v2::SparseMatrixCSC{Float64, Int64})
    sum = 0.0
    for iv2 in 1:length(v2.rowval)
        sum += v2.nzval[iv2] * m1[v2.rowval[iv2], icol]
    end
    return sum
end

function spdot(v1::SparseMatrixCSC{Float64, Int64},
               v2::Vector{Float64})
    sum = 0.0
    for iv1 in 1:length(v1.rowval)
        sum += v1.nzval[iv1] * v2[v1.rowval[iv1]]
    end
    return sum
end

function spdot(v1::SparseMatrixCSC{Float64, Int64},
               m2::Matrix{Float64}, icol::Int64)
    sum = 0.0
    for iv1 in 1:length(v1.rowval)
        sum += v1.nzval[iv1] * m2[v1.rowval[iv1], icol]
    end
    return sum
end

end

max_stars_repo_name: tsamtsam/SV-RL

max_issues_repo_name: tsamtsam/SV-RL

max_forks_repo_name: tsamtsam/SV-RL


text:
"""
    randompartition(h::Hypergraph, n::Int)::Vector{Set{Int}}

Generates a random partition for vertices of a hypergraph `h` into `n` subsets.
"""
randompartition(h::Hypergraph, n::Int) = randompartition(nhv(h), n)

"""
    randompartition(N::Int, n::Int)::Vector{Set{Int}}

Generates a random partition for graph having `N` vertices into `n` subsets.
"""
function randompartition(N::Int, n::Int)
    res = [Set{Int}() for i in 1:n]
    for i in 1:N
        push!(res[rand(1:n)], i)
    end
    res
end

"""
    HypergraphAggs(h::Hypergraph)

Precomputes vertex and edge basic stats for a hypergraph.
The stats are being used for efficiency reasons by community search algorithms.
"""
struct HypergraphAggs
    hes::Vector{Int}
    max_hes::Int
    Ed::Vector{Int}
    deg_vs::Vector{Int}
    volV::Int
    HypergraphAggs(h::Hypergraph) = begin
        hes = [length(h.he2v[i]) for i in 1:nhe(h)]
        max_hes = maximum(hes)
        Ed = zeros(Int, max_hes)
        for ii in 1:length(hes)
            if hes[ii] >= 1
                Ed[hes[ii]]+=1
            end
        end
        deg_vs = [length(h.v2he[i]) for i in 1:nhv(h)]
        volV = sum(deg_vs)
        new(hes,max_hes,Ed,deg_vs,volV)
    end
end

"""
    LightGraphs.modularity(h::Hypergraph, partition::Vector{Set{Int}},
ha::HypergraphAggs=HypergraphAggs(h))

Calculates the strict modularity of a hypergraph `h` for a given `partition` using
the precomputed aggregates `ha`.
"""
@inline function LightGraphs.modularity(h::Hypergraph, partition::Vector{Set{Int}},
        ha::HypergraphAggs=HypergraphAggs(h))

    @boundscheck sum(length.(partition)) == nhv(h)
    @boundscheck union(partition...) == Set(1:nhv(h))
    volP_volV = [sum(ha.deg_vs[i] for i in p)/ha.volV for p in partition]
    eP = [count(i-> ha.hes[i]>0 && (keys(h.he2v[i]) ⊆ p), 1:nhe(h)) for p in partition]
    (sum(eP) - sum( ha.Ed[d]*sum(volP_volV.^d) for d in 1:ha.max_hes)) / nhe(h)
end

"""
The base type for all algorithms representing various community search patterns.
"""
abstract type AbstractCommunityFinder end

"""
    CFModularityRandom(n::Int, reps::Int) <: AbstractCommunityFinder

Represents a random search over the hypergraph `h` that finds
a partition into `n` communities (subsets) having the maximum modularity value.
During the search `reps` random `n`-partitions will be evaluated.
If there are many partitions having the same value the first
one that was randomly found will be returned.
"""
struct CFModularityRandom <: AbstractCommunityFinder
    n::Int
    reps::Int
end

"""
    findcommunities(h::Hypergraph, method::CFModularityRandom)

Makes a random search over the hypergraph `h` and finds
a partition into `method.n` communities (subsets) having the maximum modularity value.
During the search `method.reps` random `n`-partitions will be evaluated.
If there are many partitions having the same value the first
one that was randomly found will be returned.

Returns a `NamedTuple` where the field `bp` contains partition
and the field `bm` contains the modularity value for that partition.
"""
function findcommunities(h::Hypergraph, method::CFModularityRandom)
    bp = [Int[]]
    bm = -Inf
    ha = HypergraphAggs(h)
    for i in 1:method.reps
        p = Vector{Vector{Int}}([[]])
        while  minimum(length.(p)) == 0
            p = randompartition(h, method.n)
        end
        m = modularity(h, p, ha)
        if m > bm
            bm = m
            bp = p
        end
    end
    (bp=bp, bm=bm)
end

function find_first(c::Array{Set{Int}}, vals)
    for i in 1:length(c)
        for v in vals
            v in c[i] && return i
        end
    end
    throw("None of values in $vals found")
end

"""
    CFModularityCNMLike(n::Int, reps::Int) <: AbstractCommunityFinder

Represents a CNM-Like algorithm for finding communities.
In the algorithm we start with a partition where each node is in its own part.
Then in each step, we randomly select a hyperedge.
Subsequently, we consider merging each set of that parts it touches.
We actually merge the parts if the new best modularity is at least as high
as the modularity from the previous step.
The algortithm iterates through `reps` of repetitions.

For more information see `Algorithm 1` at:
Clustering via Hypergraph Modularity (submitted to Plos ONE), auhtors:
Bogumił Kamiński, Valerie Poulin, Paweł Prałat, Przemysław Szufel, Francois Theberge

"""
struct CFModularityCNMLike <: AbstractCommunityFinder
    reps::Int
end

"""
    findcommunities(h::Hypergraph, method::CFModularityCNMLike)

Iterates a CNM-Like algorithm for finding communities.
In the algorithm we start with a partition where each node is in its own part.
Then in each step, we randomly select a hyperedge.
Subsequently, we consider merging each set of that parts it touches.
We actually merge the parts if the new best modularity is at least as high
as the modularity from the previous step.

Returns a `NamedTuple` where the field `bp` contains partition
and the field `bm` contains the modularity value for that partition,
finally, the fiel `mod_history` represents modularities achieved
in subsequent steps of the algorithm.

For more information see `Algorithm 1` at:
Clustering via Hypergraph Modularity (submitted to Plos ONE), authors:
Bogumił Kamiński, Valerie Poulin, Paweł Prałat, Przemysław Szufel, Francois Theberge

"""
function findcommunities(h::Hypergraph, method::CFModularityCNMLike)
    ha = HypergraphAggs(h)
    best_modularity = 0
    comms = [Set(i) for i in 1:nhv(h)]
    mod_history = Vector{Float64}(undef, method.reps)
    for rep in 1:method.reps
        he = rand(1:nhe(h))
        vers = collect(keys(getvertices(h, he)))
        if length(vers) == 0
            continue
        end;
        c = deepcopy(comms)
        i0 = find_first(c, vers)
        max_i = length(c)
        i_cur = i0
        while i_cur < max_i
            i_cur += 1
            if length(intersect(c[i_cur],vers)) > 0
                union!(c[i0], c[i_cur])
                c[i_cur]=c[max_i]
                max_i += -1
            end
        end
        resize!(c,max_i)
        m = modularity(h, c, ha)
        if m > best_modularity
            best_modularity = m
            comms = c
        end
        mod_history[rep] = best_modularity
    end
    return (bm=best_modularity, bp=comms, mod_history=mod_history)
end


text:
import Base: display, redisplay

using DomTerm: DomTermInline
using Colors: Colorant
using MathMLRepr: LATEX
using Mxma: MExpr

#
#using DataFrames: DataFrame
#display(d::DomTermInline, x::DataFrame) = display(d, MIME("text/html"), x)
#
#using GR: SVG
#display(d::DomTermInline, x::SVG) = display(d, MIME("image/svg+xml"), x)

# Register single-argument function `display(x)` for all types that DomTerm
# supports.
display(d::DomTermInline, x::LATEX) = display(d, MIME("text/mathml+xml"), x)
display(d::DomTermInline, x::MExpr) = display(d, MIME("text/mathml+xml"), x)

typeset_matvec = false
function typeset_matrix(f::Bool)
    global typeset_matvec

    typeset_matvec = f
    return typeset_matvec
end

function display(d::DomTermInline, x::Vector{<:Number})
    if typeset_matvec
        display(d, MIME("text/mathml+xml"), x)
    else
        display(d, MIME("text/plain"), x)
    end
end
function display(d::DomTermInline, x::AbstractMatrix)
    if typeset_matvec
        display(d, MIME("text/mathml+xml"), x)
    else
        display(d, MIME("text/plain"), x)
    end
end

# in-memory Image objects
function display(d::DomTermInline, x::AbstractMatrix{<:Colorant})
    n = length(x)
    # use svg for small images, use png for large images
    if n < 4*1024
        display(d, MIME("image/svg+xml"), x)
    else
        display(d, MIME("image/png"), x)
    end
end
function display(d::DomTermInline, x::Vector{C}) where C<:Colorant
    n = length(x)
    r = Matrix{C}(undef, n, 1)
    for i = 1:n r[i,1] = x[i] end
    display(d, r)
end
function display(d::DomTermInline, x::C) where C<:Colorant
    r = Matrix{C}(undef, 1, 1)
    r[1,1] = x
    display(d, r)
end

const surfacetype_to_mime = Dict(
    :png => MIME("image/png"),
    :svg => MIME("image/svg+xml"),
    :jpg => MIME("image/jpeg"),
    :jpeg => MIME("image/jpeg")
    )

function displayinline(m::Module)
    if Symbol(m) === :Luxor
        @eval begin
            function display(d::DomTermInline, x::$m.Drawing)
                display(d, surfacetype_to_mime[x.surfacetype], x)
            end
        end
    elseif Symbol(m) === :GR
        @eval begin
            display(d::DomTermInline, x::$m.SVG) = display(d, MIME("image/svg+xml"), x)
            $m.inline()
        end
    elseif Symbol(m) === :DataFrames
        @eval begin
            display(d::DomTermInline, x::$m.DataFrame) = display(d, MIME("text/html"), x)
        end
    else
        @error("inline display for module: ", Symbol(m), " is not supported.")
    end
end
#


text:
@testset "cat-light" begin
    setlight()
    x = 1:5
    @test_throws M.InterfaceError M.categorical(x)
end
@testset "cat-full" begin
    setfull()
    M.categorical(::FI, a...; kw...) = categorical(a...; kw...)
    x = 1:5
    @test M.categorical(x) == categorical(x)
end
# ------------------------------------------------------------------------
@testset "matrix-light" begin
    setlight()
    X = ones(2,3)
    @test matrix(X) === X
    @test matrix(X, transpose=true) == ones(3,2)
    X = (1,2,3)
    @test_throws ArgumentError matrix(X)
    X = (a=[1,2,3],b=[1,2,3])
    @test_throws M.InterfaceError matrix(X)
end
@testset "matrix-full" begin
    setfull()
    M.matrix(::FI, ::Val{:table}, X; kw...) = Tables.matrix(X; kw...)
    X = (a=[1,2,3],b=[1,2,3])
    @test matrix(X) == hcat([1,2,3],[1,2,3])
end
# ------------------------------------------------------------------------
@testset "int-light" begin
    setlight()
    x = categorical([1,2,3])
    @test_throws M.InterfaceError int(x)
end
@testset "int-full" begin
    setfull()
    M.int(::FI, x::CategoricalValue; kw...) =
        collect(1:length(levels(x.pool)))[x.level]
    x = categorical(['a','b','a'])
    @test int(x[1]) == 0x01
    @test int(x[2]) == 0x02
    @test_broken int(x[2]) isa UInt32
    @test int(x[1], type=Int64) == 1
    @test int(x[1], type=Int64) isa Int64
end
# ------------------------------------------------------------------------
@testset "classes-light" begin
    setlight()
    x = categorical(['a','b','a'])
    @test_throws M.InterfaceError classes(x)
end
@testset "classes-full" begin
    setfull()
    M.classes(::FI, p::CategoricalPool) =
        [p[i] for i in invperm(1:length(levels(p)))]
    M.classes(::FI, x::CategoricalValue) = classes(x.pool)
    x = categorical(['a','b','a'])
    @test classes(x[1]) == ['a', 'b']
end
# ------------------------------------------------------------------------
@testset "schema-light" begin
    # throw error for any input anyway
    setlight()
    ary = rand(10, 3)
    @test_throws M.InterfaceError M.schema(ary)
    df = DataFrame(rand(10, 3))
    @test_throws M.InterfaceError M.schema(df)
end
@testset "schema-full" begin
    setfull()
    ary = rand(10, 3)
    @test M.schema(ary) === nothing
    M.schema(::FI, ::Val{:table}, X; kw...) =
        MLJScientificTypes.schema(X; kw...)
    df = DataFrame(A = rand(10), B = categorical(rand('a':'c', 10)))
    sch = M.schema(df)
    @test sch.names == (:A, :B)
    @test sch.types[1] <: Float64
    @test sch.types[2] <: CategoricalValue
    @test sch.scitypes[1] <: Continuous
    @test sch.scitypes[2] <: Multiclass
end
# ------------------------------------------------------------------------
@testset "istable" begin
    setlight()
    X = rand(5)
    @test !M.istable(X)
    X = randn(5,5)
    @test !M.istable(X)
    X = DataFrame(A=rand(10))
    @test M.istable(X)
end
# ------------------------------------------------------------------------
@testset "decoder-light" begin
    setlight()
    x = 5
    @test_throws M.InterfaceError decoder(x)
end
@testset "decoder-full" begin
    setfull()
    # toy test because I don't want to copy the decoder logic here
    M.decoder(::FI, x) = 0
    @test decoder(nothing) == 0
end
# ------------------------------------------------------------------------
@testset "table-light" begin
    setlight()
    X = ones(3,2)
    @test_throws M.InterfaceError table(X)
end
@testset "table-full" begin
    setfull()
    function M.table(::FI, A::AbstractMatrix; names=nothing)
        _names = [Symbol(:x, j) for j in 1:size(A, 2)]
        return Tables.table(A, header=_names)
    end
    X = ones(3,2)
    T = table(X)
    @test Tables.istable(T)
    @test Tables.matrix(T) == X
end
# ------------------------------------------------------------------------
@testset "nrows-light" begin
    setlight()
    X = ones(5)
    @test nrows(X) == 5
    X = ones(5,3)
    @test nrows(X) == 5
    X = ones(5,3,2)
    @test_throws ArgumentError nrows(X)
    X = (a=[4,2,1],b=[3,2,1])
    @test_throws M.InterfaceError nrows(X)
    @test nrows(nothing) == 0
end
@testset "nrows-full" begin
    setfull()
    M.nrows(::FI, ::Val{:table}, X) = Tables.rowcount(X)
    X = (a=[4,2,1],b=[3,2,1])
    @test nrows(X) == 3
end
# ------------------------------------------------------------------------
@testset "select-light" begin
    setlight()
    X = nothing
    @test selectrows(X, 1) === nothing
    @test selectcols(X, 1) === nothing
    @test select(X, 1, 2)  === nothing

    # vector
    X = ones(5)
    @test selectrows(X, 1)   == [1.0]
    @test selectrows(X, 1:2) == ones(2,)
    @test selectrows(X, :)  === X
    @test_throws ArgumentError selectcols(X, 5)
    @test_throws ArgumentError select(X, 2, 2)

    # matrix
    X = ones(5, 3)
    @test selectrows(X, 1)    == ones(1, 3)
    @test selectrows(X, 1:2)  == ones(2, 3)
    @test selectrows(X, :)   === X
    @test selectcols(X, 1)    == ones(5,)
    @test selectcols(X, 1:2)  == ones(5, 2)
    @test selectcols(X, :)   === X
    @test select(X, 1, 1)     == 1.0
    @test select(X, 1:2, 1)   == ones(2,)
    @test select(X, 1:2, 1:2) == ones(2, 2)

    # table
    X = (x=[1,1,1],y=[2,2,2])
    @test_throws M.InterfaceError selectrows(X, 1)
    @test_throws M.InterfaceError selectcols(X, 1)

    # something else
    X = (1,2,3)
    @test_throws ArgumentError selectrows(X, 1)
    @test_throws ArgumentError selectcols(X, 1)
    @test_throws ArgumentError select(X, 1, 1)
end
@testset "select-full" begin
    setfull()
    M.selectrows(::FI, ::Val{:table}, X, ::Colon) = X
    M.selectcols(::FI, ::Val{:table}, X, ::Colon) = X
    function M.selectrows(::FI, ::Val{:table}, X, r)
        r = r isa Integer ? (r:r) : r
        cols = Tables.columntable(X)
        new_cols = NamedTuple{keys(cols)}(tuple((c[r] for c in values(cols))...))
        return Tables.materializer(X)(new_cols)
    end
    function M.selectcols(::FI, ::Val{:table}, X, c::Union{Symbol,Integer})
        cols = Tables.columntable(X) # named tuple of vectors
        return cols[c]
    end
    function M.selectcols(::FI, ::Val{:table}, X, c::AbstractArray)
        cols = Tables.columntable(X) # named tuple of vectors
        newcols = project(cols, c)
        return Tables.materializer(X)(newcols)
    end
    # project named tuple onto a tuple with only specified `labels` or indices:
    project(t::NamedTuple, labels::AbstractArray{Symbol}) =
        NamedTuple{tuple(labels...)}(t)
    project(t::NamedTuple, label::Colon) = t
    project(t::NamedTuple, label::Symbol) = project(t, [label,])
    project(t::NamedTuple, indices::AbstractArray{<:Integer}) =
        NamedTuple{tuple(keys(t)[indices]...)}(tuple([t[i] for i in indices]...))
    project(t::NamedTuple, i::Integer) = project(t, [i,])

    X = (x=[1,2,3],y=[4,5,6], z=[0,0,0])
    @test selectrows(X, 1)   == (x=[1],y=[4],z=[0])
    @test selectrows(X, 1:2) == (x=[1,2],y=[4,5],z=[0,0])
    @test selectrows(X, :)  === X
    @test selectcols(X, 1)   == [1,2,3]
    @test selectcols(X, 1:2) == (x = [1, 2, 3], y = [4, 5, 6])
    @test selectcols(X, :)  === X
    @test select(X, 1, 1)    == 1
    @test select(X, 1:2, 1)  == [1,2]
    @test select(X, :, 1)    == [1,2,3]
    @test selectcols(X, :x)  == [1,2,3]
    @test select(X, 1:2, :z) == [0,0]
    #
    # extra tests by Anthony
    X = (x=[1,2,3], y=[10, 20, 30], z= [:a, :b, :c])
    @test select(X, 2, :y)         == 20
    @test select(X, 2, [:x, :y])   == (x=[2,], y=[20,])
    @test select(X, 2:3, :x)       == [2, 3]
    @test select(X, 2:3, [:x, :y]) == (x=[2, 3], y=[20, 30])
    @test select(X, :, [:x, :y])   == select(X, 1:3, [:x, :y])
    @test select(X, 2, :)          == select(X, 2, 1:3)
    @test select(X, 2:3, :)        == select(X, 2:3, 1:3)

    @test select(X, 2, 2)        == 20
    @test select(X, 2, [1, 2])   == (x=[2,], y=[20,])
    @test select(X, 2:3, 1)      == [2, 3]
    @test select(X, 2:3, [1, 2]) == (x=[2, 3], y=[20, 30])
    @test select(X, :, [1, 2])   == select(X, 1:3, [1, 2])
    @test select(X, 2, :)        == select(X, 2, 1:3)
    @test select(X, 2:3, :)      == select(X, 2:3, 1:3)
end
# ------------------------------------------------------------------------
@testset "univ-finite" begin
    setlight()
    @test_throws M.InterfaceError UnivariateFinite(Dict(2=>3,3=>4))
    @test_throws M.InterfaceError UnivariateFinite(randn(2), randn(2))
end


text:
function render!(
    rendermodel::RenderModel,
    boundary::LaneBoundary,
    pts::Matrix{Float64},
    lane_marking_width  :: Real=0.15, # [m]
    lane_dash_len       :: Real=0.91, # [m]
    lane_dash_spacing   :: Real=2.74, # [m]
    lane_dash_offset    :: Real=0.00  # [m]
    )

    marker_color = boundary.color == :yellow ? COLOR_LANE_MARKINGS_YELLOW : COLOR_LANE_MARKINGS_WHITE
    if boundary.style == :broken
        add_instruction!(rendermodel, render_dashed_line, (pts, marker_color, lane_marking_width, lane_dash_len, lane_dash_spacing, lane_dash_offset))
    else
        add_instruction!(rendermodel, render_line, (pts, marker_color, lane_marking_width))
    end
    rendermodel
end

function render!(rendermodel::RenderModel, lane::Lane, roadway::Roadway;
    color_asphalt       :: Colorant=COLOR_ASPHALT,
    )

    n = length(lane.curve)
    pts = Array{Float64}(2, n + has_next(lane))
    for (i,pt) in enumerate(lane.curve)
        pts[1,i] = pt.pos.x
        pts[2,i] = pt.pos.y
    end
    if has_next(lane)
        pt = next_lane_point(lane, roadway)
        pts[1,end] = pt.pos.x
        pts[2,end] = pt.pos.y
    end

    add_instruction!(rendermodel, render_line, (pts, color_asphalt, lane.width))
    rendermodel
end

function render!(rendermodel::RenderModel, roadway::Roadway;
    color_asphalt       :: Colorant=COLOR_ASPHALT,
    lane_marking_width  :: Real=0.15, # [m]
    lane_dash_len       :: Real=0.91, # [m]
    lane_dash_spacing   :: Real=2.74, # [m]
    lane_dash_offset    :: Real=0.00  # [m]
    )

    # render the asphalt between the leftmost and rightmost lane markers
    for seg in roadway.segments
        if !isempty(seg.lanes)
            laneR = seg.lanes[1]
            laneL = seg.lanes[end]

            pts = Array{Float64}(2, length(laneL.curve) + has_next(laneL) +
                                        length(laneR.curve) + has_next(laneR) +
                                        2*length(seg.lanes))
            pts_index = 0
            for pt in laneL.curve
                edgept = pt.pos + polar(laneL.width/2, pt.pos.θ + π/2)
                pts_index += 1
                pts[1, pts_index] = edgept.x
                pts[2, pts_index] = edgept.y
            end
            if has_next(laneL)
                pt = next_lane_point(laneL, roadway)
                edgept = pt.pos + polar(laneL.width/2, pt.pos.θ + π/2)
                pts_index += 1
                pts[1, pts_index] = edgept.x
                pts[2, pts_index] = edgept.y
            end
            for i in reverse(1:length(seg.lanes))
                lane = seg.lanes[i]
                if has_next(lane)
                    pt = next_lane_point(lane, roadway).pos
                else
                    pt = lane.curve[end].pos
                end
                pts_index += 1
                pts[1, pts_index] = pt.x
                pts[2, pts_index] = pt.y
            end

            if has_next(laneR)
                pt = next_lane_point(laneR, roadway)
                edgept = pt.pos + polar(laneR.width/2, pt.pos.θ - π/2)
                pts_index += 1
                pts[1, pts_index] = edgept.x
                pts[2, pts_index] = edgept.y
            end
            for j in length(laneR.curve) : -1 : 1
                pt = laneR.curve[j]
                edgept = pt.pos + polar(laneR.width/2, pt.pos.θ - π/2)
                pts_index += 1
                pts[1, pts_index] = edgept.x
                pts[2, pts_index] = edgept.y
            end
            for i in 1:length(seg.lanes)
                lane = seg.lanes[i]
                pt = lane.curve[1].pos
                pts_index += 1
                pts[1, pts_index] = pt.x
                pts[2, pts_index] = pt.y
            end

            add_instruction!(rendermodel, render_fill_region, (pts, color_asphalt))
        end
        # for lane in seg.lanes
        #     render!(rendermodel, lane, roadway)
        # end
    end

    # render the lane edges
    for seg in roadway.segments
        for lane in seg.lanes

            N = length(lane.curve)
            halfwidth = lane.width/2

            # always render the left lane marking
            pts_left = Array{Float64}(2, N)
            for (i,pt) in enumerate(lane.curve)
                p_left = pt.pos + polar(halfwidth, pt.pos.θ + π/2)

                pts_left[1,i] = p_left.x
                pts_left[2,i] = p_left.y
            end
            if has_next(lane)
                lane2 = next_lane(lane, roadway)
                pt = lane2.curve[1]
                p_left = pt.pos + polar(lane2.width/2, pt.pos.θ + π/2)
                pts_left = hcat(pts_left, [p_left.x, p_left.y])
            end

            render!(rendermodel, lane.boundary_left, pts_left, lane_marking_width, lane_dash_len, lane_dash_spacing, lane_dash_offset)

            # only render the right lane marking if this is the first lane
            if lane.tag.lane == 1
                pts_right = Array{Float64}(2, N)

                for (i,pt) in enumerate(lane.curve)
                    p_right = pt.pos - polar(halfwidth, pt.pos.θ + π/2)

                    pts_right[1,i] = p_right.x
                    pts_right[2,i] = p_right.y
                end

                if has_next(lane)
                    lane2 = next_lane(lane, roadway)
                    pt = lane2.curve[1]
                    p_right = pt.pos - polar(lane2.width/2, pt.pos.θ + π/2)
                    pts_right = hcat(pts_right, [p_right.x, p_right.y])
                end

                render!(rendermodel, lane.boundary_right, pts_right, lane_marking_width, lane_dash_len, lane_dash_spacing, lane_dash_offset)
            end
        end
    end

    rendermodel
end

function render(roadway::Roadway;
    canvas_width::Int=DEFAULT_CANVAS_WIDTH,
    canvas_height::Int=DEFAULT_CANVAS_HEIGHT,
    rendermodel = RenderModel(),
    cam::Camera = FitToContentCamera(),
    )

    s = CairoRGBSurface(canvas_width, canvas_height)
    ctx = creategc(s)
    clear_setup!(rendermodel)

    render!(rendermodel, roadway)

    camera_set!(rendermodel, cam, canvas_width, canvas_height)
    render(rendermodel, ctx, canvas_width, canvas_height)
    s
end

Base.show(io::IO, ::MIME"image/png", roadway::Roadway) = show(io, MIME"image/png"(), render(roadway))


text:
import Base: size, getindex, setindex!, iterate, length, eachindex, CartesianIndices, tail, IndexStyle, copyto!
using TupleTools

function symarrlength(Nts,Nsyms)
    len = 1
    for (Nt,Nsym) in zip(Nts,Nsyms)
        len *= prod(Nt+ii for ii=0:Nsym-1) ÷ factorial(Nsym)
    end
    len
end

struct SymArray{Nsyms,T,N,M,VecType<:AbstractVector} <: AbstractArray{T,N}
    data::VecType
    size::NTuple{N,Int}
    Nts::NTuple{M,Int}
    function SymArray{Nsyms,T}(size::Vararg{Int,N}) where {Nsyms,T,N}
        @assert sum(Nsyms)==N
        ii::Int = 0
        f = i -> (ii+=Nsyms[i]; size[ii])
        Nts = ntuple(f,length(Nsyms))
        len = symarrlength(Nts,Nsyms)
        data = Vector{T}(undef,len)
        new{Nsyms,T,N,length(Nsyms),Vector{T}}(data,size,Nts)
    end
    # this creates a SymArray that serves as a view on an existing vector
    function SymArray{Nsyms}(data::VecType,size::Vararg{Int,N}) where {Nsyms,N,VecType<:AbstractVector{T}} where T
        @assert sum(Nsyms)==N
        ii::Int = 0
        f = i -> (ii+=Nsyms[i]; size[ii])
        Nts = ntuple(f,length(Nsyms))
        len = symarrlength(Nts,Nsyms)
        @assert length(data) == len
        new{Nsyms,T,N,length(Nsyms),VecType}(data,size,Nts)
    end
end

size(A::SymArray) = A.size
length(A::SymArray) = length(A.data)

copyto!(S::SymArray,A::AbstractArray) = begin
    @assert size(S) == size(A)
    @inbounds for (i,I) in enumerate(eachindex(S))
        S[i] = A[I]
    end
    S
end

SymArray{Nsyms}(A::AbstractArray{T}) where {Nsyms,T} = (S = SymArray{Nsyms,T}(size(A)...); copyto!(S,A))
# to avoid ambiguity with Vararg "view" constructor above
SymArray{(1,)}(A::AbstractVector{T}) where {T} = (S = SymArray{(1,),T}(size(A)...); copyto!(S,A))

"`SymArr_ifsym(A,Nsyms)` make a SymArray if there is some symmetry (i.e., any of the Nsyms are not 1)"
SymArr_ifsym(A,Nsyms) = all(Nsyms.==1) ? A : SymArray{(Nsyms...,)}(A)

@generated sub2ind(A::SymArray{Nsyms,T,N}, I::Vararg{Int,N}) where {Nsyms,T,N} = begin
    body = quote
        stride::Int = 1
        ind::Int = 1
        Nt::Int = 0
    end
    ii::Int = 0
    for (iN,Nsym) in enumerate(Nsyms)
        isyms = Symbol.(:i,1:Nsym)
        Ilocs = ( :( I[$(ii+=1)] ) for _=1:Nsym)
        indexpr = :( $(isyms[1])-1 )
        for (inum,isym) in enumerate(isyms[2:end])
            prodterms = (:( $isym + $(j-1) ) for j=0:inum)
            indexpr = :( $indexpr + *($(prodterms...)) ÷ $(factorial(inum+1)) )
        end
        strideterms = ( :( Nt+$ii ) for ii=0:Nsym-1)
        expr = quote
            Nt = A.Nts[$iN]
            ($(isyms...),) = TupleTools.sort(($(Ilocs...),))
            ind += $indexpr * stride
            stride *= *($(strideterms...)) ÷ $(factorial(Nsym))
        end
        push!(body.args,expr)
    end
    push!(body.args,:( ind ))
    body
end

IndexStyle(::Type{<:SymArray}) = IndexCartesian()
getindex(A::SymArray, i::Int) = A.data[i]
getindex(A::SymArray{Nsyms,T,N}, I::Vararg{Int,N}) where {Nsyms,T,N} = A.data[sub2ind(A,I...)]
setindex!(A::SymArray, v, i::Int) = A.data[i] = v
setindex!(A::SymArray{Nsyms,T,N}, v, I::Vararg{Int,N}) where {Nsyms,T,N} = (A.data[sub2ind(A,I...)] = v);

@generated function lessnexts(::SymArray{Nsyms}) where Nsyms
    lessnext = ones(Bool,sum(Nsyms))
    istart = 1
    for Nsym in Nsyms
        istart += Nsym
        lessnext[istart-1] = false
    end
    Tuple(lessnext)
end

struct SymArrayIter{N}
    lessnext::NTuple{N,Bool}
    sizes::NTuple{N,Int}
    "create an iterator that gives i1<=i2<=i3 etc for each index group"
    SymArrayIter(A::SymArray{Nsyms,T,N,M}) where {Nsyms,T,N,M} = new{N}(lessnexts(A),A.size)
end
eachindex(A::SymArray) = CartesianIndices(A);
CartesianIndices(A::SymArray) = SymArrayIter(A);

first(iter::SymArrayIter) = CartesianIndex(map(one, iter.sizes))

@inline function iterate(iter::SymArrayIter)
    iterfirst = first(iter)
    iterfirst, iterfirst
end
@inline function iterate(iter::SymArrayIter, state)
    valid, I = __inc(state.I, iter.sizes, iter.lessnext)
    valid || return nothing
    return CartesianIndex(I...), CartesianIndex(I...)
end
# increment post check to avoid integer overflow
@inline __inc(::Tuple{}, ::Tuple{}, ::Tuple{}) = false, ()
@inline function __inc(state::Tuple{Int}, size::Tuple{Int}, lessnext::Tuple{Int})
    valid = state[1] < size[1]
    return valid, (state[1]+1,)
end

@inline function __inc(state, sizes, lessnext)
    smax = lessnext[1] ? state[2] : sizes[1]
    if state[1] < smax
        return true, (state[1]+1, tail(state)...)
    end
    valid, I = __inc(tail(state), tail(sizes), tail(lessnext))
    return valid, (1, I...)
end


text:
export pooling

const MAXPOOLING2D_F32 = Libdl.dlsym(libmerlin, :maxpooling2d_f32)
const AVGPOOLING2D_F32 = Libdl.dlsym(libmerlin, :avgpooling2d_f32)
const ∇MAXPOOLING2D_F32 = Libdl.dlsym(libmerlin, :maxpooling2d_f32_grad)
const ∇AVGPOOLING2D_F32 = Libdl.dlsym(libmerlin, :avgpooling2d_f32_grad)

maxpooling2d_handle(::Type{Float32}) = MAXPOOLING2D_F32
avgpooling2d_handle(::Type{Float32}) = AVGPOOLING2D_F32
∇maxpooling2d_handle(::Type{Float32}) = ∇MAXPOOLING2D_F32
∇avgpooling2d_handle(::Type{Float32}) = ∇AVGPOOLING2D_F32

"""
    pooling(mode, x::Var, dims::Tuple, [pads], [stride])

N-dimensional pooling function.

* mode: `:max` or `:average`
* dims: kernel size.
* [pads=(0,0,...)]: spatial padding size.
* [strides=(1,1,...)]: kernel strides.

```julia
x = Var(rand(Float32,5,4,3,2))
y = pooling(:max, x, (2,2), pads=(0,0), strides=(1,1))
```
"""
function pooling{N}(mode::Symbol, x::Var, dims::NTuple{N,Int}; pads=nothing, strides=nothing)
    pads == nothing && (pads = ntuple(_ -> 0, N))
    strides == nothing && (strides = ntuple(_ -> 1, N))
    if mode == :max
        f = maxpooling
    elseif mode == :average
        f = avgpooling
    else
        throw("Invalid mode is specified: $(mode).")
    end
    isvoid(x.data) && return Var(nothing, f, (x,dims,pads,strides))
    iscuda(x.data) && return CUDA.f(x, dims, pads, strides)
    f(x, dims, pads, strides)
end

function maxpooling(x::Var, dims::NTuple{2,Int}, pads, strides)
    T = eltype(x.data)
    h = maxpooling2d_handle(T)
    spdims = ntuple(i -> (size(x.data,i) + 2pads[i] - dims[i]) ÷ strides[i] + 1, 2)
    y = Array(T, spdims..., size(x,3), size(x,4))
    inds = Array(Cint, length(y))
    function f{T}(x::Array{T})
        ccall(h, Void, (Ptr{T},Ptr{T},Ptr{Cint},Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint),
            x, y, inds, size(x,1), size(x,2), size(x,3)*size(x,4),
            dims[1], dims[2], pads[1], pads[2], strides[1], strides[2])
    end
    f(x.data)

    function df{T}(gy::Array{T})
        isvoid(x.grad) && return
        ccall(∇maxpooling2d_handle(T), Void, (Ptr{T},Ptr{T},Ptr{Cint},Cint),
            gy, x.grad, inds, length(gy))
    end
    Var(y, df, (x,))
end

function avgpooling(x::Var, dims::NTuple{2,Int}, pads, strides)
    T = eltype(x.data)
    h = avgpooling2d_handle(T)
    spdims = ntuple(i -> (size(x.data,i) + 2pads[i] - dims[i]) ÷ strides[i] + 1, 2)
    y = Array(T, spdims..., size(x,3), size(x,4))
    function f{T}(x::Array{T})
        ccall(h, Void, (Ptr{T},Ptr{T},Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint),
            x, y, size(x,1), size(x,2), size(x,3)*size(x,4),
            dims[1], dims[2], pads[1], pads[2], strides[1], strides[2])
    end
    f(x.data)

    function df{T}(gy::Array{T})
        isvoid(x.grad) && return
        gx = x.grad
        h = ∇avgpooling2d_handle(T)
        ccall(h, Void, (Ptr{T},Ptr{T},Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint,Cint),
            gy, gx, size(gx,1), size(gx,2), size(gx,3)*size(gx,4),
            dims[1], dims[2], pads[1], pads[2], strides[1], strides[2])
    end
    Var(y, df, (x,))
end


text:
abstract type AbstractTransitionSystem end

include("systems/transitions.jl")

include("systems/common.jl")
include("systems/arc_standard.jl")
include("systems/arc_eager.jl")
include("systems/arc_hybrid.jl")
include("systems/arc_swift.jl")
include("systems/listbased.jl")


text:
# ================================================================================
# IO
coefnames(model::MixedModel, anova::Val{:anova}) = vectorize(coefnames(model.formula.rhs[1], anova))

# anovatable api
function _anovatable(aov::AnovaResult{<: LinearMixedModel, FTest}; kwargs...)
    AnovaTable(hcat(vectorize.((
                    dof(aov),
                    dof_residual(aov),
                    teststat(aov),
                    pval(aov)
                    ))...),
              ["DOF", "Res.DOF", "F value", "Pr(>|F|)"],
              ["x$i" for i in eachindex(pval(aov))], 4, 3)
end

function _anovatable(aov::AnovaResult{<: Tuple, LRT},
                    modeltype1::Union{Type{<: LinearMixedModel}, Type{<: TableRegressionModel{<: GLM.LinearModel}}},
                    modeltype2::Type{<: LinearMixedModel};
                    kwargs...)
    if last(aov.model).optsum.REML
        AnovaTable(hcat(vectorize.((
                        dof(aov),
                        [NaN, _diff(dof(aov))...],
                        dof_residual(aov),
                        deviance(aov),
                        teststat(aov),
                        pval(aov)
                        ))...),
                ["DOF", "ΔDOF", "Res.DOF", "-2 logLik", "χ²", "Pr(>|χ²|)"],
                ["$i" for i in eachindex(pval(aov))], 6, 5)
    else
        AnovaTable(hcat(vectorize.((
                        dof(aov),
                        [NaN, _diff(dof(aov))...],
                        dof_residual(aov),
                        aic.(aov.model),
                        bic.(aov.model),
                        deviance(aov),
                        teststat(aov),
                        pval(aov)
                        ))...),
                ["DOF", "ΔDOF", "Res.DOF", "AIC", "BIC", "-2 logLik", "χ²", "Pr(>|χ²|)"],
                ["$i" for i in eachindex(pval(aov))], 8, 7)
    end
end


text:
using StatsBase

"""
    AbstractPolicy
"""
abstract type AbstractPolicy end

# function Base.get(π::AbstractPolicy, state_t, action_t)
#     throw(DomainError("get(PolicyType, args...) not defined!"))
# end

function Base.get(π::AbstractPolicy; kwargs...)
    throw(DomainError("get(PolicyType, args...) not defined!"))
end

"""
    NullPolicy <: AbstractPolicy

Always returns 1.0
"""
struct NullPolicy <: AbstractPolicy
end
# Base.get(π::NullPolicy, state_t, action_t) = 1.0
Base.get(π::NullPolicy; kwargs...) = 1.0

"""
    PersistentPolicy <: AbstractPolicy

Returns 1 of `action_t == action`
"""
struct PersistentPolicy <: AbstractPolicy
    action::Int64
end

# Base.get(π::PersistentPolicy, state_t, action_t) =
#     π.action == action_t ? 1 : 0

Base.get(π::PersistentPolicy; kwargs...) =
        π.action == kwargs[:action_t] ? 1 : 0

"""
    RandomPolicy{T<:AbstractFloat} <: AbstractPolicy

"""
struct RandomPolicy{T<:AbstractFloat} <: AbstractPolicy
    probabilities::Array{T,1}
    weight_vec::Weights{T, T, Array{T, 1}}
    RandomPolicy(probabilities::Array{T,1}) where {T<:AbstractFloat} = new{T}(probabilities, Weights(probabilities))
end

# Base.get(π::RandomPolicy, state_t, action_t) =
#     π.probabilities[action_t]
Base.get(π::RandomPolicy; kwargs...) =
        π.probabilities[kwargs[:action_t]]

StatsBase.sample(π::RandomPolicy) = StatsBase.sample(Random.GLOBAL_RNG, π)
StatsBase.sample(rng::Random.AbstractRNG, π::RandomPolicy) = StatsBase.sample(rng, π.weight_vec)
StatsBase.sample(rng::Random.AbstractRNG, π::RandomPolicy, state) = StatsBase.sample(rng, π.weight_vec)

"""
    FunctionalPolicy{F <: Function} <: AbstractPolicy

"""
struct FunctionalPolicy{F <: Function} <: AbstractPolicy
    func::F
end

# Base.get(π::FunctionalPolicy, state_t, action_t) =
#     π.func(state_t, action_t)
Base.get(π::FunctionalPolicy; kwargs...) =
        π.func(;kwargs...)

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
struct Cycle{T<:Integer}
    n::T
end
nvertices(g::Cycle) = getfield(g, :n)

function AdjMap(g::Cycle{T}) where {T <: Integer}
    n = Int(nvertices(g))
    if n <= 1
        return AdjMap{T}(n)
    elseif n == 2
        g = AdjMap{T}(n)
        add_edge!(g, 1, 2)
        return g
    else
        NP = n + 1
        NM = n - 1
        cptr = Vector{T}(undef, NP)
        @inbounds(setindex!(cptr, 1, 1))
        @inbounds(setindex!(cptr, NP, NP))
        @inbounds for i in OneTo(NM)
            cptr[i + 1] = i
        end
        rval = Vector{T}(undef, n)
        @inbounds(setindex!(rval, 1, NM))
        @inbounds(setindex!(rval, NM, n))
        @inbounds for i in OneTo(n - 2)
            rval[i] = i
        end
        return AdjMap(n, cptr, rval, no_metadata)
    end
end

gs = SimpleGraph(5)
g = AdjMap(5)
add_edge!(g, 1, 2)
add_edge!(g, 1, 5)
add_edge!(g, 2, 3)
add_edge!(g, 3, 4)
add_edge!(g, 4, 5)

julia> add_edge!(g, 1, 2, 1)
true

julia> add_edge!(g, 1, 5, 1)
true

julia> add_edge!(g, 2, 3, 1)
true

julia> add_edge!(g, 3, 4, 1)
true

julia> add_edge!(g, 4, 5, 1)
true

julia> g
{5, 5} undirected simple Int64 graph with Float64 weights

julia>

julia> g.weights
5×5 SparseArrays.SparseMatrixCSC{Float64, Int64} with 10 stored entries:
  ⋅   1.0   ⋅    ⋅   1.0
 1.0   ⋅   1.0   ⋅    ⋅
  ⋅   1.0   ⋅   1.0   ⋅
  ⋅    ⋅   1.0   ⋅   1.0
 1.0   ⋅    ⋅   1.0   ⋅

julia> g.^C

julia> w = weights(g)
5×5 SparseArrays.SparseMatrixCSC{Float64, Int64} with 10 stored entries:
  ⋅   1.0   ⋅    ⋅   1.0
 1.0   ⋅   1.0   ⋅    ⋅
  ⋅   1.0   ⋅   1.0   ⋅
  ⋅    ⋅   1.0   ⋅   1.0
 1.0   ⋅    ⋅   1.0   ⋅

julia> w.colptr
6-element Vector{Int64}:
  1
  3
  5
  7
  9
 11

julia> w.rowval
10-element Vector{Int64}:
 2
 5
 1
 3
 2
 4
 3
 5
 1
 4

julia> outneighbors(g, 2)
2-element view(::Vector{Int64}, 3:4) with eltype Int64:
 1
 3

julia> inneighbors(g, 2)
2-element view(::Vector{Int64}, 3:4) with eltype Int64:
 1
 3


text:
mean(d::LogitNormal; kwargs...) = estimateMean(d,kwargs...) # moved to NormalTransforms

function mode(d::LogitNormal)
    (μ, σ) = params(d)
    # if mu<0 then maximum is left of median, if mu>0 right of median
    # if mu=0 the maximum is either at mu or there are two maxima of the same height
    # here, report the left mode
    interval = μ <= 0.0 ? (0.0,logistic(μ)) : (logistic(μ),1.0)
    resOpt = optimize(x -> -pdf(d, x), interval[1], interval[2])
    Optim.minimizer(resOpt)
end

function var(d::LogitNormal; mean=missing, kwargs...)
    # estimateVariance does not pass kwargs to mean, need to do it here
    m = ismissing(mean) ? estimateMean(d,kwargs...) : mean
    estimateVariance(d; kwargs..., mean=m)
end

### Helpers

# opjective function minimized in mode(d)
# function objectiveModeLogitNormal(mode::Real, mu::Real, sd2::Real)
#     mode <= 0.0 || mode >= 1.0 && return prevfloat(T(Inf))
#     diff =  mu - sd2*(1 - 2*mode) - logit(mode)
#     diff^2
# end


text:
using DelayDiffEq, DiffEqProblemLibrary.DDEProblemLibrary
DDEProblemLibrary.importddeproblems()

# load QS problem
using DDEProblemLibrary: prob_dde_qs

using PGFPlotsX

function plot_qs_comparison()
    # time points at which solution is computed and plotted
    ts = range(prob_dde_qs.tspan[1]; stop = prob_dde_qs.tspan[2], length = 10_000)

    # create empty group plot
    @pgf gp = GroupPlot(
        {
            group_style = {
                group_size = "2 by 2",
                xlabels_at = "edge bottom",
                ylabels_at = "edge left",
                horizontal_sep = "1cm",
                vertical_sep = "1.25cm"
            },
            xlabel = raw"$t$", ylabel = raw"$A(t)$",
            x_unit = raw"\hour", y_unit = raw"\mol\per\litre",
            xmin = prob_dde_qs.tspan[1], xmax = prob_dde_qs.tspan[2],
            width = "5.5cm", height = "3.5cm",
            grid = "major",
            label_style = { font = raw"\footnotesize" },
            tick_label_style = { font = raw"\footnotesize" },
            title_style = { font = raw"\footnotesize\bfseries" },
        }
    )

    # for different ODE solvers
    for alg in (BS3(), Tsit5(), Rosenbrock23(), Rodas4())
        # compute solution of third component
        sol = solve(prob_dde_qs, MethodOfSteps(alg); save_idxs = 3)

        # plot AHL dynamics
        @pgf push!(
            gp,
            {
                title = string(nameof(typeof(alg)))
            },
            PlotInc(
                {
                    no_markers
                },
                Table(; x = ts, y = sol(ts))),
            PlotInc(
                {
                    only_marks, mark = "o"
                },
                Table(; x = sol.t, y = sol.u)))
    end

    gp
end

pgfsave(joinpath(@__DIR__, "..", "figures", "generated", "qs_comparison.tex"),
        plot_qs_comparison();
        include_preamble = false)


text:
using ThreeBodyDecay
using Test

rands = rand(3)
σs = Invariants(rands...)

@testset "Invariants structure" begin
    #
    @test σs.σ1 == σs[1] == rands[1]
    @test σs.σ2 == σs[2] == rands[2]
    @test σs.σ3 == σs[3] == rands[3]
    #
    @test_throws ErrorException σs[5]
end

@testset "operations and interate" begin
    @test sum(σs) == sum(rands)
    @test length(σs) == 3
end

ms = ThreeBodyMasses(1.1, 3.3,5.5; m0=20.0)
σs = Invariants(ms; σ3=4.5^2, σ2=10.1^2)

@testset "Creation from two given" begin
    @test sum(σs .≈ Invariants(ms; σ1=σs.σ1, σ2=σs.σ2))==3
    @test sum(σs .≈ Invariants(ms; σ2=σs.σ2, σ3=σs.σ3))==3
    @test sum(σs .≈ Invariants(ms; σ3=σs.σ3, σ1=σs.σ1))==3
    #
    @test nt(σs).σ1 == σs.σ1
    @test_throws ErrorException Invariants(ms; σ1=1.0, σ2=1.0, σ3=1.0)
    #
    @test Kibble(σs,ms^2) < 0
end

@testset "Random points" begin
    σs = randomPoint(ms)
    @test Kibble(σs,ms^2) < 0
end

tbs = ThreeBodySystem(ms)

@testset "Three-body system" begin
    @test tbs == ThreeBodySystem(ms=ms,
        two_js = ThreeBodySpins(0,0,0;two_h0=0))
    @test tbs == ThreeBodySystem(ms,
        ThreeBodySpins(0,0,0;two_h0=0))
    @test tbs == ThreeBodySystem(ms.m1, ms.m2, ms.m3; m0=ms.m0,
        two_js = ThreeBodySpins(0,0,0;two_h0=0))
end

tbsS = ThreeBodySystem(ms,ThreeBodySpins(1,3,5;two_h0=1))
rp = randomPoint(tbsS)

@testset "Random points" begin
    @test Kibble(rp.σs, tbsS.ms^2) < 0
    @test Tuple(rp.two_λs) ∈ collect(itr(tbsS.two_js))
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
module ConnectionRequest

import ..Layer, ..request
using ..URIs, ..Sockets
using ..Messages
using ..IOExtras
using ..ConnectionPool
using MbedTLS: SSLContext
import ..@debug, ..DEBUG_LEVEL

"""
    request(ConnectionPoolLayer, ::URI, ::Request, body) -> HTTP.Response

Retrieve an `IO` connection from the [`ConnectionPool`](@ref).

Close the connection if the request throws an exception.
Otherwise leave it open so that it can be reused.

`IO` related exceptions from `Base` are wrapped in `HTTP.IOError`.
See [`isioerror`](@ref).
"""
abstract type ConnectionPoolLayer{Next <: Layer} <: Layer end
export ConnectionPoolLayer

function request(::Type{ConnectionPoolLayer{Next}}, url::URI, req, body;
                 proxy=nothing,
                 socket_type::Type=TCPSocket,
                 reuse_limit::Int=ConnectionPool.nolimit, kw...) where Next

    if proxy !== nothing
        target_url = url
        url = URI(proxy)
        if target_url.scheme == "http"
            req.target = string(target_url)
        end
    end

    IOType = ConnectionPool.Transaction{sockettype(url, socket_type)}
    local io
    try
        io = getconnection(IOType, url.host, url.port;
                           reuse_limit=reuse_limit, kw...)
    catch e
        rethrow(isioerror(e) ? IOError(e, "during request($url)") : e)
    end

    if io.sequence >= reuse_limit
        defaultheader!(req, "Connection" => "close")
    end

    try
        if proxy !== nothing && target_url.scheme == "https"
            return tunnel_request(Next, io, target_url, req, body; kw...)
        end

        return request(Next, io, req, body; kw...)

    catch e
        @debug 1 "❗  ConnectionLayer $e. Closing: $io"
        close(io)
        rethrow(isioerror(e) ? IOError(e, "during request($url)") : e)
    finally
        if (io.sequence >= reuse_limit
        || (proxy !== nothing && target_url.scheme == "https"))
            close(io)
        end
    end

end

sockettype(url::URI, default) = url.scheme in ("wss", "https") ? SSLContext :
                                                                 default

function tunnel_request(Next, io, target_url, req, body; kw...)
    r = connect_tunnel(io, target_url, req)
    if r.status != 200
        return r
    end
    io = ConnectionPool.sslupgrade(io, target_url.host; kw...)
    req.headers = filter(x->x.first != "Proxy-Authorization", req.headers)
    return request(Next, io, req, body; kw...)
end

function connect_tunnel(io, target_url, req)
    target = "$(URIs.hoststring(target_url.host)):$(target_url.port)"
    @debug 1 "📡  CONNECT HTTPS tunnel to $target"
    headers = Dict(filter(x->x.first == "Proxy-Authorization", req.headers))
    request = Request("CONNECT", target, headers)
    writeheaders(io, request)
    startread(io)
    readheaders(io, request.response)
    return request.response
end

end # module ConnectionRequest


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""SecurityContext holds security configuration that will be applied to a container. Some fields are present in both SecurityContext and PodSecurityContext.  When both are set, the values in SecurityContext take precedence.

    IoK8sApiCoreV1SecurityContext(;
        allowPrivilegeEscalation=nothing,
        capabilities=nothing,
        privileged=nothing,
        procMount=nothing,
        readOnlyRootFilesystem=nothing,
        runAsGroup=nothing,
        runAsNonRoot=nothing,
        runAsUser=nothing,
        seLinuxOptions=nothing,
        windowsOptions=nothing,
    )

    - allowPrivilegeEscalation::Bool : AllowPrivilegeEscalation controls whether a process can gain more privileges than its parent process. This bool directly controls if the no_new_privs flag will be set on the container process. AllowPrivilegeEscalation is true always when the container is: 1) run as Privileged 2) has CAP_SYS_ADMIN
    - capabilities::IoK8sApiCoreV1Capabilities : The capabilities to add/drop when running containers. Defaults to the default set of capabilities granted by the container runtime.
    - privileged::Bool : Run container in privileged mode. Processes in privileged containers are essentially equivalent to root on the host. Defaults to false.
    - procMount::String : procMount denotes the type of proc mount to use for the containers. The default is DefaultProcMount which uses the container runtime defaults for readonly paths and masked paths. This requires the ProcMountType feature flag to be enabled.
    - readOnlyRootFilesystem::Bool : Whether this container has a read-only root filesystem. Default is false.
    - runAsGroup::Int64 : The GID to run the entrypoint of the container process. Uses runtime default if unset. May also be set in PodSecurityContext.  If set in both SecurityContext and PodSecurityContext, the value specified in SecurityContext takes precedence.
    - runAsNonRoot::Bool : Indicates that the container must run as a non-root user. If true, the Kubelet will validate the image at runtime to ensure that it does not run as UID 0 (root) and fail to start the container if it does. If unset or false, no such validation will be performed. May also be set in PodSecurityContext.  If set in both SecurityContext and PodSecurityContext, the value specified in SecurityContext takes precedence.
    - runAsUser::Int64 : The UID to run the entrypoint of the container process. Defaults to user specified in image metadata if unspecified. May also be set in PodSecurityContext.  If set in both SecurityContext and PodSecurityContext, the value specified in SecurityContext takes precedence.
    - seLinuxOptions::IoK8sApiCoreV1SELinuxOptions : The SELinux context to be applied to the container. If unspecified, the container runtime will allocate a random SELinux context for each container.  May also be set in PodSecurityContext.  If set in both SecurityContext and PodSecurityContext, the value specified in SecurityContext takes precedence.
    - windowsOptions::IoK8sApiCoreV1WindowsSecurityContextOptions : The Windows specific settings applied to all containers. If unspecified, the options from the PodSecurityContext will be used. If set in both SecurityContext and PodSecurityContext, the value specified in SecurityContext takes precedence.
"""
mutable struct IoK8sApiCoreV1SecurityContext <: SwaggerModel
    allowPrivilegeEscalation::Any # spec type: Union{ Nothing, Bool } # spec name: allowPrivilegeEscalation
    capabilities::Any # spec type: Union{ Nothing, IoK8sApiCoreV1Capabilities } # spec name: capabilities
    privileged::Any # spec type: Union{ Nothing, Bool } # spec name: privileged
    procMount::Any # spec type: Union{ Nothing, String } # spec name: procMount
    readOnlyRootFilesystem::Any # spec type: Union{ Nothing, Bool } # spec name: readOnlyRootFilesystem
    runAsGroup::Any # spec type: Union{ Nothing, Int64 } # spec name: runAsGroup
    runAsNonRoot::Any # spec type: Union{ Nothing, Bool } # spec name: runAsNonRoot
    runAsUser::Any # spec type: Union{ Nothing, Int64 } # spec name: runAsUser
    seLinuxOptions::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SELinuxOptions } # spec name: seLinuxOptions
    windowsOptions::Any # spec type: Union{ Nothing, IoK8sApiCoreV1WindowsSecurityContextOptions } # spec name: windowsOptions

    function IoK8sApiCoreV1SecurityContext(;allowPrivilegeEscalation=nothing, capabilities=nothing, privileged=nothing, procMount=nothing, readOnlyRootFilesystem=nothing, runAsGroup=nothing, runAsNonRoot=nothing, runAsUser=nothing, seLinuxOptions=nothing, windowsOptions=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1SecurityContext, Symbol("allowPrivilegeEscalation"), allowPrivilegeEscalation)
        setfield!(o, Symbol("allowPrivilegeEscalation"), allowPrivilegeEscalation)
        validate_property(IoK8sApiCoreV1SecurityContext, Symbol("capabilities"), capabilities)
        setfield!(o, Symbol("capabilities"), capabilities)
        validate_property(IoK8sApiCoreV1SecurityContext, Symbol("privileged"), privileged)
        setfield!(o, Symbol("privileged"), privileged)
        validate_property(IoK8sApiCoreV1SecurityContext, Symbol("procMount"), procMount)
        setfield!(o, Symbol("procMount"), procMount)
        validate_property(IoK8sApiCoreV1SecurityContext, Symbol("readOnlyRootFilesystem"), readOnlyRootFilesystem)
        setfield!(o, Symbol("readOnlyRootFilesystem"), readOnlyRootFilesystem)
        validate_property(IoK8sApiCoreV1SecurityContext, Symbol("runAsGroup"), runAsGroup)
        setfield!(o, Symbol("runAsGroup"), runAsGroup)
        validate_property(IoK8sApiCoreV1SecurityContext, Symbol("runAsNonRoot"), runAsNonRoot)
        setfield!(o, Symbol("runAsNonRoot"), runAsNonRoot)
        validate_property(IoK8sApiCoreV1SecurityContext, Symbol("runAsUser"), runAsUser)
        setfield!(o, Symbol("runAsUser"), runAsUser)
        validate_property(IoK8sApiCoreV1SecurityContext, Symbol("seLinuxOptions"), seLinuxOptions)
        setfield!(o, Symbol("seLinuxOptions"), seLinuxOptions)
        validate_property(IoK8sApiCoreV1SecurityContext, Symbol("windowsOptions"), windowsOptions)
        setfield!(o, Symbol("windowsOptions"), windowsOptions)
        o
    end
end # type IoK8sApiCoreV1SecurityContext

const _property_map_IoK8sApiCoreV1SecurityContext = Dict{Symbol,Symbol}(Symbol("allowPrivilegeEscalation")=>Symbol("allowPrivilegeEscalation"), Symbol("capabilities")=>Symbol("capabilities"), Symbol("privileged")=>Symbol("privileged"), Symbol("procMount")=>Symbol("procMount"), Symbol("readOnlyRootFilesystem")=>Symbol("readOnlyRootFilesystem"), Symbol("runAsGroup")=>Symbol("runAsGroup"), Symbol("runAsNonRoot")=>Symbol("runAsNonRoot"), Symbol("runAsUser")=>Symbol("runAsUser"), Symbol("seLinuxOptions")=>Symbol("seLinuxOptions"), Symbol("windowsOptions")=>Symbol("windowsOptions"))
const _property_types_IoK8sApiCoreV1SecurityContext = Dict{Symbol,String}(Symbol("allowPrivilegeEscalation")=>"Bool", Symbol("capabilities")=>"IoK8sApiCoreV1Capabilities", Symbol("privileged")=>"Bool", Symbol("procMount")=>"String", Symbol("readOnlyRootFilesystem")=>"Bool", Symbol("runAsGroup")=>"Int64", Symbol("runAsNonRoot")=>"Bool", Symbol("runAsUser")=>"Int64", Symbol("seLinuxOptions")=>"IoK8sApiCoreV1SELinuxOptions", Symbol("windowsOptions")=>"IoK8sApiCoreV1WindowsSecurityContextOptions")
Base.propertynames(::Type{ IoK8sApiCoreV1SecurityContext }) = collect(keys(_property_map_IoK8sApiCoreV1SecurityContext))
Swagger.property_type(::Type{ IoK8sApiCoreV1SecurityContext }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1SecurityContext[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1SecurityContext }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1SecurityContext[property_name]

function check_required(o::IoK8sApiCoreV1SecurityContext)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1SecurityContext }, name::Symbol, val)
end


text:
#=

Implements the empirical cumulative distribution function given an array
of observations.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-08-25

=#

@deprecate ECDF(x) StatsBase.ecdf(x)


text:
"""

  Construct the sylvester matrix of degree k from polynomials f and g = f'.
  The matrix has 2k+1 columns. Each column contains a shifted version of f or g.
  Order of columns: k+1 * g, k * f
"""
function sylves1(f::AbstractVector{T}, g::AbstractVector{T}, k::Int) where {T<:Number}
#
#  The k-th sylvester matrix of f(x), f'(x) = g(x)
#
    l = length(f)
    m = l + k - 1
    n = 2 * k + 1
    A = zeros(T, m, n)
	for j = 1:k+1
        A[j:j+l-2,j] = g
    end
    for j = 1:k
        A[j:j+l-1,j+k+1] = f
    end
    A
end

"""
Extract result polynomials u and v if x is a solution of sylves1(f,g) * x == 0.
  u * g + v * f == 0, degree of v is k.
"""
function extract_sylves1(x::AbstractVector{T}) where {T<:Number}
  k = length(x)
  k > 2 && k % 2 == 1 || error("vector too short or not odd length")
  k = (k + 1) ÷ 2
  x[1:k] / x[1], x[k+1:end] / x[k+1]  # u, v
end

"""

  Construct the sylvester matrix of degree k from polynomials f and g = f'.
  The matrix has 2k+1 columns. Each column contains a shifted version of f or g.
  Order of columns: g, k * [f g]
"""
function sylves2(f::AbstractVector{T}, g::AbstractVector{T}, k::Int) where {T<:Number}
#
#  The k-th sylvester matrix of f(x), f'(x) = g(x)
#
    l = length(f)
    m = l + k - 1
    n = 2 * k + 1
    A = zeros(T, m, n)
	for j = 1:k+1
        A[j:j+l-2,j*2-1] = g
    end
    for j = 1:k
        A[j:j+l-1,j*2] = f
    end
    A
end

"""
Extract result polynomials u and v if x is a solution of sylves2(f,g) * x == 0.
  u * g + v * f == 0, degree of v is k.
"""
function extract_sylves2(x::AbstractVector{T}) where {T<:Number}
  k = length(x)
  k > 2 && k % 2 == 1 || error("vector too short or not odd length")
  x[1:2:k] / x[1], x[2:2:k] / x[2]  # u, v
end


text:
module Blockfield

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Strawberry Blockfield" => Ahorn.EntityPlacement(
        Maple.StrawberryBlockField,
        "rectangle"
    ),
)

Ahorn.minimumSize(entity::Maple.StrawberryBlockField) = 8, 8
Ahorn.resizable(entity::Maple.StrawberryBlockField) = true, true

Ahorn.selection(entity::Maple.StrawberryBlockField) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.StrawberryBlockField, room::Maple.Room)
    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (0.4, 0.4, 1.0, 0.4), (0.4, 0.4, 1.0, 1.0))
end

end


text:
@testset "IndexedReader" begin
    buf = UInt8[]
    writer = BSONWriter(buf)
    writer["x"] = Int64(1)
    writer["y"] = w -> begin
        w["x"] = Int64(2)
    end
    close(writer)
    reader = BSONReader(buf)
    ireader = IndexedBSONReader(BSONIndex(128), reader)
    @test ireader["x"][Int64] == 1
    @test ireader["y"]["x"][Int64] == 2
end


text:
import Distributions: Normal, pdf, cdf

###############################
# Type Definition/Constructors
###############################
mutable struct Grid
    points::Vector{Float64}
    weights::Vector{Float64}
    scale::Float64
    Grid(points, weights, scale) = sum(weights) ≈ scale ? new(points, weights, scale) : error("scaled weights do not sum up properly")
end

# Constructor utilizing a custom weight calculation function
function Grid(quadrature::Function,
              lower_bound::T, upper_bound::T,
              n_points::Int; scale::Float64 = 1.) where {T<:Real}
    grid, weights = quadrature(lower_bound, upper_bound, n_points)
    Grid(reshape(grid, size(grid)[1]), weights, scale)
end

####################
# Quadrature rules
####################
# Enforce that quadrature rules with additional keyword arguments get mapped to
# A pre-populated version of that rule that only accepts 3 arguments:
# Lower bound, upper bound, and number of points
function uniform_quadrature(lower_bound::T, upper_bound::T, n_points::Int;
                            scale::T = 1) where {T<:Real}
    grid = collect(range(lower_bound, stop = upper_bound, length = n_points))
    weights = fill(scale/n_points, n_points)
    return grid, weights
end

function uniform_quadrature(scale::T = 1) where {T<:Real}
    return (lb, ub, n) -> uniform_quadrature(lb, ub, n, scale = scale)
end

# NOTE: Should modify chebpts function to return cleaner types...
# Kind indicates the "kind" of the chebyshev polynomial grid
# The terminology is generally, "Chebyshev polynomials of the Nth kind"
function curtis_clenshaw_quadrature(lower_bound::T, upper_bound::T,
                                    n_points::Int; kind::Int = 2) where {T<:Real}
    grid, weights = chebpts(n_points, lower_bound, upper_bound, kind)
    return grid', dropdims(weights', dims=2)
end

function curtis_clenshaw_quadrature(kind::Int = 2)
    return (lb, ub, n) -> curtis_clenshaw_quadrature(lb, ub, n, kind = kind)
end

# Tauchen86 developed a method for "finding a discrete-valued Markov chain
# whose sample paths approximate well those of a vector autoregression"
# Hence the method below is the Tauchen86 method applied to a general AR(1).
function tauchen86(μ::AbstractFloat,ρ::AbstractFloat,σ::AbstractFloat,n::Int64,λ::AbstractFloat)
    #output is xgrid, xprob
    # x_t+1 = μ + ρ x_t + σ e_{t+1}, e_t+1 ∼ N(0,1)
    xhi = μ/(1-ρ) + λ*sqrt(σ^2/(1-ρ)^2)
    xlo = μ/(1-ρ) - λ*sqrt(σ^2/(1-ρ)^2)
    xgrid = zeros(n);
    xscale =(xhi-xlo)/(n-1)

    for i=1:n
        xgrid[i] = xlo + xscale*(i-1)
    end
    m=zeros(n-1);
    for i=1:n-1
        m[i] = (xgrid[i]+xgrid[i+1])/2
    end
    xprob = zeros(n,n) # xprob[i,j] = Pr(x_t+1=xgrid[j]|x_t=xgrid[i])
    for i=1:n # this is the state today
        normie = Normal(μ+ρ*xgrid[i],σ)
        normpdf(x) = pdf.(normie,x)
        normcdf(x) = cdf.(normie,x)
        for j=2:n-1
            xprob[i,j] = normcdf(m[j]) - normcdf(m[j-1])
        end
        xprob[i,1] = normcdf(m[1])
        xprob[i,n] = 1 - normcdf(m[n-1])
    end
    xprob = xprob./sum(xprob, dims = 2) # make sure the rows sum to 1
    return ( xgrid,xprob, xscale )
end

# However, Tauchen86 can also be applied to "degenerate" AR(1)s where ρ = 0.
# i.e. i.i.d noise with mean μ, std. σ, but the way the returned
# Markov transition matrix, xprob, is used in the code, assumes it to be
# 1-dimensional, since the second-dimension of the Markov transition matrix is trivial
# (since the process is i.i.d).
function tauchen86(μ::AbstractFloat, σ::AbstractFloat, n::Int64, λ::AbstractFloat)
    xgrid, xprob, xscale = tauchen86(μ, 0., σ, n, λ)
    return xgrid, vec(xprob[:, 1]), xscale
end

####################
# Grid-based utils
####################
function get_grid(m::AbstractDSGEModel, grid_name::Symbol)
    return m.grids[grid_name]
end

function quadrature_sum(x::Vector{T}, grid::Grid) where {T<:Real}
    return sum(grid.weights .* x .* grid.points)
end

function quadrature_sum(grid::Grid, x::Vector{T}) where {T<:Real}
    return sum(x, grid)
end

# # Defining arithmetic on/standard function evaluation of grids
# for op in (:(Base.:+),
           # :(Base.:-),
           # :(Base.:*),
           # :(Base.:/))
           # # Exponentiation also omitted because of broadcasting issues

    # @eval ($op)(g::Grid, x::Integer)           = ($op)(g.points, x)
    # @eval ($op)(g::Grid, x::Number)            = ($op)(g.points, x)
    # if op == :(Base.:/)
        # op = :(Base.:./) # Will not work because the syntax for broadcasting has been
                         # # updated to broadcast(/, ...)?
        # @eval ($op)(x::Integer, g::Grid)       = ($op)(x, g.points)
        # @eval ($op)(x::Number, g::Grid)        = ($op)(x, g.points)
    # end
# end

# for f in (:(Base.:-),
          # :(Base.:<),
          # :(Base.:>),
          # :(Base.:<=),
          # :(Base.:>=))

    # @eval ($f)(g::Grid) = ($f)(g.points)

    # # Similar issues with broadcasting...
    # if f != :(Base.:-)
        # @eval ($f)(g::Grid, x::Number) = ($f)(g.points, x)
    # end
# end


text:
# Some segments appear good when checking the quality,
const IDENTIFIERS_BAD = Set([
    # These three come from the same lab, and I have a strong suspicion they
    # are chimeric assemblies - there is recombination at pos ~230
    "MK168420",
    "MK168525",
    "MK168545",
    # This one has a 75 bp deletion that does not seem plausible
    "KY926084",
])

"""
    isok_all(::SegmentData)

Check if the `SegmentData` passes all criteria. Checks:
* It's not a known bad segment
* Segment has less than 5 ambiguous bases
* `isok_segment_label`
* `isok_minimum_proteins`
* `isok_seq_length` (inf A only, skipped for B)
* `isok_orf_length` (inf A only, skipped for B)
* `isok_translatable`
"""
function isok_all(data::SegmentData)::Bool
    (!in(data.id, IDENTIFIERS_BAD)) &&
    count(BioSequences.isambiguous, data.seq) < 5 &&
    isok_segment_label(data) &&
    isok_minimum_proteins(data) &&
    isok_seq_length(data) &&
    isok_orf_length(data) &&
    isok_translatable(data)
end

"""
    isok_segment_label(::SegmentData)

Check whether a `SegmentData`'s DNA and protein variant is consistent, i.e. no PB2
protein in a PB1 segment. This may happen when segments gets mislabeled, e.g.
KT853292.1 is labeled "segment 2", but is PB2 (segment 1).
"""
function isok_segment_label(data::SegmentData)::Bool
    all(p -> source(p.var) == data.segment, data.proteins)
end

"""
    isok_minimum_proteins(::SegmentData)

Checks whether a segment has the minimum number of ORFS to be considered "complete",
i.e. whether it does not miss any non-auxilliary proteins.
"""
function isok_minimum_proteins(data::SegmentData)::Bool
    segment = data.segment
    serotype = data.serotype
    minimum = if segment == Segments.PB2
        (Proteins.PB2,)
    elseif segment == Segments.PB1
        (Proteins.PB1,)
    elseif segment == Segments.PA
        (Proteins.PA,)
    elseif segment == Segments.HA
        (Proteins.HA,)
    elseif segment == Segments.NP
        (Proteins.NP,)
    elseif segment == Segments.NA
        # NB is auxilliary
        (Proteins.NA,)
    elseif segment == Segments.MP
        is_error(serotype) ? (Proteins.M1, Proteins.BM2) : (Proteins.M1, Proteins.M2)
    elseif segment == Segments.NS
        (Proteins.NS1, Proteins.NEP)
    else
        @assert false "Unreachable!" # This helps inference
    end
    return !any(minimum) do protein
        isnothing(findfirst(i -> i.var == protein, data.proteins))
    end
end

# This has been empirially determined by looking at the distributions
# in order to find outliers
const ACCEPTABLE_SEGMENT_LENGTHS = Dict(
    Segments.NP => 1490:1570,
    Segments.HA => 1680:1780,
    Segments.MP => 980:1030,
    Segments.PB1 => 2270:2345,
    Segments.PA => 2150:2240,
    Segments.NA => 1350:1470,  # very wide distribution??
    Segments.NS => 820:895,
    Segments.PB2 => 2280:2345
)

"""
    isok_seq_length(::SegmentData)

Check whether a segment's sequence length falls within a reasonable length interval.
Always returns `true` for Influenza B segments.
"""
function isok_seq_length(data::SegmentData)::Bool
    # The sizes for influenza B are a little different.
    is_error(data.serotype) && return true
    length(data.seq) in ACCEPTABLE_SEGMENT_LENGTHS[data.segment]
end

# This is pretty lax, so can be tightened later. It's just to remove
# obvious annotation errors.
const ACCEPTABLE_ORF_LENGTHS = Dict(
    Proteins.PB2 => 2250:2400,
    Proteins.PB1 => 2200:2400,
    Proteins.N40 => nothing, # nothings mean too little data to know range
    Proteins.PB1F2 => 0:350, # often truncated in real viruses
    Proteins.PA => 2100:2250,
    Proteins.PAX => 600:800,
    Proteins.HA => 1600:1800,
    Proteins.NP => 1450:1550,
    Proteins.NA => 1300:1450,
    Proteins.NB => nothing,
    Proteins.M1 => 720:780,
    Proteins.M2 => 250:325,
    Proteins.BM2 => nothing,
    Proteins.M42 => nothing,
    Proteins.NS1 => 600:750,
    Proteins.NEP => 325:400,
    Proteins.NS3 => nothing
)

# Yeah, some of these ORFs are _clearly_ not good enough.
"""
    isok_orf_length(::SegmentData)

Check whether a segment's ORF length falls within a reasonable length interval.
If the segment does not have a SeroType, returns `true.`.
"""
function isok_orf_length(data::SegmentData)::Bool
    is_error(data.serotype) && return true

    for protein in data.proteins
        len = sum(length, protein.orfs, init=0)
        range = ACCEPTABLE_ORF_LENGTHS[protein.var]
        if range !== nothing && !(len in range)
            return false
        end
    end
    return true
end

"""
    isok_translatable(::SegmentData)

Check whether a segment's ORFs are translatable:
* If the segment does not have a sequence or any ORFS, return `false`.
* If any ORF exceed past the segment sequence, return `false`
* If the joined ORFs' length is not divisible by 3, return `false`
* If any stops in sequence, return `false`.
* If sequence is not followed by a stop codon, return false
"""
function isok_translatable(data::SegmentData)::Bool
    nt_sequence = LongDNASeq()
    for protein in data.proteins
        if any(last(orf) > length(data.seq) for orf in protein.orfs)
            return false
        end
        join!(nt_sequence, (@view(data.seq[orf]) for orf in protein.orfs))

        # Must have a nonzerolength divisible by 3
        (length(nt_sequence) > 2 || iszero(length(nt_sequence) % 3)) || return false
        aa_sequence = BioSequences.translate(nt_sequence)

        # Must not contain stop
        findfirst(AA_Term, aa_sequence) === nothing || return false

        # There must be a stop just after the seq
        lastpos = protein.orfs[end][end]
        stoprange = lastpos+1:lastpos+3
        checkbounds(Bool, eachindex(data.seq), stoprange) || return false
        stopseq = data.seq[stoprange]
        if any(BioSequences.isambiguous, stopseq) || !Influenza.is_stop(BioSequences.DNACodon(stopseq))
            return false
        end
    end
    return true
end

function join!(seq::BioSequence, it)
    len = sum(length, it, init=0)
    length(seq) != len && resize!(seq, len)
    offset = 0
    for i in it
        seq[offset+1:offset+length(i)] = i
        offset += length(i)
    end
    @assert offset == len
    return seq
end


text:
# 2D FWI gradient test with 4 sources
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using PyCall, PyPlot, JUDI.TimeModeling, Test, LinearAlgebra

## Set up model structure
n = (120,100)	# (x,y,z) or (x,z)
d = (10.,10.)
o = (0.,0.)

# Velocity [km/s]
v = ones(Float32,n) .* 2.0f0
v[:,Int(round(end/2)):end] .= 3.0f0
v0 = smooth10(v,n)
rho = ones(Float32, n)
rho[:, Int(round(end/2)):end] .= 1.5f0

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2
m0 = (1f0 ./ v0).^2
dm = m0 - m

# Setup info and model structure
nsrc = 1	# number of sources
ntComp = 250
info = Info(prod(n), nsrc, ntComp)	# number of gridpoints, number of experiments, number of computational time steps
model = Model(n,d,o,m;rho=rho)
model0 = Model(n,d,o,m0;rho=rho)

## Set up receiver geometry
nxrec = 81
xrec = range(200f0,stop=1000f0,length=nxrec)
yrec = 0.
zrec = range(100f0,stop=100f0,length=nxrec)

# receiver sampling and recording time
timeR = 1000f0	# receiver recording time [ms]
dtR = 4f0	# receiver sampling interval

# Set up receiver structure
recGeometry = Geometry(xrec,yrec,zrec;dt=dtR,t=timeR,nsrc=nsrc)

## Set up source geometry (cell array with source locations for each shot)
xsrc = convertToCell([600f0])
ysrc = convertToCell([0f0])
zsrc = convertToCell([50f0])

# source sampling and number of time steps
timeS = 1000f0
dtS = 4f0	# receiver sampling interval

# Set up source structure
srcGeometry = Geometry(xsrc,ysrc,zsrc;dt=dtS,t=timeS)

# setup wavelet
f0 = 0.01f0
wavelet = ricker_wavelet(timeS,dtS,f0)

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

# Gradient test
h = 1f0
iter = 8
error1 = zeros(iter)
error2 = zeros(iter)
h_all = zeros(iter)
srcnum = 1:nsrc
modelH = deepcopy(model0)

# Observed data
opt = Options(sum_padding=true)
F = judiModeling(info,model,srcGeometry,recGeometry;options=opt)
q = judiVector(srcGeometry,wavelet)
d = F*q

# FWI gradient and function value for m0
Jm0, grad = fwi_objective(model0,q,d;options=opt)

for j=1:iter
	# FWI gradient and function falue for m0 + h*dm
	modelH.m = model0.m + h*dm
	Jm, gradm = fwi_objective(modelH,q,d;options=opt)

	dJ = dot(grad,vec(dm))

	# Check convergence
	error1[j] = abs(Jm - Jm0)
	error2[j] = abs(Jm - (Jm0 + h*dJ))

	println(h, " ", error1[j], " ", error2[j])
	h_all[j] = h
	global h = h/2f0
end

# Check error decay
rate_0th_order = 2^(iter - 1)   # error decays w/ factor 2
rate_1st_order = 4^(iter - 1)   # error decays w/ factor 4

@test error1[end] <= error1[1] / rate_0th_order
@test error2[end] <= error2[1] / rate_1st_order

# Plot errors
if isinteractive()
    loglog(h_all, error1); loglog(h_all, 1e2*h_all)
    loglog(h_all, error2); loglog(h_all, 1e2*h_all.^2)
    legend([L"$\Phi(m) - \Phi(m0)$", "1st order", L"$\Phi(m) - \Phi(m0) - \nabla \Phi \delta m$", "2nd order"], loc="lower right")
    xlabel("h")
    ylabel(L"Error $||\cdot||^\infty$")
    title("FWI gradient test")
end


text:
##################################################
#   Override default handling (default = inline) #
##################################################

###############################################
#         Overload environment functions      #
###############################################

function latexraw(x::SymEngine.Basic; kwargs...)
    str = string(x)
    ex = Meta.parse(str)
    latexraw(ex; kwargs...)
end

###############################################
#         Overload utilities                  #
###############################################

add_brackets(syms::SymEngine.Basic, vars) = add_brackets(Meta.parse("$syms"), vars)


text:
import NeuroCore: EncodingDirection

@testset "Encoding directions" begin
    @test EncodingDirection(1) == EncodingDirection("ipos") == EncodingDirection("i")
    @test EncodingDirection(-1) == EncodingDirection("ineg") == EncodingDirection("i-")
    @test EncodingDirection(2) == EncodingDirection("jpos") == EncodingDirection("j")
    @test EncodingDirection(-2) == EncodingDirection("jneg") == EncodingDirection("j-")
    @test EncodingDirection(3) == EncodingDirection("kpos") == EncodingDirection("k")
    @test EncodingDirection(-3) == EncodingDirection("kneg") == EncodingDirection("k-")
    @test String(EncodingDirection(-3)) == "kneg"

    @test_throws ErrorException EncodingDirection(:xyz)
end


text:
function discount!(rs::AbstractVector{T}, γ::T) where T <: AbstractFloat
    @inbounds for i in length(rs)-1:-1:1
        rs[i] += γ * rs[i+1]
    end
end

function nstepr!(rs::AbstractVector{T}, brs::AbstractVector{T}, γ::T, n::Int64) where T <: AbstractFloat
    @inbounds for t in 1:length(rs)
        runningadd = zero(T)
        l = min(n, length(rs)-t)
        for i in 0:l
            runningadd += γ^(i) * (i == n ? brs[t+i] : rs[t+i])
        end
        rs[t] = runningadd
    end
end

function GAEλ!(vs::AbstractVector{T}, rs::AbstractVector{T}, γ::T, λ::T) where T <: AbstractFloat
    vs[1:end-1] .= rs[1:end-1] .+ γ .* vs[2:end] .- vs[1:end-1]
    vs[end] = rs[end] - vs[end]
    discount!(vs, γ*λ)
end

function norm!(vs::AbstractVector{T}) where T <: AbstractFloat
    θmean = mean(vs)
    θstd = std(vs, mean=θmean)
    vs .= (vs .- θmean) ./ (θstd + T(1f-12))
end


text:
export readjats

using EzXML

struct UnsupportedException <: Exception
    message::String
end

"""
Read JATS file and returns `Tree` instance.
"""
function readjats(path::String)
    try
        xml_article = root(readxml(path))
        @assert nodename(xml_article) == "article"
        if countelements(xml_article) < 3
            warn("#xml elements < 3")
            return
        end

        article = Tree("article")
        xml_front = findfirst(xml_article, "front")
        push!(article, parse_front(xml_front))
        # push!(article[end], Tree("pdf-link",Tree("http://www.aclweb.org/anthology/P12-1046")))
        # push!(article[end], Tree("xml-link",Tree("http://example.xml")))

        body = findfirst(xml_article, "body")
        push!(article, parse_body(body))

        back = find(xml_article, "back")
        if !isempty(back)
            push!(article, parse_back(back[1]))
        end

        push!(article, Tree("floats-group"))
        append!(article[end], findfloats(article))
        floats = find(xml_article, "floats-group")
        if !isempty(floats)
            append!(article[end], parse_body(floats[1]).children)
        end
        isempty(article[end]) && deleteat!(article,length(article)) # no floats

        postprocess!(article)
        return article
    catch e
        if isa(e, UnsupportedException)
            println(e.message)
            return
        else
            rethrow(e)
        end
    end
end

"""
Convert EzXML.Node into `Tree`.

* dict: dictionary of keep nodes
"""
function Base.convert(::Type{Tree}, enode::EzXML.Node, dict=nothing)
    elements = filter(nodes(enode)) do n
        iselement(n) && return true
        istext(n) && !ismatch(r"^\s*$",nodecontent(n))
    end
    elements = collect(elements)
    tempnodes = map(elements) do e
        istext(e) ? Tree(nodecontent(e)) : convert(Tree,e,dict)
    end
    deletable = begin
        if any(istext, elements)
            false
        elseif dict == nothing || any(e -> haskey(dict,e), elements)
            true
        elseif any(n -> any(!isempty,n.children), tempnodes)
            true
        else
            false
        end
    end

    children = Tree[]
    for i = 1:length(elements)
        e = elements[i]
        t = tempnodes[i]
        if istext(e) || dict == nothing || haskey(dict,e)
            if isempty(children)
                push!(children, t)
            elseif isempty(t) && isempty(children[end])
                children[end].name *= t.name
            else
                push!(children, t)
            end
        elseif any(!isempty, t.children) || !deletable
            for c in t.children
                if isempty(children)
                    push!(children, c)
                elseif isempty(c) && isempty(children[end])
                    children[end].name *= c.name
                else
                    push!(children, c)
                end
            end
        end
    end
    Tree(nodename(enode), children)
end

function parse_front(front::EzXML.Node)
    xpath = """
        journal-meta/journal-title-group/journal-title
        | article-meta/title-group/article-title
        | article-meta/abstract
        """
    tree = Tree(nodename(front), map(parse_body,find(front,xpath)))

    contrib = "article-meta/contrib-group/contrib[@contrib-type=\"author\"]"
    xpath = """
        $contrib
        | $contrib/name/prefix
        | $contrib/name/given-names
        | $contrib/name/surname
        | $contrib/name/suffix
        | $contrib/collab
        """
    dict = Dict(n => n for n in find(front,xpath))
    authors = convert(Tree, front, dict).children
    foreach(a -> a.name = "author", authors)
    append!(tree, authors)
    tree
end

function parse_body(body::EzXML.Node)
    xpath = """
        //boxed-text
        | //boxed-text/label
        | //boxed-text/caption
        | //boxed-text/caption/title
        | //code
        | //def-list
        | //def-list/label
        | //def-list/title
        | //def-list/def-item
        | //def-list/def-item/term
        | //def-list/def-item/def
        | //disp-formula
        | //disp-formula//mml:math
        | //disp-formula//mml:math//*
        | //disp-formula/label
        | //disp-formula-group/label
        | //disp-formula-group/caption
        | //disp-formula-group/caption/title
        | //fig
        | //fig/label
        | //fig/caption
        | //fig/caption/title
        | //fig-group
        | //fig-group/label
        | //fig-group/caption
        | //fig-group/caption/title
        | //list
        | //list/label
        | //list/title
        | //list/list-item
        | //p
        | //sec
        | //sec/label
        | //sec/title
        | //statement
        | //statement/label
        | //statement/title
        | //table
        | //table/thead
        | //table/tbody
        | //table/tfoot
        | //table//tr
        | //table//th
        | //table//td
        | //table-wrap
        | //table-wrap/label
        | //table-wrap/caption
        | //table-wrap/caption/title
        | //table-wrap-group
        | //table-wrap-group/label
        | //table-wrap-group/caption
        | //table-wrap-group/caption/title
        """
    dict = Dict(n => n for n in find(body,xpath))
    convert(Tree, body, dict)
end

function parse_back(back::EzXML.Node)
    tree = Tree(nodename(back))
    for node in find(back, ".//element-citation | .//mixed-citation")
        setnodename!(node, "citation")
    end
    xpath = """
        ref-list
        | ref-list/label
        | ref-list/title
        | ref-list/ref
        | ref-list/ref/citation
        | ref-list/ref/citation/article-title
        | ref-list/ref/citation/name
        | ref-list/ref/citation/collab
        | ref-list/ref/citation/day
        | ref-list/ref/citation/month
        | ref-list/ref/citation/year
        | ref-list/ref/citation/fpage
        | ref-list/ref/citation/lpage
        | ref-list/ref/citation/issue
        | ref-list/ref/citation/pub-id
        | ref-list/ref/citation/publisher-loc
        | ref-list/ref/citation/publisher-name
        | ref-list/ref/citation/source
        | ref-list/ref/citation/edition
        | ref-list/ref/citation/volume
        """
    dict = Dict(n => n for n in find(back,xpath))
    convert(Tree, back, dict)
end

function findfloats(tree::Tree)
    floatset = Set(["boxed-text","code","fig","fig-group","table-wrap","table-wrap-group"])
    floats = Tree[]
    topdown_while(tree) do t
        if t.name in floatset
            push!(floats, t)
            false
        else
            true
        end
    end
    floats
end

function postprocess!(tree::Tree)
    # merge <label> with <title>
    for i = length(tree)-1:-1:1
        if tree[i].name == "label" && tree[i+1].name == "title"
            prepend!(tree[i+1], tree[i].children)
            deleteat!(tree, i)
        end
    end
    foreach(postprocess!, tree.children)
end

max_stars_repo_name: paperai/LightJATS

max_issues_repo_name: paperai/LightJATS

max_forks_repo_name: paperai/LightJATS


text:
@testset "Testing Path Types" begin
    using Mapper2.MapperCore
    # Path Constructors
    let
        # Empty constructors
        ref = Path{Nothing}(String[])
        @test ref == Path{Nothing}()
        @test ref == Path{Nothing}("")

        # Various constructors.
        ref = Path{Nothing}(["a", "b", "c"])
        @test ref == Path{Nothing}("a.b.c")
        @test ref == Path{Nothing}("a", "b", "c")
        # Test fallback inequality
        @test Path{Nothing}("a.b.c") != Path{Int}("a.b.c")

        @test first(Path{Nothing}("a.b.c")) == "a"
        @test first(Path{Nothing}("a")) == "a"
        @test last(Path{Nothing}("a.b.c")) == "c"
        @test last(Path{Nothing}("c")) == "c"
    end
    # Length tests
    let
        @test length(Path{Nothing}("a.b.c")) == 3
        @test length(Path{Nothing}("a.b")) == 2
        @test length(Path{Nothing}("a")) == 1
        @test length(Path{Nothing}()) == 0
    end
    # Test dictionaries
    let
        paths = [Path{Nothing}("a.b.c"), Path{Int}("a.b.c"), Path{Nothing}("r")]

        d = Dict(paths[i] => i for i in 1:length(paths))

        for i in 1:length(paths)
            @test d[paths[i]] == i
        end
    end
end

@testset "Testing Ports" begin
    name = "test"
    meta = Dict("bob" => 10)
    p = Port(name, Input; metadata = meta)
    @test checkclass(p, MapperCore.Source) == true
    @test checkclass(p, MapperCore.Sink) == false

    pi = MapperCore.invert(p)
    @test pi.name == name
    @test pi.class == Output
    @test pi.metadata == meta

    # Error with incorrect class
    @test_throws Exception Port(name, :not_a_class)

    # Invert an output
    @test MapperCore.invert(MapperCore.invert(p)) == p
end

@testset "Testing Architecture Modeling" begin
    # Use Example1 for this testing.
    # Test General primitive
    let
        p = build_general_primitive()
        @test sort(collect(keys(p.ports))) == ["in[0]", "in[1]", "out[0]", "out[1]"]
        @test assert_no_children(p)
        @test assert_no_intrarouting(p)
    end
    let
        t = build_io_tile()
        @test assert_no_children(t) == false
        @test assert_no_intrarouting(t) == false
        # Do a recusrive search to find the correct task
        @test search_metadata!(t, "task", "input")
        @test search_metadata!(t, "task", "output")
    end

    expected_ports = sort([
        "north_out",
        "north_in",
        "south_out",
        "south_in",
        "east_out",
        "east_in",
        "west_out",
        "west_in",
    ])

    let
        t = build_general_tile()
        sorted_ports = sort(collect(keys(t.ports)))
        @test sorted_ports == expected_ports
        @test assert_no_children(t) == false
        @test assert_no_intrarouting(t) == false
        # Do a recusrive search to find the correct task
        @test search_metadata!(t, "task", "general")
    end
    let
        t = build_super_tile()
        sorted_ports = sort(collect(keys(t.ports)))
        @test sorted_ports == expected_ports
        @test assert_no_children(t) == false
        @test assert_no_intrarouting(t) == false
        # Do a recusrive search to find the correct task
        @test search_metadata!(t, "task", "general")
    end
    let
        t = build_double_general_tile()
        sorted_ports = sort(collect(keys(t.ports)))
        @test sorted_ports == expected_ports
        @test assert_no_children(t) == false
        @test assert_no_intrarouting(t) == false
        # Do a recusrive search to find the correct task
        @test search_metadata!(t, "task", "general")
    end
    let
        t = build_routing_tile()
        sorted_ports = sort(collect(keys(t.ports)))
        @test sorted_ports == expected_ports
        @test assert_no_children(t) == false
        @test assert_no_intrarouting(t) == false
        # Do a recusrive search to find the correct task
        @test search_metadata!(t, "task", "general") == false
    end
    let
        t = build_test_arch()
        # Do a recusrive search to find the correct task
        @test search_metadata!(t, "task", "general") == true
    end
end