text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/UMAT.jl/blob/master/LICENSE

using UMAT, Materials, Test, Tensors, DelimitedFiles

material = DruckerPragerMaterial(Sigma_t = 200.,
                               Sigma_c = 200.,
                               K_Zero = 100.,
                               C_const = 100.,
                               N_Power = 1.,
                               Eta = 5.,
                               A_Const = 200.,
                               B_Const = 0.,
                               E = 2.e5,
                               NU = 0.3,
                               DT = 0.,
                               PRNTYP = false,
                               FTOLER = 1.e-13,
                               MAXITER = 20,
                               ELNUM = 1,
                               IPNUM = 1,
                               MTYP = 2,
                               HARDTYP = 1,
                               DPC="DPC_LH1")

dtime = 0.25
stresses = [copy(tovoigt(material.variables.stress))]
dstrain11 = 1e-3*dtime
dtimes = [dtime, dtime, dtime, dtime, 1.0]
dstrains11 = [dstrain11, dstrain11, dstrain11, -dstrain11, -4*dstrain11]

expected_S11 = [50. 100. 150. 100. -100.]
expected_stress = zeros(6)

expected_strains = [[0.00025, -7.5e-5, -7.5e-5, 0.0, 0.0, 0.0],
                    [0.0005, -0.00015, -0.00015, 0.0, 0.0, 0.0],
                    [0.00075, -0.000225, -0.000225, 0.0, 0.0, 0.0],
                    [0.0005, -0.00015, -0.00015, 0.0, 0.0, 0.0],
                    [-0.0005, 0.00015, 0.00015, 0.0, 0.0, 0.0]]

for i in 1:length(dtimes)
    dstrain11 = dstrains11[i]
    dtime = dtimes[i]
    uniaxial_increment!(material, dstrain11, dtime)
    update_material!(material)
    expected_stress[1] = expected_S11[i]
    @test isapprox(tovoigt(material.variables.stress), expected_stress)
    @test isapprox(tovoigt(material.drivers.strain; offdiagscale=2.0), expected_strains[i])
end

#=
# Define strain history
e11 = 0.028.*vcat(Array(range(0, stop=1.5e-2, length=30)), Array(range(1.5e-2, stop=-1.5e-2, length=30)), Array(range(-1.5e-2, stop=3e-2, length=40)))
e22 = -e11/2
e33 = e22
strains = [[e11[i], e22[i], e33[i], 0.0, 0.0, 0.0] for i in 1:100]

s11 = [material.variables.STRESS[1]]
s22 = [material.variables.STRESS[2]]
s33 = [material.variables.STRESS[3]]

dtime = 0.01

for i=2:100
    dstrain = strains[i]-strains[i-1]
    material.ddrivers = UmatDriverState(NTENS=6, STRAN=dstrain)
    UMAT.integrate_material!(material)
    push!(s11, material.variables.STRESS[1])
    push!(s22, material.variables.STRESS[2])
    push!(s33, material.variables.STRESS[3])
    Materials.update_material!(material)
end

for sig in ["s11", "s22", "s33"]
    ref_stress = readdlm("test_druckerprager/ref_stresses_" * sig *".txt")
    for (comp,ref) in zip(eval(Symbol(sig)), ref_stress)
        @test isapprox(comp,ref,atol=sqrt(eps()))
    end
end
=#


text:
function CategoricalPool{T, R, V}(index::Vector{T},
                                  invindex::Dict{T, R},
                                  order::Vector{R},
                                  ordered::Bool) where {T, R, V}
    levels = similar(index)
    levels[order] = index
    pool = CategoricalPool{T, R, V}(index, invindex, order, levels, V[], ordered)
    buildvalues!(pool)
    return pool
end

function CategoricalPool(index::Vector{S},
                         invindex::Dict{S, T},
                         order::Vector{R},
                         ordered::Bool=false) where {S, T <: Integer, R <: Integer}
    invindex = convert(Dict{S, R}, invindex)
    C = catvaluetype(S, R)
    V = leveltype(C) # might be different from S (e.g. S == SubString, V == String)
    CategoricalPool{V, R, C}(index, invindex, order, ordered)
end

CategoricalPool{T, R, C}(ordered::Bool=false) where {T, R, C} =
    CategoricalPool{T, R, C}(T[], Dict{T, R}(), R[], ordered)
CategoricalPool{T, R}(ordered::Bool=false) where {T, R} =
    CategoricalPool(T[], Dict{T, R}(), R[], ordered)
CategoricalPool{T}(ordered::Bool=false) where {T} =
    CategoricalPool{T, DefaultRefType}(ordered)

function CategoricalPool{T, R}(index::Vector,
                               ordered::Bool=false) where {T, R}
    invindex = buildinvindex(index, R)
    order = Vector{R}(1:length(index))
    CategoricalPool(index, invindex, order, ordered)
end

function CategoricalPool(index::Vector, ordered::Bool=false)
    invindex = buildinvindex(index)
    order = Vector{DefaultRefType}(1:length(index))
    return CategoricalPool(index, invindex, order, ordered)
end

function CategoricalPool(invindex::Dict{S, R},
                         ordered::Bool=false) where {S, R <: Integer}
    index = buildindex(invindex)
    order = Vector{DefaultRefType}(1:length(index))
    return CategoricalPool(index, invindex, order, ordered)
end

# TODO: Add tests for this
function CategoricalPool(index::Vector{S},
                         invindex::Dict{S, R},
                         ordered::Bool=false) where {S, R <: Integer}
    order = Vector{DefaultRefType}(1:length(index))
    return CategoricalPool(index, invindex, order, ordered)
end

function CategoricalPool(index::Vector{T},
                         levels::Vector{T},
                         ordered::Bool=false) where {T}
    invindex = buildinvindex(index)
    order = buildorder(invindex, levels)
    return CategoricalPool(index, invindex, order, ordered)
end

function CategoricalPool(invindex::Dict{S, R},
                         levels::Vector{S},
                         ordered::Bool=false) where {S, R <: Integer}
    index = buildindex(invindex)
    order = buildorder(invindex, levels)
    return CategoricalPool(index, invindex, order, ordered)
end

Base.convert(::Type{T}, pool::T) where {T <: CategoricalPool} = pool

Base.convert(::Type{CategoricalPool{S}}, pool::CategoricalPool{T, R}) where {S, T, R <: Integer} =
    convert(CategoricalPool{S, R}, pool)

function Base.convert(::Type{CategoricalPool{S, R}}, pool::CategoricalPool) where {S, R <: Integer}
    if length(levels(pool)) > typemax(R)
        throw(LevelsException{S, R}(levels(pool)[typemax(R)+1:end]))
    end

    indexS = convert(Vector{S}, pool.index)
    invindexS = convert(Dict{S, R}, pool.invindex)
    order = convert(Vector{R}, pool.order)
    return CategoricalPool(indexS, invindexS, order, pool.ordered)
end

Base.copy(pool::CategoricalPool{T, R, V}) where {T, R, V} =
    CategoricalPool{T, R, V}(copy(pool.index), copy(pool.invindex), copy(pool.order),
                             copy(pool.levels), copy(pool.valindex), pool.ordered)

function Base.show(io::IO, pool::CategoricalPool{T, R}) where {T, R}
    @printf(io, "%s{%s,%s}([%s])", typeof(pool).name, T, R,
                join(map(repr, levels(pool)), ","))

    pool.ordered && print(io, " with ordered levels")
end

Base.length(pool::CategoricalPool) = length(pool.index)

Base.getindex(pool::CategoricalPool, i::Integer) = pool.valindex[i]
Base.get(pool::CategoricalPool, level::Any) = pool.invindex[level]
Base.get(pool::CategoricalPool, level::Any, default::Any) = get(pool.invindex, level, default)

"""
add the returned value to pool.invindex, this function doesn't do this itself to
avoid doing a dict lookup twice
"""
@inline function push_level!(pool::CategoricalPool{T, R}, level) where {T, R}
    x = convert(T, level)
    n = length(pool)
    if n >= typemax(R)
        throw(LevelsException{T, R}([level]))
    end

    i = R(n + 1)
    push!(pool.index, x)
    push!(pool.order, i)
    push!(pool.levels, x)
    push!(pool.valindex, catvalue(i, pool))
    i
end

function mergelevels(ordered, levels...)
    T = Base.promote_eltype(levels...)
    res = Vector{T}(undef, 0)

    nonempty_lv = Compat.findfirst(!isempty, levels)
    if nonempty_lv === nothing
        # no levels
        return res, ordered
    elseif all(l -> isempty(l) || l == levels[nonempty_lv], levels)
        # Fast path if all non-empty levels are equal
        append!(res, levels[nonempty_lv])
        return res, ordered
    end

    for l in levels
        levelsmap = indexin(l, res)

        i = length(res)+1
        for j = length(l):-1:1
            @static if VERSION >= v"0.7.0-DEV.3627"
                if levelsmap[j] === nothing
                    insert!(res, i, l[j])
                else
                    i = levelsmap[j]
                end
            else
                if levelsmap[j] == 0
                    insert!(res, i, l[j])
                else
                    i = levelsmap[j]
                end
            end
        end
    end

    # Check that result is ordered
    if ordered
        levelsmaps = [Compat.indexin(res, l) for l in levels]

        # Check that each original order is preserved
        for m in levelsmaps
            issorted(Iterators.filter(x -> x != nothing, m)) || return res, false
        end

        # Check that all order relations between pairs of subsequent elements
        # are defined in at least one set of original levels
        pairs = fill(false, length(res)-1)
        for m in levelsmaps
            @inbounds for i in eachindex(pairs)
                pairs[i] |= (m[i] != nothing) & (m[i+1] != nothing)
            end
            all(pairs) && return res, true
        end
    end

    res, false
end

@inline function Base.get!(pool::CategoricalPool, level::Any)
    get!(pool.invindex, level) do
        if isordered(pool)
            throw(OrderedLevelsException(level, pool.levels))
        end

        push_level!(pool, level)
    end
end

@inline function Base.get!(pool::CategoricalPool, level::CatValue)
    pool === level.pool && return level.level
    # Use invindex for O(1) lookup
    # TODO: use a global table to cache this information for all pairs of pools
    if level.pool.levels ⊈ keys(pool.invindex)
        if isordered(pool)
            throw(OrderedLevelsException(level, pool.levels))
        end
        newlevs, ordered = mergelevels(isordered(pool), pool.levels, level.pool.levels)
        # Exception: empty pool marked as ordered if new value is ordered
        if length(pool) == 0 && isordered(level.pool)
            ordered!(pool, true)
        end
        levels!(pool, newlevs)
    end
    get!(pool, get(level))
end

@inline function Base.push!(pool::CategoricalPool, level)
    get!(pool.invindex, level) do
        push_level!(pool, level)
    end
    return pool
end

# TODO: optimize for multiple additions
function Base.append!(pool::CategoricalPool, levels)
    for level in levels
        push!(pool, level)
    end
    return pool
end

function Base.delete!(pool::CategoricalPool{S}, levels...) where S
    for level in levels
        levelS = convert(S, level)
        if haskey(pool.invindex, levelS)
            ind = pool.invindex[levelS]
            delete!(pool.invindex, levelS)
            splice!(pool.index, ind)
            ord = splice!(pool.order, ind)
            splice!(pool.levels, ord)
            splice!(pool.valindex, ind)
            for i in ind:length(pool)
                pool.invindex[pool.index[i]] -= 1
                pool.valindex[i] = catvalue(i, pool)
            end
            for i in 1:length(pool)
                pool.order[i] > ord && (pool.order[i] -= 1)
            end
        end
    end
    return pool
end

function levels!(pool::CategoricalPool{S, R}, newlevels::Vector) where {S, R}
    levs = convert(Vector{S}, newlevels)
    if !allunique(levs)
        throw(ArgumentError(string("duplicated levels found in levs: ",
                                   join(unique(filter(x->sum(levs.==x)>1, levs)), ", "))))
    end

    n = length(levs)

    if n > typemax(R)
        throw(LevelsException{S, R}(setdiff(levs, levels(pool))[typemax(R)-length(levels(pool))+1:end]))
    end

    # No deletions: can preserve position of existing levels
    # equivalent to issubset but faster due to JuliaLang/julia#24624
    if isempty(setdiff(pool.index, levs))
        append!(pool, setdiff(levs, pool.index))
    else
        empty!(pool.invindex)
        resize!(pool.index, n)
        resize!(pool.valindex, n)
        resize!(pool.order, n)
        resize!(pool.levels, n)
        for i in 1:n
            v = levs[i]
            pool.index[i] = v
            pool.invindex[v] = i
            pool.valindex[i] = catvalue(i, pool)
        end
    end

    buildorder!(pool.order, pool.invindex, levs)
    for (i, x) in enumerate(pool.order)
        pool.levels[x] = pool.index[i]
    end
    return pool
end

index(pool::CategoricalPool) = pool.index
DataAPI.levels(pool::CategoricalPool) = pool.levels
order(pool::CategoricalPool) = pool.order

isordered(pool::CategoricalPool) = pool.ordered
ordered!(pool::CategoricalPool, ordered) = (pool.ordered = ordered; pool)

# LevelsException
function Base.showerror(io::IO, err::LevelsException{T, R}) where {T, R}
    levs = join(repr.(err.levels), ", ", " and ")
    print(io, "cannot store level(s) $levs since reference type $R can only hold $(typemax(R)) levels. Use the decompress function to make room for more levels.")
end

# OrderedLevelsException
function Base.showerror(io::IO, err::OrderedLevelsException)
    print(io, "cannot add new level $(err.newlevel) since ordered pools cannot be extended implicitly. Use the levels! function to set new levels, or the ordered! function to mark the pool as unordered.")
end

ext: jl


text:
B_data = [  -2.5           0.5           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   3.0          -2.0           1.0            0.0   0.0   0.0   0.0     0.0   0.0    0.0          0.0           0.0            0.0            0.0            0.0           0.0       0.0        0.0       0.0        0.0
1.0          -3.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   3.0          -2.0           1.0            0.0   0.0   0.0   0.0     0.0   0.0    0.0          0.0           0.0            0.0            0.0            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0          -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0           0.0           0.0            0.0   0.0   0.0   0.0     0.0   0.0    0.5          0.0           0.0            0.0            0.0            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0          -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           2.0           0.0            0.0   0.0   0.0   0.0     0.0   0.0    0.0          0.5           0.0            0.0            0.0            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           2.0            0.0   0.0   0.0   0.0     0.0   0.0    0.0          0.0           0.5            0.0            0.0            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0            2.0   0.0   0.0   0.0     0.0   0.0    0.0          0.0           0.0            0.5            0.0            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0            0.0   2.0   0.0   0.0     0.0   0.0    0.0          0.0           0.0            0.0            0.5            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0            0.0   0.0   2.0   0.0     0.0   0.0    0.0          0.0           0.0            0.0            0.0            0.5           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0            0.0   0.0   0.0   2.0     0.0   0.0    0.0          0.0           0.0            0.0            0.0            0.0           0.5       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0            0.0   0.0   0.0   0.0     2.0   0.0    0.0          0.0           0.0            0.0            0.0            0.0           0.0       0.5        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0            0.0   0.0   0.0   0.0     0.0   2.0    0.0          0.0           0.0            0.0            0.0            0.0           0.0       0.0        0.5       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0           0.0           0.0            0.0   0.0   0.0   0.0     0.0   0.0    0.5          0.0           0.0            0.0            0.0            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           2.0           0.0            0.0   0.0   0.0   0.0     0.0   0.0    0.0          0.5           0.0            0.0            0.0            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           2.0            0.0   0.0   0.0   0.0     0.0   0.0    0.0          0.0           0.5            0.0            0.0            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0           0.0           0.0            2.0   0.0   0.0   0.0     0.0   0.0    0.0          0.0           0.0            0.5            0.0            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0           0.0           0.0            0.0   2.0   0.0   0.0     0.0   0.0    0.0          0.0           0.0            0.0            0.5            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0           0.0           0.0            0.0   0.0   2.0   0.0     0.0   0.0    0.0          0.0           0.0            0.0            0.0            0.5           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0           0.0           0.0            0.0   0.0   0.0   2.0     0.0   0.0    0.0          0.0           0.0            0.0            0.0            0.0           0.5       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0           0.0           0.0            0.0   0.0   0.0   0.0     2.0   0.0    0.0          0.0           0.0            0.0            0.0            0.0           0.0       0.5        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0           0.0           0.0            0.0   0.0   0.0   0.0     0.0   2.0    0.0          0.0           0.0            0.0            0.0            0.0           0.0       0.0        0.5       0.0        0.0
0.0           0.0           1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -6.0           6.66667      -2.66667        0.0   0.0   0.0   0.0     0.0   0.0    1.0          0.0           0.0            0.0            0.0            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.666667     -2.0           0.666667       0.0   0.0   0.0   0.0     0.0   0.0    0.0          1.0           0.0            0.0            0.0            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.666667     -2.66667      -4.0            6.0  -4.0   2.0   0.0     0.0   0.0    0.0          0.0           1.0            0.0            0.0            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0          -14.0  20.0  -8.0   0.0     0.0   0.0    0.0          0.0           0.0            1.0            0.0            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0           -2.0  -2.0   2.0   0.0     0.0   0.0    0.0          0.0           0.0            0.0            1.0            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0            2.0  -8.0  -8.0  18.0   -12.0   6.0    0.0          0.0           0.0            0.0            0.0            1.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0            0.0   0.0   0.0  -3.5     2.5  -1.0    0.0          0.0           0.0            0.0            0.0            0.0           1.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0            0.0   0.0   0.0  -0.25   -2.0   0.25   0.0          0.0           0.0            0.0            0.0            0.0           0.0       1.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0            0.0   0.0   0.0   0.25   -1.0  -2.75   0.0          0.0           0.0            0.0            0.0            0.0           0.10592  -0.211839   1.31776   0.218262   1.0699
4.35288       0.910678      0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   3.10467      -0.736445      0.368223       0.0   0.0   0.0   0.0     0.0   0.0   -6.0         -4.0           1.0            0.0            0.0            0.0           0.0       0.0        0.0       0.0        0.0
1.81225e-16   3.79146e-17   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   4.5991e-17    2.0           1.53303e-17    0.0   0.0   0.0   0.0     0.0   0.0    1.0         -3.0          -1.0            0.0            0.0            0.0           0.0       0.0        0.0       0.0        0.0
-1.20817e-15  -2.52764e-16   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0  -3.06607e-16   2.04405e-16   2.0            0.0   0.0   0.0   0.0     0.0   0.0   -4.0         10.0          -9.0            0.0            0.0            0.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0            2.0   0.0   0.0   0.0     0.0   0.0    9.0        -18.0          27.0          -12.0          -12.0            3.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0            0.0   2.0   0.0   0.0     0.0   0.0    3.747e-16   -7.49401e-16   1.1241e-15     3.0           -3.0           -3.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0           0.0           0.0            0.0   0.0   2.0   0.0     0.0   0.0   -2.498e-15    4.996e-15    -7.49401e-15  -12.0           30.0          -21.0           0.0       0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0           0.0           0.0            0.0   0.0   0.0   2.0     0.0   0.0    0.0          0.0           0.0            1.125         -2.25           3.375        -4.125    -1.5        0.375     0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0           0.0           0.0            0.0   0.0   0.0   0.0     2.0   0.0    0.0          0.0           0.0            4.68375e-17   -9.36751e-17    1.40513e-16   0.375    -3.0       -0.375     0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0           0.0           0.0            0.0   0.0   0.0   0.0     0.0   2.0    0.0          0.0           0.0           -3.1225e-16     6.245e-16     -9.36751e-16  -1.5       3.75      -5.25      0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0            0.0   0.0   0.0   0.0     0.0   0.0    0.0          0.0           0.0            0.0            0.0            0.0           0.25     -0.5        0.75     -2.5        2.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0            0.0   0.0   0.0   0.0     0.0   0.0    0.0          0.0           0.0            0.0            0.0            0.0           0.5      -1.0        1.5       1.0       -2.0]


text:
export Astrobee3D

mutable struct Astrobee3D{T<:AbstractFloat} <: Robot
  mass::T
  J::Matrix{T}
  Jinv::Matrix{T}
  n_thrusters::Int
  r::T
  hard_limit_vel::T
  hard_limit_accel::T
  hard_limit_ω::T
  hard_limit_α::T
  btCollisionObject
end
function Astrobee3D{T}() where T
  n_thrusters = 12
  s = 0.5*0.305   # each side of cube is 30.5cm
  r = sqrt(3)*s   # inflate to sphere

  # flight robot param: freeflyer/astrobee/config/world/iss.config
  mass = 7.0
  hard_limit_vel = 0.5      # Actual limit: 0.5
  hard_limit_accel = 0.1
  hard_limit_ω = 45*π/180
  hard_limit_α = 50*π/180

  J = 0.1083*Eye(3)
  Jinv = inv(J)

  btCollisionObject = BulletCollision.sphere(SVector{3}(zeros(T,3)), r)

  # new astrobee instance
  return Astrobee3D{T}(mass, J, Jinv, n_thrusters, r, hard_limit_vel, hard_limit_accel, hard_limit_ω, hard_limit_α,btCollisionObject)
end
Astrobee3D(::Type{T} = Float64; kwargs...) where {T} = Astrobee3D{T}(; kwargs...)


text:
abstract type Analysis end
abstract type CallAnalysis <: Analysis end

struct StaticAnalysis <: CallAnalysis
    reach::Dict
    sites::Vector{Variable}
    addrs::Vector{QuoteNode}
    map::Dict
end

struct FallbackAnalysis <: CallAnalysis
end

mutable struct CallGraph <: Analysis
    dependencies::Dict{Address, Set{Address}}
    CallGraph(map) = new(map)
    CallGraph() = new(Dict{Address, Set{Address}}())
end

# ---- Analysis ---- #

function lower_to_ir(call, type...)
    sig = Tuple{typeof(call), type...}
    m = meta(sig)
    ir = IR(m)
    return ir
end

function control_flow_check(ir)
    length(ir.blocks) > 1 && return false
    return true
end

# Reaching analysis.
function reaching!(reach::Vector{Variable}, p, var, ir)
    for (v, st) in ir
        st.expr isa Expr && begin
            if var in st.expr.args
                push!(reach, v)
                reaching!(reach, p, v, ir)
            end
        end
    end
end

function reaching(var::Variable, ir)
    reach = Variable[]
    for (v, st) in ir
        st.expr isa Expr && begin
            if var in st.expr.args
                push!(reach, v)
                reaching!(reach, var, v, ir)
            end
        end
    end
    return reach
end

function flow_analysis(ir)
    sites = Variable[]
    addrs = QuoteNode[]
    var_sym_map = Dict{Variable, QuoteNode}()
    reach = Dict{Variable, Any}()
    for (v, st) in ir
        MacroTools.postwalk(st) do e
            @capture(e, call_(sym_, args__))
            call isa GlobalRef && call.name == :rand && begin
                push!(sites, v)
                push!(addrs, sym)
                var_sym_map[v] = sym
                r = reaching(v, ir)
                !isempty(r) && begin
                    reach[v] = r
                end
            end
            e
        end
    end
    return StaticAnalysis(reach, sites, addrs, var_sym_map)
end

function dependency(a::Analysis)
    map = Dict{Symbol, Set{Symbol}}()
    for (k, vars) in a.reach
        depends = Symbol[]
        for v in vars
            haskey(a.map, v) && push!(depends, a.map[v].value)
        end
        map[a.map[k].value] = Set(depends)
    end
    return CallGraph(map)
end

# ---- Driver ---- #

# Returns the dependency analysis in call graph (tree) form.
function construct_graph!(parent, addr, call, type)
    ir = lower_to_ir(call, type)
    if control_flow_check(ir)
        graph = CallGraph()
    else
        analysis = flow_analysis(ir)
        dependencies = dependency(analysis)
        graph = CallGraph(dependencies)
    end
    parent[addr] = graph
end

# Toplevel analysis driver.
function construct_graph(call, type)
    ir = lower_to_ir(call, type)
    if !control_flow_check(ir)
        graph = CallGraph()
    else
        analysis = flow_analysis(ir)
        dependencies = dependency(analysis)
        graph = dependencies
    end
    return graph
end

function construct_graph(ir)
    if !control_flow_check(ir)
        graph = CallGraph()
    else
        analysis = flow_analysis(ir)
        dependencies = dependency(analysis)
        graph = dependencies
    end
    return graph
end


text:
using BlossomV

if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

@testset "Trivial Case" begin
    differences = [
        0 1 500
    ]

    m = Matching(2) # defaults to Int32
    for row_ii in 1:size(differences,1)
        n1,n2, c = differences[row_ii,:]
        add_edge(m,n1,n2,c)
    end
    solve(m)
    @test get_match(m,0) == 1
    @test get_match(m,1) == 0

    differences = [
        0 1 1.3
    ]

    m = Matching(Float64, 2)
    for row_ii in 1:size(differences,1)
        n1,n2, c = differences[row_ii,:]
        add_edge(m,Int(n1),Int(n2),c)
    end
    solve(m)
    @test get_match(m,0) == 1
    @test get_match(m,1) == 0
end

@testset "Base Case -- greedy will work" begin
    differences = [
        0 1 500
        0 2 600
        1 2 700
        2 3 100
        1 3 1000
    ]

    m = Matching(Int, 4)
    for row_ii in 1:size(differences,1)
        n1,n2, c = differences[row_ii,:]
        add_edge(m,n1,n2,c)
    end
    solve(m)
    @test get_match(m,0) == 1
    @test get_match(m,1) == 0
    @test get_match(m,2) == 3
    @test get_match(m,3) == 2
end

@testset "Intermediate Case -- no greedy solution" begin
    differences = [
        0 1 500
        0 2 400
        1 2 300
        2 3 1000
        1 3 1000
    ]

    m = Matching(Int, 4)
    for row_ii in 1:size(differences,1)
        n1,n2, c = differences[row_ii,:]
        add_edge(m,Int(n1),Int(n2),c)
    end
    solve(m)
    @test get_match(m,0) == 2
    @test get_match(m,1) == 3
    @test get_match(m,2) == 0
    @test get_match(m,3) == 1
end

@testset "Error handling" begin
    m = Matching(Int32, 2)
    @test_throws ErrorException add_edge(m, 0, 0, 1)
    @test_throws ErrorException add_edge(m, -1, 0, 1)
    @test_throws ErrorException add_edge(m, 0, -1, 1)
    @test_throws ErrorException add_edge(m, 0, 1, -1)
end


text:
using DataPackage
using Base.Test

import DataPackage: get_resource, read

include("read.jl")


text:
function eventloop(socket)
    task_local_storage(:IJulia_task, "write task")
    try
        while true
            msg = recv_ipython(socket)
            try
                send_status("busy", msg)
                invokelatest(handlers[msg.header["msg_type"]], socket, msg)
            catch e
                # Try to keep going if we get an exception, but
                # send the exception traceback to the front-ends.
                # (Ignore SIGINT since this may just be a user-requested
                #  kernel interruption to interrupt long calculations.)
                if !isa(e, InterruptException)
                    content = error_content(e, msg="KERNEL EXCEPTION")
                    map(s -> println(orig_stderr[], s), content["traceback"])
                    send_ipython(publish[], msg_pub(execute_msg, "error", content))
                end
            finally
                # @async no longer implicitly captures local variables (see PR #618)
                let msg = msg
                    @async begin
                        flush_all()
                        send_status("idle", msg)
                    end
                end
            end
        end
    catch e
        # the Jupyter manager may send us a SIGINT if the user
        # chooses to interrupt the kernel; don't crash on this
        if isa(e, InterruptException)
            eventloop(socket)
        else
            rethrow()
        end
    end
end

function waitloop()
    @async eventloop(control[])
    requests_task = @async eventloop(requests[])
    while true
        try
            wait()
        catch e
            # send interrupts (user SIGINT) to the code-execution task
            if isa(e, InterruptException)
                @async Base.throwto(requests_task, e)
            else
                rethrow()
            end
        end
    end
end


text:
"""
A tabular data type.

Fields:

* `index::Dict{Symbol, Int}`: A mapping from symbols to numeric column indices
* `columns::Vector{NullableVector}`: A vector of NullableVector columns
* `fields::Vector{Symbol}`: A vector of fields (column names)
<!-- * `allowsnulls::Vector{Bool}`: A vector of flags indicating whether the
respective column allows storage of null values -->
* `hasnulls::Vector{Bool}`: A vector of flags indicating whether the
respective column currently stores any null values

Notes: The ordering of the columns in `columns` respects that of the numeric
indices given by `index`.

<!-- If `allowsnulls[i] == false` then `hasnulls[i]` must also `== false`.

Columns are assumed to be either `Array` or `NullableArray` objects. If
`isa(columns[i], Array)` then `allowsnulls[i] == false`. -->

The types of columns are not coerced upon `Table` initialization.

`'.'` is not a valid character within a field. All instances of `'.'` will be
replaced by `'_'` upon constructing a `Table` or setting a column.
"""
type Table <: AbstractTables.AbstractTable
    index::Dict{Symbol, Int}
    columns::Vector{NullableVector}
    fields::Vector{Symbol}
    # allowsnulls::Vector{Bool}
    hasnulls::Vector{Bool}
    metadata::Dict{Symbol, Any}

    function Table(index, columns)
        ncols = length(columns)
        # allowsnulls = Vector{Bool}(ncols)
        hasnulls = Vector{Bool}(ncols)
        if ncols > 1
            nrows = length(columns[1])
            equallengths = true
            for (j, col) in enumerate(columns)
                equallengths &= length(col) == nrows
                # allowsnulls[j] = _allowsnulls(col)
                hasnulls[j] = _hasnulls(col)
            end
            if !equallengths
                msg = "All columns in a Table must be the same length"
                throw(ArgumentError(msg))
            end
        end
        fields = Array{Symbol}(ncols)
        for key in keys(index)
            key = Symbol(replace_dots(string(key)))
            fields[index[key]] = key
        end
        length(index) == length(columns) || error()
        new(index, columns, fields, hasnulls, Dict{Symbol, Any}())
    end
end


text:
LinAlg.Ac_ldiv_B!{T}(D::Diagonal{T}, B) = A_ldiv_B!(D, B)

function LinAlg.A_ldiv_B!{T}(D::Diagonal{T}, B::Diagonal{T})
    if size(D) ≠ size(B)
        throw(DimensionMismatch("size(D) ≠ size(B)"))
    end
    map!(/, B.diag, B.diag, D.diag)
    B
end

function LinAlg.A_ldiv_B!{T}(D::Diagonal{T}, B::SparseMatrixCSC{T})
    dd = D.diag
    if length(dd) ≠ size(B, 1)
        throw(DimensionMismatch("size(D,2) ≠ size(B,1)"))
    end
    vals = nonzeros(B)
    rows = rowvals(B)
    @inbounds for k in eachindex(vals)
        vals[k] /= dd[rows[k]]
    end
    B
end

function rowlengths(L::LowerTriangular)
    ld = L.data
    [(sl = Compat.view(ld, i, 1:i); sqrt(dot(sl, sl))) for i in 1:size(L, 1)]
end


text:
module PolyLib
	const polylib = "libpolylib64"

	include("common.jl")
	include("functions.jl")
	include("extra.jl")
end


text:
########################################
# Timezone utilities
########################################

function isfixedtz(s)
    match(r"^UTC[+-][0-9]{1,2}$", s) !== nothing
end

function getsendertz(db, sender_id)
    sender = select(db, Vector{Sender}, (:id => sender_id, ))
    isempty(sender) ? "UTC+0" : sender[1].tz
end

function process_timezone(obj::ZulipRequest, db, opts)
    @debug "timezone"

    m = match(r"timezone\s*(.*)", obj.data)
    if isempty(m[1])
        # Show stored value
        tz = getsendertz(db, obj.message.sender_id)
        return "Your timezone is $tz"
    else
        if m[1] in TZS || isfixedtz(m[1])
            upsert(db, Sender(obj.message.sender_id, m[1]), (:tz, ))
            return "New timezone saved"
        else
            return "Unknown timezone format. Please refer [TimeZones.jl documentation](https://juliatime.github.io/TimeZones.jl/stable/types/#TimeZone-1) for the list of available timezones"
        end
    end
end

function process_remove(obj::ZulipRequest, db, channel, ts, opts)
    @debug "remove"
    m = match(r"remove\s+(.*)", obj.data)
    m === nothing && return "No scheduled messages were removed"
    ids = strip.(split(m[1], ","))
    for str_id in ids
        id = tryparse(Int, str_id)
        id === nothing && continue
        delete(db, TimedMessage, (:id => id, :msg_sender_id => obj.message.sender_id))
        tmsg = TimedMessage(id, ts, ts, Message("", "", "", obj.message.sender_id, ""))
        put!(channel, (tmsg, -1))
    end
    return "Messages were removed from the schedule"
end

function process_list(obj::ZulipRequest, db, opts)
    @debug "list"
    tmsgs = select(db, Vector{TimedMessage}, (:msg_sender_id => obj.message.sender_id, ))
    isempty(tmsgs) && return "No messages scheduled"
    sort!(tmsgs, by = x -> x.exects)

    iob = IOBuffer()
    tz0 = getsendertz(db, obj.message.sender_id)
    tz = tz0 in TZS ? TimeZone(tz0) : FixedTimeZone(tz0)
    isfirst = true
    for tmsg in tmsgs
        exects = astimezone(ZonedDateTime(unix2datetime(round(Int, tmsg.exects/1000)), TimeZone("GMT")), tz)
        if !isfirst
            print(iob, "\n---\n")
        end
        print(iob, "**id:** ", tmsg.id, "\n")
        print(iob, "**scheduled:** ", Dates.format(exects, "yyyy-mm-dd HH:MM:SS z"), "\n")
        print(iob, String(base64decode(tmsg.msg.content)))
        isfirst = false
    end

    return String(take!(iob))
end

function process_help(obj::ZulipRequest, db, opts)
    @debug "help"
    return """
Currently following commands are supported

1. `<where> <when> <what>`: set a reminder.
    - `where` is optional, should be either `me` or `here`. In latter case bot sends a message to the topic where reminder was set. If `where` is omitted reminder bot sends message privately to the person who set the reminder.
    - `when` can be either in relative form `X days Y hours Z minutes` or in an absolute form `2020-10-01 23:15:00`. In relative forms single or plural form of `month`, `week`, `day`, `hour`, `minute`, `second` are allowed as well as positive and negative values of `X`, `Y`, `Z`. In absolute form date is mandatory, but hours, minute or second part can be omitted.
    - `what` is a message that should be shown by reminder bot.
2. `list`: show all current reminders of a user.
3. `remove <id>`: remove your reminder with the id `<id>`. Multiple `<id>` can be given comma separated.
4. `timezone <value>`: set timezone for current user. If `<value>` is omitted, then current setting is used. Value should be in a form `Europe/Amsterdam`, `America/New_York` and the like.
5. `help`: this message

Examples of usage:
- `me 2 days drink coffee` (send private message in two hours)
- `1 day -2 hour redo things` (send message tomorrow, two hours earlier than today's current time)
- `2021-12-31 12:00 Happy New Year` (send private message on the midnight of 31 December 2021)
- `here 1 hour Say something` (send message to the stream in 1 hour)
    """
end

function process_reminder(obj::ZulipRequest, db, channel, ts, opts)
    @debug "reminder"
    gde, status, msg, exects = zparse(obj.data, ts)
    status == :unknown && return "Unable to process message. Please refer to `help` for the list of possible commands."

    tz0 = getsendertz(db, obj.message.sender_id)
    tz = tz0 in TZS ? TimeZone(tz0) : FixedTimeZone(tz0)
    if status == :absolute
        exects0 = ZonedDateTime(exects, tz)
        exects = astimezone(exects0, localzone())
    else
        exects = ZonedDateTime(exects, localzone())
        exects0 = astimezone(exects, tz)
    end
    exects = DateTime(exects)

    content = ((gde == :here) & (obj.message.type == "stream")) ? "On behalf of @**$(obj.message.sender_full_name)**\n" : ""
    content *= if obj.message.type == "stream"
        startswithnarrow(msg) ? "" : (narrow(obj.message) * "\n")
    else
        ""
    end
    content *= gde == :me ? "**Reminder**: " * msg : msg

    content = base64encode(content)
    msg = if obj.message.type == "stream" && gde == :here
        stream = base64encode(obj.message.display_recipient)
        topic = base64encode(obj.message.subject)
        Message(stream, topic, "stream", obj.message.sender_id, content)
    else
        Message("", "", "private", obj.message.sender_id, content)
    end

    tmsg = TimedMessage(ts, exects, msg)
    tmsg = insert(db, tmsg)
    @debug tmsg
    put!(channel, (tmsg, 1))
    "Message is scheduled on $(Dates.format(exects0, "yyyy-mm-dd HH:MM:SS z"))"
end


text:
module SpringCollab2020CustomSandwichLavaSettingsTrigger

using ..Ahorn, Maple

@mapdef Trigger "SpringCollab2020/CustomSandwichLavaSettingsTrigger" CustomSandwichLavaSettingsTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
    onlyOnce::Bool=false, direction::String="CoreModeBased", speed::Number=20.0)

const directions = String["AlwaysUp", "AlwaysDown", "CoreModeBased"]

const placements = Ahorn.PlacementDict(
    "Custom Sandwich Lava Settings (Spring Collab 2020)" => Ahorn.EntityPlacement(
        CustomSandwichLavaSettingsTrigger,
        "rectangle",
    ),
)

Ahorn.editingOptions(entity::CustomSandwichLavaSettingsTrigger) = Dict{String, Any}(
    "direction" => directions
)

end


text:
import Statistics: cov, mean
import LinearAlgebra: inv, det, norm
import MultivariateStats

"""
    circlepoints(r, n)

Generate `n` evenly spaced points on circle of radius `r`
"""
function circlepoints(n, r; noise = 0.0)
    t = LinRange(0.0,2π,n)
    x = r .* cos.(t) .+ rand(n).*noise
    y = r .* sin.(t) .+ rand(n).*noise
    return x, y
end

function gaussian(X::AbstractMatrix{<:Real}; dims=2)
    k = size(X, 1)
    Σ = cov(X, dims=dims)
    Σ⁻¹ = inv(Σ)
    μ = mean(X, dims=dims)
    Xc = X.-μ
    exp.(-0.5*Xc'*Σ⁻¹*Xc)./sqrt(det(Σ)*(2π)^k)
end

function density(X::AbstractMatrix{<:Real}; ε=1.0)
    D = MultivariateStats.pairwise((x,y)->exp(-(1/ε)*norm(x-y)^2.0), X)
    SD = sum(D, dims=2)
    vec(SD./sum(SD))
end

function eccentricity(X::AbstractMatrix{<:Real}; p=2.0)
    N = size(X, 2)
    D = MultivariateStats.pairwise((x,y)->norm(x-y,p), X)
    vec(sum(D, dims=2)./N).^(1/p)
end

"""
    balanced(x[; intervals=10, overlap=0.2])

Balanced cover with `intervals` and `overlap` fraction.
The interval boundaries are distributed so that each patch covers the same fraction of the data set.
"""
function balanced(x::AbstractVector{<:Real}; kwargs...)

    # setup parameters
    intervals=10
    overlap=0.5
    for (p,v) in kwargs
        if p == :intervals
            intervals = v
        elseif p == :overlap
            overlap = v
        end
    end

    xmin, xmax = extrema(x)
    xrng = xmax - xmin
    ilen = 1.0 / (intervals - (intervals-1)*overlap)
    istep = ilen*(1-overlap)

    # compute ranges of intervals with an overlap
    irngs = [let rmin = (i-1)*istep;
                (rmin, (rmin+ilen)) .* xrng .+ xmin
            end
            for i in 1:intervals]

    # generate indexes for cover elements
    patches = map(i->findall(e->i[1]<=e<=i[2], x), irngs)
    return patches
end


text:
import LinearAlgebra.generic_lufact!

const lu_mt_threshold = Ref(64.0)
mt_thresholds[:lu] = lu_mt_threshold;

function generic_lufact!(A::StridedMatrix{DoubleFloat{T}}, ::Val{Pivot} = Val(true);
                         check::Bool = true) where {T<:AbstractFloat,Pivot}
    m, n = size(A)
    minmn = min(m,n)
    info = 0
    ipiv = Vector{BlasInt}(undef, minmn)
    use_threads = (nthreads() > 1) &&
        (Float64(m)*Float64(n) > lu_mt_threshold[])
    liA = LinearIndices(A)

    @inbounds begin
        for k = 1:minmn
            # find index max
            kp = k
            if Pivot
                amax = abs(zero(T))
                for i = k:m
                    absi = abs(A[i,k])
                    if absi > amax
                        kp = i
                        amax = absi
                    end
                end
            end
            ipiv[k] = kp
            if !iszero(A[kp,k])
                if k != kp
                    # Interchange
                    for i = 1:n
                        tmp = A[k,i]
                        A[k,i] = A[kp,i]
                        A[kp,i] = tmp
                    end
                end
                # Scale first column
                Akkinv = inv(A[k,k])
                for i = k+1:m
                    A[i,k] *= Akkinv
                end
            elseif info == 0
                info = k
            end
            # Update the rest
            if use_threads
                _mt_lu_loop(m,n,k,A,Vec{Npref,T})
            else
                for j = k+1:n
                    #axpy!(-A[k,j],view(A,k+1:m,k),view(A,k+1:m,j),Vec{Npref,T})
                    _axpy!(m-k,-A[k,j],A,liA[k+1,k],A,liA[k+1,j],Vec{Npref,T})
                end
            end
        end
    end
    check && checknonsingular(info)
    return LU{DoubleFloat{T},typeof(A)}(A, ipiv, convert(BlasInt, info))
end

# Someday @threads will not stupefy the inference engine.

@noinline function _mt_lu_loop(m,n,k,A,::Type{Vec{N,T}}) where {N,T}
    liA = LinearIndices(A)
    @threads for j = k+1:n
        @inbounds begin
            _axpy!(m-k,-A[k,j],A,liA[k+1,k],A,liA[k+1,j],Vec{N,T})
        end
    end
end

max_stars_repo_name: RalphAS/DoubleBLAS


text:
import PredictMD
import Test

a = Union{String, Missing}["foo", "bar", "foo"]
b = PredictMD.disallowmissing(a)
Test.@test isa(b, Vector{String})
Test.@test b == String["foo", "bar", "foo"]
c = Union{String, Missing}["foo", "bar", missing]
Test.@test_throws MethodError PredictMD.disallowmissing(c)


text:
include("temperature_render.jl")

# angles = 75.0:1.0:85.0
# angles = [-179.0:1.0:-1.0;1.0:1.0:179.0]
masses = 1.0:5.0:200
n_frames = length(masses)
anim = @animate for (i, mass) in enumerate(masses)
    print("Frame: $i / $n_frames complete")
    # fangle = Float64(angle)
    fmass = Float64(mass)
    temperature_render(
        mass = mass,
        tolerance = 1e-10,
        resolution = 1080,
        size_multiplier = 4,
        dtmax = 2.5,
    )
end
gif(anim, "massgif.gif", fps = 20)


text:
module EvaluateEquationModule

import ..CoreModule: Node, Options
import ..UtilsModule: @return_on_false, is_bad_array

macro return_on_check(val, T, n)
    # This will generate the following code:
    # if !isfinite(val)
    #     return (Array{T, 1}(undef, n), false)
    # end

    :(
        if !isfinite($(esc(val)))
            return (Array{$(esc(T)),1}(undef, $(esc(n))), false)
        end
    )
end

macro return_on_nonfinite_array(array, T, n)
    :(
        if is_bad_array($(esc(array)))
            return (Array{$(esc(T)),1}(undef, $(esc(n))), false)
        end
    )
end

"""
    eval_tree_array(tree::Node, cX::AbstractMatrix{T}, options::Options)

Evaluate a binary tree (equation) over a given input data matrix. The
options contain all of the operators used. This function fuses doublets
and triplets of operations for lower memory usage.

This function can be represented by the following pseudocode:

```
function eval(current_node)
    if current_node is leaf
        return current_node.value
    elif current_node is degree 1
        return current_node.operator(eval(current_node.left_child))
    else
        return current_node.operator(eval(current_node.left_child), eval(current_node.right_child))
```
The bulk of the code is for optimizations and pre-emptive NaN/Inf checks,
which speed up evaluation significantly.

# Returns

- `(output, complete)::Tuple{AbstractVector{T}, Bool}`: the result,
    which is a 1D array, as well as if the evaluation completed
    successfully (true/false). A `false` complete means an infinity
    or nan was encountered, and a large loss should be assigned
    to the equation.
"""
function eval_tree_array(
    tree::Node, cX::AbstractMatrix{T}, options::Options
)::Tuple{AbstractVector{T},Bool} where {T<:Real}
    n = size(cX, 2)
    result, finished = _eval_tree_array(tree, cX, options)
    @return_on_false finished result
    @return_on_nonfinite_array result T n
    return result, finished
end

function _eval_tree_array(
    tree::Node, cX::AbstractMatrix{T}, options::Options
)::Tuple{AbstractVector{T},Bool} where {T<:Real}
    if tree.degree == 0
        deg0_eval(tree, cX, options)
    elseif tree.degree == 1
        # TODO: We could all do Val(tree.l.degree) here, instead of having
        # different kernels for const vs data.

        # We fuse (and compile) the following:
        #  - op(op2(x, y)), where x, y, z are constants or variables.
        #  - op(op2(x)), where x is a constant or variable.
        #  - op(x), for any x.
        if tree.l.degree == 2 && tree.l.l.degree == 0 && tree.l.r.degree == 0
            deg1_l2_ll0_lr0_eval(tree, cX, Val(tree.op), Val(tree.l.op), options)
        elseif tree.l.degree == 1 && tree.l.l.degree == 0
            deg1_l1_ll0_eval(tree, cX, Val(tree.op), Val(tree.l.op), options)
        else
            deg1_eval(tree, cX, Val(tree.op), options)
        end
    else
        # We fuse (and compile) the following:
        #  - op(x, y), where x, y are constants or variables.
        #  - op(x, y), where x is a constant or variable but y is not.
        #  - op(x, y), where y is a constant or variable but x is not.
        #  - op(x, y), for any x or y
        # TODO - add op(op2(x, y), z) and op(x, op2(y, z))
        if tree.l.degree == 0 && tree.r.degree == 0
            deg2_l0_r0_eval(tree, cX, Val(tree.op), options)
        elseif tree.l.degree == 0
            deg2_l0_eval(tree, cX, Val(tree.op), options)
        elseif tree.r.degree == 0
            deg2_r0_eval(tree, cX, Val(tree.op), options)
        else
            deg2_eval(tree, cX, Val(tree.op), options)
        end
    end
end

function deg2_eval(
    tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, options::Options
)::Tuple{AbstractVector{T},Bool} where {T<:Real,op_idx}
    n = size(cX, 2)
    (cumulator, complete) = _eval_tree_array(tree.l, cX, options)
    @return_on_false complete cumulator
    @return_on_nonfinite_array cumulator T n
    (array2, complete2) = _eval_tree_array(tree.r, cX, options)
    @return_on_false complete2 cumulator
    @return_on_nonfinite_array array2 T n
    op = options.binops[op_idx]

    # We check inputs (and intermediates), not outputs.
    @inbounds @simd for j in 1:n
        x = op(cumulator[j], array2[j])::T
        cumulator[j] = x
    end
    # return (cumulator, finished_loop) #
    return (cumulator, true)
end

function deg1_eval(
    tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, options::Options
)::Tuple{AbstractVector{T},Bool} where {T<:Real,op_idx}
    n = size(cX, 2)
    (cumulator, complete) = _eval_tree_array(tree.l, cX, options)
    @return_on_false complete cumulator
    @return_on_nonfinite_array cumulator T n
    op = options.unaops[op_idx]
    @inbounds @simd for j in 1:n
        x = op(cumulator[j])::T
        cumulator[j] = x
    end
    return (cumulator, true) #
end

function deg0_eval(
    tree::Node, cX::AbstractMatrix{T}, options::Options
)::Tuple{AbstractVector{T},Bool} where {T<:Real}
    n = size(cX, 2)
    if tree.constant
        return (fill(convert(T, tree.val), n), true)
    else
        return (cX[tree.feature, :], true)
    end
end

function deg1_l2_ll0_lr0_eval(
    tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, ::Val{op_l_idx}, options::Options
)::Tuple{AbstractVector{T},Bool} where {T<:Real,op_idx,op_l_idx}
    n = size(cX, 2)
    op = options.unaops[op_idx]
    op_l = options.binops[op_l_idx]
    if tree.l.l.constant && tree.l.r.constant
        val_ll = convert(T, tree.l.l.val)
        val_lr = convert(T, tree.l.r.val)
        @return_on_check val_ll T n
        @return_on_check val_lr T n
        x_l = op_l(val_ll, val_lr)::T
        @return_on_check x_l T n
        x = op(x_l)::T
        @return_on_check x T n
        return (fill(x, n), true)
    elseif tree.l.l.constant
        val_ll = convert(T, tree.l.l.val)
        @return_on_check val_ll T n
        feature_lr = tree.l.r.feature
        cumulator = Array{T,1}(undef, n)
        @inbounds @simd for j in 1:n
            x_l = op_l(val_ll, cX[feature_lr, j])::T
            x = isfinite(x_l) ? op(x_l)::T : T(Inf) # These will get discovered by _eval_tree_array at end.
            cumulator[j] = x
        end
        return (cumulator, true)
    elseif tree.l.r.constant
        feature_ll = tree.l.l.feature
        val_lr = convert(T, tree.l.r.val)
        @return_on_check val_lr T n
        cumulator = Array{T,1}(undef, n)
        @inbounds @simd for j in 1:n
            x_l = op_l(cX[feature_ll, j], val_lr)::T
            x = isfinite(x_l) ? op(x_l)::T : T(Inf)
            cumulator[j] = x
        end
        return (cumulator, true)
    else
        feature_ll = tree.l.l.feature
        feature_lr = tree.l.r.feature
        cumulator = Array{T,1}(undef, n)
        @inbounds @simd for j in 1:n
            x_l = op_l(cX[feature_ll, j], cX[feature_lr, j])::T
            x = isfinite(x_l) ? op(x_l)::T : T(Inf)
            cumulator[j] = x
        end
        return (cumulator, true)
    end
end

# op(op2(x)) for x variable or constant
function deg1_l1_ll0_eval(
    tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, ::Val{op_l_idx}, options::Options
)::Tuple{AbstractVector{T},Bool} where {T<:Real,op_idx,op_l_idx}
    n = size(cX, 2)
    op = options.unaops[op_idx]
    op_l = options.unaops[op_l_idx]
    if tree.l.l.constant
        val_ll = convert(T, tree.l.l.val)
        @return_on_check val_ll T n
        x_l = op_l(val_ll)::T
        @return_on_check x_l T n
        x = op(x_l)::T
        @return_on_check x T n
        return (fill(x, n), true)
    else
        feature_ll = tree.l.l.feature
        cumulator = Array{T,1}(undef, n)
        @inbounds @simd for j in 1:n
            x_l = op_l(cX[feature_ll, j])::T
            x = isfinite(x_l) ? op(x_l)::T : T(Inf)
            cumulator[j] = x
        end
        return (cumulator, true)
    end
end

function deg2_l0_r0_eval(
    tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, options::Options
)::Tuple{AbstractVector{T},Bool} where {T<:Real,op_idx}
    n = size(cX, 2)
    op = options.binops[op_idx]
    if tree.l.constant && tree.r.constant
        val_l = convert(T, tree.l.val)
        @return_on_check val_l T n
        val_r = convert(T, tree.r.val)
        @return_on_check val_r T n
        x = op(val_l, val_r)::T
        @return_on_check x T n
        return (fill(x, n), true)
    elseif tree.l.constant
        cumulator = Array{T,1}(undef, n)
        val_l = convert(T, tree.l.val)
        @return_on_check val_l T n
        feature_r = tree.r.feature
        @inbounds @simd for j in 1:n
            x = op(val_l, cX[feature_r, j])::T
            cumulator[j] = x
        end
    elseif tree.r.constant
        cumulator = Array{T,1}(undef, n)
        feature_l = tree.l.feature
        val_r = convert(T, tree.r.val)
        @return_on_check val_r T n
        @inbounds @simd for j in 1:n
            x = op(cX[feature_l, j], val_r)::T
            cumulator[j] = x
        end
    else
        cumulator = Array{T,1}(undef, n)
        feature_l = tree.l.feature
        feature_r = tree.r.feature
        @inbounds @simd for j in 1:n
            x = op(cX[feature_l, j], cX[feature_r, j])::T
            cumulator[j] = x
        end
    end
    return (cumulator, true)
end

function deg2_l0_eval(
    tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, options::Options
)::Tuple{AbstractVector{T},Bool} where {T<:Real,op_idx}
    n = size(cX, 2)
    (cumulator, complete) = _eval_tree_array(tree.r, cX, options)
    @return_on_false complete cumulator
    @return_on_nonfinite_array cumulator T n
    op = options.binops[op_idx]
    if tree.l.constant
        val = convert(T, tree.l.val)
        @return_on_check val T n
        @inbounds @simd for j in 1:n
            x = op(val, cumulator[j])::T
            cumulator[j] = x
        end
    else
        feature = tree.l.feature
        @inbounds @simd for j in 1:n
            x = op(cX[feature, j], cumulator[j])::T
            cumulator[j] = x
        end
    end
    return (cumulator, true)
end

function deg2_r0_eval(
    tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, options::Options
)::Tuple{AbstractVector{T},Bool} where {T<:Real,op_idx}
    n = size(cX, 2)
    (cumulator, complete) = _eval_tree_array(tree.l, cX, options)
    @return_on_false complete cumulator
    @return_on_nonfinite_array cumulator T n
    op = options.binops[op_idx]
    if tree.r.constant
        val = convert(T, tree.r.val)
        @return_on_check val T n
        @inbounds @simd for j in 1:n
            x = op(cumulator[j], val)::T
            cumulator[j] = x
        end
    else
        feature = tree.r.feature
        @inbounds @simd for j in 1:n
            x = op(cumulator[j], cX[feature, j])::T
            cumulator[j] = x
        end
    end
    return (cumulator, true)
end

# Evaluate an equation over an array of datapoints
# This one is just for reference. The fused one should be faster.
function differentiable_eval_tree_array(
    tree::Node, cX::AbstractMatrix{T}, options::Options
)::Tuple{AbstractVector{T},Bool} where {T<:Real}
    n = size(cX, 2)
    if tree.degree == 0
        if tree.constant
            return (ones(T, n) .* tree.val, true)
        else
            return (cX[tree.feature, :], true)
        end
    elseif tree.degree == 1
        return deg1_diff_eval(tree, cX, Val(tree.op), options)
    else
        return deg2_diff_eval(tree, cX, Val(tree.op), options)
    end
end

function deg1_diff_eval(
    tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, options::Options
)::Tuple{AbstractVector{T},Bool} where {T<:Real,op_idx}
    (left, complete) = differentiable_eval_tree_array(tree.l, cX, options)
    @return_on_false complete left
    op = options.unaops[op_idx]
    out = op.(left)
    no_nans = !any(x -> (!isfinite(x)), out)
    return (out, no_nans)
end

function deg2_diff_eval(
    tree::Node, cX::AbstractMatrix{T}, ::Val{op_idx}, options::Options
)::Tuple{AbstractVector{T},Bool} where {T<:Real,op_idx}
    (left, complete) = differentiable_eval_tree_array(tree.l, cX, options)
    @return_on_false complete left
    (right, complete2) = differentiable_eval_tree_array(tree.r, cX, options)
    @return_on_false complete2 left
    op = options.binops[op_idx]
    out = op.(left, right)
    no_nans = !any(x -> (!isfinite(x)), out)
    return (out, no_nans)
end

end


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

JLLWrappers.@generate_main_file_header("Libcerf")
JLLWrappers.@generate_main_file("Libcerf", UUID("af83a40a-c4c4-57a0-81df-2309fbd279e3"))
end  # module Libcerf_jll


text:
using ApplicationBuilder
using CircuitscapeUI

dir(mod) = normpath(joinpath(dirname(pathof(mod)),".."))
const BUILD_FILE = joinpath(dir(CircuitscapeUI), "build", "Circuitscape.jl")
const SNOOP_FILE = joinpath(dir(CircuitscapeUI), "build", "snoop.jl")

function build_cs_binary(; build_path=pwd(), snoop = false, verbose = false)

    resources = [
                 joinpath(dir(CircuitscapeUI.Blink), "src", "AtomShell", "main.js"),
                 joinpath(dir(CircuitscapeUI.Blink), "src", "content", "main.html"),
                 joinpath(dir(CircuitscapeUI.Blink), "res"),
                 joinpath(dir(CircuitscapeUI.Blink), "res", "webio_setup.js"),
                 joinpath(dir(CircuitscapeUI.Tachyons), "assets", "tachyons.min.css"),
                 joinpath(dir(CircuitscapeUI.WebIO), "assets", "webio", "dist", "bundle.js"),
                 joinpath(dir(CircuitscapeUI.InteractBase), "assets", "all.js"),
                 joinpath(dir(CircuitscapeUI.InteractBase), "assets", "style.css"),
                 joinpath(dir(CircuitscapeUI.InteractBulma), "assets", "main.css"),
                 joinpath(dir(CircuitscapeUI.InteractBulma), "assets", "bulma.min.css"),
                 joinpath(dir(CircuitscapeUI.InteractBulma), "assets", "bulma-slider.min.css"),
                 joinpath(dir(CircuitscapeUI.InteractBulma), "assets", "bulma-switch.min.css"),
                 joinpath(dir(CircuitscapeUI.InteractBulma), "assets", "bulma-checkradio.min.css"),
                 joinpath(dir(CircuitscapeUI.InteractBulma), "assets", "bulma-tooltip.min.css"),
                 joinpath(dir(CircuitscapeUI.InteractBulma), "assets", "bulma-accordion.min.css"),
                 joinpath(dir(CircuitscapeUI.Knockout), "assets", "knockout.js"),
                 joinpath(dir(CircuitscapeUI.Knockout), "assets", "knockout_punches.js"),
                 joinpath(dir(CircuitscapeUI.Circuitscape), "test"),
                 joinpath(dir(CircuitscapeUI), "assets", "cs_logo.ico")
                ]

    libs = readdir(joinpath(dir(CircuitscapeUI.MbedTLS), "deps", "usr", "lib"))
    lib_paths = joinpath.(dir(CircuitscapeUI.MbedTLS), "deps", "usr", "lib", libs)

    electron_app = joinpath(dir(CircuitscapeUI.Blink), "deps/Julia.app")
    if !ispath(electron_app)
        CircuitscapeUI.Blink.AtomShell.install()
    end
    libraries = vcat(
                     electron_app,
                     lib_paths
                    )

    julia_path = joinpath(split(string(Base.julia_cmd()))[1][2:end])

    if Sys.isapple()
        system_path = joinpath("Circuitscape.app", "Contents", "MacOS")
        bin_path = joinpath(build_path, system_path)
        mkpath(bin_path, mode=0o755)
        run(`cp $julia_path $bin_path`)
        run(Cmd(`install_name_tool -rpath '@executable_path/../lib' '@executable_path' julia`, dir=bin_path))
    end

    if snoop
        build_app_bundle(BUILD_FILE,
            resources = resources,
            builddir = build_path,
            snoopfile = SNOOP_FILE,
            libraries = libraries,
            icns_file = joinpath(dir(CircuitscapeUI), "assets", "circuitscape.icns"),
            verbose = verbose)
    else
        build_app_bundle(BUILD_FILE,
            resources = resources,
            builddir = build_path,
            libraries = libraries,
            icns_file = joinpath(dir(CircuitscapeUI), "assets", "circuitscape.icns"),
            verbose = verbose)
    end

end


text:
import YAML

const WRONG_FORMAT_CONFIG_FILE = joinpath(dirname(pathof(PowerSystems)),
                                            "descriptors", "config.yml")
include(joinpath(DATA_DIR,"data_5bus_pu.jl"))

@testset "Test reading in config data" begin
    data = IS.read_validation_descriptor(PSY.POWER_SYSTEM_STRUCT_DESCRIPTOR_FILE)
    @test data isa Vector
    @test !isempty(data)
    function find_struct()
        for item in data
            if item["struct_name"] == "TechThermal"
                return true
            end
        end
        return false
    end
    @test find_struct()
    @test_throws(ErrorException, IS.read_validation_descriptor("badfile.toml"))
end

@testset "Test adding custom validation YAML file to System" begin
    sys_no_config = System(nodes5, thermal_generators5, loads5, nothing, nothing,
                            100.0, nothing, nothing, nothing; runchecks=true)
    @test !isempty(sys_no_config.data.validation_descriptors)

    sys_no_runchecks = System(nodes5, thermal_generators5, loads5, nothing, nothing,
                            100.0, nothing, nothing, nothing; runchecks=false)
    @test isempty(sys_no_runchecks.data.validation_descriptors)
end

@testset "Test extracting struct info from validation_descriptor vector" begin
    data =  [Dict("fields"=>Dict{Any,Any}[
            Dict("name"=>"curtailpenalty","valid_range"=>Dict{Any,Any}("max"=>nothing,"min"=>0.0)),
            Dict("name"=>"variablecost","valid_range"=>Dict{Any,Any}("max"=>nothing,"min"=>0.0)),
            Dict("name"=>"internal")],
                "struct_name"=>"EconHydro"),
            Dict("fields"=>Dict{Any,Any}[
            Dict("name"=>"curtailpenalty","valid_range"=>Dict{Any,Any}("max"=>nothing,"min"=>0.0)),
            Dict("name"=>"variablecost","valid_range"=>Dict{Any,Any}("max"=>nothing,"min"=>0.0)),
            Dict("name"=>"internal")],
                "struct_name"=>"EconLoad")]
    struct_name = "EconHydro"
    descriptor = IS.get_config_descriptor(data, struct_name)
    @test descriptor isa Dict{String,Any}
    @test haskey(descriptor, "struct_name")
    @test haskey(descriptor, "fields")
    @test descriptor["struct_name"] == struct_name

end

@testset "Test extracting field info from struct descriptor dictionary" begin
    config = Dict{Any,Any}("fields"=>Dict{Any,Any}[
                            Dict("name"=>"name","data_type"=>"String"),
                            Dict("name"=>"available","data_type"=>"Bool"),
                            Dict("name"=>"bus","data_type"=>"Bus"),
                            Dict("name"=>"tech","data_type"=>"Union{Nothing, TechThermal}"),
                            Dict("name"=>"econ","data_type"=>"Union{Nothing, EconThermal}"),
                            Dict("name"=>"internal","data_type"=>"IS.InfrastructureSystemsInternal")],
                            "struct_name"=>"ThermalStandard")
    field_name = "econ"
    field_descriptor = IS.get_field_descriptor(config, field_name)
    @test field_descriptor isa Dict{Any, Any}
    @test haskey(field_descriptor, "name")
    @test field_descriptor["name"] == field_name
end

@testset "Test retrieving validation action" begin
    warn_descriptor = Dict{Any,Any}("name"=>"ramplimits",
                        "valid_range"=>Dict{Any,Any}("max"=>5,"min"=>0),
                        "validation_action"=>"warn")
    error_descriptor = Dict{Any,Any}("name"=>"ramplimits",
                        "valid_range"=>Dict{Any,Any}("max"=>5,"min"=>0),
                        "validation_action"=>"error")
    typo_descriptor = Dict{Any,Any}("name"=>"ramplimits",
                        "valid_range"=>Dict{Any,Any}("max"=>5,"min"=>0),
                        "validation_action"=>"asdfasdfsd")
    @test IS.get_validation_action(warn_descriptor) == IS.validation_warning
    @test IS.get_validation_action(error_descriptor) == IS.validation_error
    @test_throws(ErrorException, IS.get_validation_action(typo_descriptor))
end

@testset "Test field validation" begin
    #test recursive call of validate_fields and a regular valid range
    bad_therm_gen_rating = deepcopy(thermal_generators5)
    bad_therm_gen_rating[1].tech.rating = -10
    @test_logs((:error, r"Invalid range"),
        @test_throws(PSY.InvalidRange,
                        System(nodes5, bad_therm_gen_rating, loads5, nothing, nothing,
                               100.0, nothing, nothing, nothing; runchecks=true)
        )
    )

    #test custom range (activepower and activepowerlimits)
    bad_therm_gen_act_power = deepcopy(thermal_generators5)
    bad_therm_gen_act_power[1].activepower = 10
    @test_logs (:warn, r"Invalid range") System(nodes5, bad_therm_gen_act_power, loads5,
            nothing, nothing, 100.0, nothing, nothing, nothing; runchecks=true)

    #test validating named tuple
    bad_therm_gen_ramp_lim = deepcopy(thermal_generators5)
    bad_therm_gen_ramp_lim[1].tech.ramplimits = (up = -10, down = -3)
    @test_logs((:error, r"Invalid range"), match_mode=:any,
        @test_throws(PSY.InvalidRange,
                     System(nodes5, bad_therm_gen_ramp_lim, loads5, nothing, nothing, 100.0,
                            nothing, nothing, nothing; runchecks=true)
        )
    )
end

@testset "Test field validation" begin
    sys = System(nodes5, thermal_generators5, loads5, nothing, nothing,
                               100.0, nothing, nothing, nothing; runchecks=true)

    add_component!(sys,Bus(11,"11",PSY.PQ,1,1,(min=.9,max=1.1),123))
    B = get_components(Bus,sys) |> collect
    a = Arc(B[1],B[6])
    badline = Line("badline",true,0.01,0.01,a,0.002,0.014,(from = 0.015, to = 0.015),5.0,(min = -1, max = 1))
    @test_logs((:error, r"cannot create Line"), match_mode=:any,
        @test_throws(PSY.InvalidValue,
                     add_component!(sys, badline)
        )
    )
end

@testset "Test field validation after deserialization" begin
    sys = System(nodes5, thermal_generators5, loads5, nothing, nothing,
                 100.0, nothing, nothing, nothing)

    add_component!(sys, Bus(11, "11", PSY.PQ, 1, 1, (min=.9, max=1.1), 123))
    path, io = mktemp()
    PSY.to_json(sys, path)
    try
        to_json(io, sys)
    catch
        close(io)
        rm(path)
        rethrow()
    end
    close(io)

    try
        sys2 = PSY.System(path)

        B = get_components(Bus, sys2) |> collect
        a = Arc(B[1], B[6])
        badline = Line("badline", true, 0.01, 0.01, a, 0.002, 0.014,
                       (from = 0.015, to = 0.015), 5.0, (min = -1, max = 1))
        @test_logs((:error, r"cannot create Line"), match_mode=:any,
            @test_throws(PSY.InvalidValue,
                         add_component!(sys2, badline)
            )
        )
    finally
        rm(path)
    end
end


text:
import StateSpaceReconstruction:
	minima_and_stepsizes,
    splitaxes,
    Rectangle

using RecipesBase

@recipe function f(r::RectangularInvariantMeasure,
    boxfillfactor::Int = 3, linesegments = true,
    lw = 0.8, lc = :black, lα = 0.5, ls = :do,
    ms = 0.2, mc = :black, mα = 0.3)

    # All bins are assigned a bin origin.
    unique_visited = unique(r.visited_bins_coordinates, dims = 2)

    # Bins may be visited by several times, so we'll get rid
    # of the repetitions.
    n_visited_boxes = size(unique_visited, 2)

    D = size(unique_visited, 1)

    if D != 3
        throw(ArgumentError("Plot recipe not implemented for dimension $D."))
    end

    # Plot the partition grid over the points of the reconstructed
    # orbit.
    legend --> false

    for i = 1:n_visited_boxes
        origin = unique_visited[:, i]
        rect = Rectangle(
            SVector{3, Float64}(unique_visited[:, i]),
            SVector{3, Float64}(r.ϵ_absolute))

        @series begin
            rect
        end

        # Get the measure of the box, and fill the box with a
        # number of points scaled to the measure.
        #μ = riv.measure.dist[i]
        #fillpts = fill_hyperrectangle_3D(origin,
        #            stepsizes,
        #            ceil(Int, μ*100)*boxfillfactor)
    end

    @series begin
        seriestype := :scatter3d
        splitaxes(r.points)
    end

end

@recipe function f(r::InducedRectangularInvariantMeasure, plot_boxes = true)

    legend --> false

    # Plot all boxes at the target resolution with nonzero measure.
    bins_with_nonzero_measure_ϵF = r.bins_ϵF[:, findall(r.measure_induced.dist .> 0)]
    n_visited_boxes_ϵF = maximum(size(bins_with_nonzero_measure_ϵF))

    for i = 1:n_visited_boxes_ϵF
        origin = r.bins_ϵF[:, i]
        rect = Rectangle(
            SVector{3, Float64}(bins_with_nonzero_measure_ϵF[:, i]),
            SVector{3, Float64}(r.ϵF_absolute))

        @series begin
            lc --> :black
            lw --> 2
            rect
        end
    end

    # Plot all boxes visited at the source resolution.
    n_visited_boxes_ϵj = maximum(size(r.bins_visited_ϵj))

    for i = 1:n_visited_boxes_ϵj
        origin = r.bins_visited_ϵj[:, i]
        rect = Rectangle(
            SVector{3, Float64}(r.bins_visited_ϵj[:, i]),
            SVector{3, Float64}(r.ϵj_absolute))

        @series begin
            lc --> :green
            rect
        end
    end

    @series begin
        seriestype := :scatter3d
        mc --> :red
        splitaxes(r.points)
    end
end

@recipe function f(r::AverageRectangularInvariantMeasure, boxes = true, measure = false)

    if boxes && !measure
        #####################
        # TARGET RESOLUTION
        #####################

        # All bins are assigned a bin origin.
        unique_visited = unique(r.measure_ϵF.visited_bins_coordinates, dims = 2)

        # Bins may be visited by several times, so we'll get rid
        # of the repetitions.
        n_visited_boxes = size(unique_visited, 2)
          D = size(unique_visited, 1)

        if D != 3
            throw(ArgumentError("Plot recipe not implemented for dimension $D."))
        end

        # Plot the partition grid over the points of the reconstructed
        # orbit.
        legend --> false

        for i = 1:n_visited_boxes
            origin = unique_visited[:, i]
            rect = Rectangle(
                SVector{3, Float64}(unique_visited[:, i]),
                SVector{3, Float64}(r.measure_ϵF.ϵ_absolute))

            @series begin
                lw --> 2
                rect
            end
        end

        ######################################################
        # RESOLUTIONS FROM WHICH THE MEASURE IS INDUCED
        ######################################################

        cols = [:red, :purple, :green, :blue, :yellow, :magenta, :pink]
        for (k, measure) in enumerate(r.induced_measures)
            # Plot all boxes visited at the source resolution.
            n_visited_boxes_ϵj = maximum(size(measure.bins_visited_ϵj))

            for i = 1:n_visited_boxes_ϵj
                origin = measure.bins_visited_ϵj[:, i]
                rect = Rectangle(
                    SVector{3, Float64}(measure.bins_visited_ϵj[:, i]),
                    SVector{3, Float64}(measure.ϵj_absolute))

                @series begin
                    lc --> cols[k]
                    rect
                end
            end
        end

        ######################################################
        # THE POINTS FROM WHICH THE MEASURE IS INDUCED
        ######################################################
        @series begin
            seriestype := :scatter3d
            splitaxes(r.measure_ϵF.points)
        end

    end
end


text:
function SemiVariance()
    return SemiVariance(())
end

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

function SemiVariance(arg0::SemiVariance_Direction)
    return SemiVariance((SemiVariance_Direction,), arg0)
end

function SemiVariance(arg0::jboolean)
    return SemiVariance((jboolean,), arg0)
end

function SemiVariance(arg0::jboolean, arg1::SemiVariance_Direction)
    return SemiVariance((jboolean, SemiVariance_Direction), arg0, arg1)
end

function copy(obj::SemiVariance)
    return jcall(obj, "copy", SemiVariance, ())
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function evaluate(obj::AbstractUnivariateStatistic)
    return jcall(obj, "evaluate", jdouble, ())
end

function evaluate(obj::SemiVariance, arg0::Vector{jdouble}, arg1::SemiVariance_Direction)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, SemiVariance_Direction), arg0, arg1)
end

function evaluate(obj::SemiVariance, arg0::Vector{jdouble}, arg1::jdouble)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jdouble), arg0, arg1)
end

function evaluate(obj::SemiVariance, arg0::Vector{jdouble}, arg1::jdouble, arg2::SemiVariance_Direction)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jdouble, SemiVariance_Direction), arg0, arg1, arg2)
end

function evaluate(obj::SemiVariance, arg0::Vector{jdouble}, arg1::jdouble, arg2::SemiVariance_Direction, arg3::jboolean, arg4::jint, arg5::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jdouble, SemiVariance_Direction, jboolean, jint, jint), arg0, arg1, arg2, arg3, arg4, arg5)
end

function evaluate(obj::SemiVariance, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function get_data(obj::AbstractUnivariateStatistic)
    return jcall(obj, "getData", Vector{jdouble}, ())
end

function get_variance_direction(obj::SemiVariance)
    return jcall(obj, "getVarianceDirection", SemiVariance_Direction, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function is_bias_corrected(obj::SemiVariance)
    return jcall(obj, "isBiasCorrected", jboolean, ())
end

function set_data(obj::AbstractUnivariateStatistic, arg0::Vector{jdouble})
    return jcall(obj, "setData", void, (Vector{jdouble},), arg0)
end

function set_data(obj::AbstractUnivariateStatistic, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(obj, "setData", void, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function with_bias_corrected(obj::SemiVariance, arg0::jboolean)
    return jcall(obj, "withBiasCorrected", SemiVariance, (jboolean,), arg0)
end

function with_variance_direction(obj::SemiVariance, arg0::SemiVariance_Direction)
    return jcall(obj, "withVarianceDirection", SemiVariance, (SemiVariance_Direction,), arg0)
end


text:
#
# This file is part of the YAActL.jl Julia package, MIT license
#
# Paul Bayer, 2020
#

function __init__()
    if myid() == 1
        global USR = Link(
            RemoteChannel(()->Channel{Message}(32)),
            myid(), :remote)
        global _REG = Link(
            RemoteChannel(()->Actor(_reg, Dict{Symbol, Link}()).chn),
            1, :remote)
    else
        tmp = Actor(1, ()->YAActL.USR)
        global USR = call!(tmp)
        become!(tmp, ()->YAActL._REG)
        global _REG = call!(tmp)
        exit!(tmp)
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, SuiteSparse_jll

@testset "SuiteSparse_jll" begin
    @test ccall((:SuiteSparse_version, libsuitesparseconfig), Cint, (Ptr{Cint},), C_NULL) == 5008
end

max_stars_repo_name: jmkuhn/julia

max_issues_repo_name: simsurace/julia

max_forks_repo_name: simsurace/julia


text:
function dfs_path(netw::SimpleWeightedGraph, startnode::Int64)
    n_vertex = nv(netw)
    path = zeros(Int, n_vertex)
    visited = zeros(Int, n_vertex)
    parents = zeros(Int, n_vertex)

    function _dfs(at)
        if (visited[at] == 1)
            nothing
        else
            visited[at] = 1
            if parents[at] > 0
                pa = parents[at]
            end
            path[at] = pa
            neig = neighbors(netw, at)
            for p in neig
                parents[p] = at
                _dfs(p)
            end
        end
    end

    at = startnode
    pa = startnode
    parents[at] = startnode
    _dfs(at)

    return path
end

function path_to_adj_matrix(path::Vector{Int64})
    dim = maximum([maximum(path), length(path)])
    adj_mat = zeros(dim, dim)

    for i in 1:length(path)
        if (path[i] != 0) && (path[i] != i)
            adj_mat[path[i], i] = 1
        end
    end

    return SparseMatrixCSC(adj_mat)
end

function node_color(netw::AbstractSimpleWeightedGraph)
    adj = netw.weights
    count = 0
    n = nv(netw)
    components = zeros(Int, n)
    visited = zeros(Int, n)

    function findComponents()
        for i = 1:n
            if (visited[i] == 0)
                count += 1
                dfs(i)
            end
        end
        return count, components
    end

    function dfs(at)
        visited[at] = 1
        components[at] = count
        neigh = neighbors(netw, at)
        for next in neigh
            if (visited[next] == 0)
                dfs(next)
            end
        end
    end

    count, components = findComponents()
    @show count
    return components

end

function bfs_path(g::AbstractSimpleWeightedGraph, startnode::Int64, endnode::Int64)
    n = nv(g)

    q = Queue{Int64}()
    enqueue!(q, startnode)

    visitied = zeros(Int64, 1, n)
    visitied[startnode] = 1

    prev = zeros(Int64, 1, n)

    while (isempty(q.store) == false)
        node = dequeue!(q)
        neigh = neighbors(g, node)
        for next in neigh

            if (visitied[next] != 1)
                enqueue!(q, next)
                visitied[next] = 1
                prev[next] = node
            end
        end
    end

    if prev[endnode] == 0
        return []
    end

    at = endnode
    path = []
    while true
        at = prev[at]
        if at == 0
            break
        else
            push!(path, at)
        end
    end

    path = reverse(path)
    push!(path, endnode)
    return path

end

# shortest_path using dijkstra_shortest_paths form `Graphs.jl`:
function shortest_path(g::AbstractSimpleWeightedGraph, source::Int64, sink::Int64)
    ds = Graphs.dijkstra_shortest_paths(g, sink)
    try
        return reverse(spath(ds, source, sink))
    catch e
        return []
    end

end

# convert precedessor list in path:
spath(ds, source, sink) = source == sink ? source : [spath(ds, ds.parents[source], sink) source]

max_stars_repo_name: tragisch/NetworkOptimization

max_issues_repo_name: tragisch/NetworkOptimization

max_forks_repo_name: tragisch/NetworkOptimization


text:
using IRTools
using IRTools: IR, Variable, Statement, Lambda, @dynamo, functional, cond,
  argument!, arguments, isexpr, stmt, xcall, return!, returnvalue

struct Func
  f # Avoid over-specialising on the continuation object.
end

(f::Func)(args...) = f.f(args...)

function captures(ir, vs)
  us = Set()
  for v in vs
    isexpr(ir[v].expr) || continue
    foreach(x -> x isa Variable && push!(us, x), ir[v].expr.args)
  end
  push!(us, returnvalue(IRTools.block(ir, 1)))
  return setdiff(us, vs)
end

rename(env, x) = x
rename(env, x::Variable) = env[x]
rename(env, x::Expr) = Expr(x.head, rename.((env,), x.args)...)
rename(env, x::Statement) = stmt(x, expr = rename(env, x.expr))

excluded = [GlobalRef(Base, :getindex)]

function continuation!(bl, ir, env, vs, ret)
  rename(x) = Effects.rename(env, x)
  local v, st
  while true
    isempty(vs) && return return!(bl, rename(Expr(:call, ret, returnvalue(IRTools.block(ir, 1)))))
    v = popfirst!(vs)
    st = ir[v]
    isexpr(st.expr, :call) && !(st.expr.args[1] ∈ excluded) && break
    isexpr(st.expr, :lambda) &&
      (st = stmt(st, expr = Expr(:lambda, cpslambda(st.expr.args[1]), st.expr.args[2:end]...)))
    env[v] = push!(bl, rename(st))
  end
  cs = [ret, setdiff(captures(ir, vs), [v])...]
  next = push!(bl, Expr(:lambda, continuation(ir, vs, cs, v, ret), rename.(cs)...))
  next = xcall(Effects, :Func, next)
  ret = push!(bl, stmt(st, expr = xcall(Effects, :cps, next, rename(st.expr).args...)))
  return!(bl, ret)
end

function continuation(ir, vs, cs, in, ret)
  bl = empty(ir)
  env = Dict()
  self = argument!(bl)
  env[in] = argument!(bl)
  for (i, c) in enumerate(cs)
    env[c] = pushfirst!(bl, xcall(:getindex, self, i))
  end
  continuation!(bl, ir, env, vs, ret)
end

cpslambda(ir) = cpstransform(ir, true)

function cpstransform(ir, lambda = false)
  lambda || (ir = functional(ir))
  k = argument!(ir, at = lambda ? 2 : 1)
  bl = empty(ir)
  env = Dict()
  for arg in arguments(ir)
    env[arg] = argument!(bl)
  end
  continuation!(bl, ir, env, keys(ir), k)
end

cps_lambda(f::Lambda{S,I}, k) where {S,I} =
  Lambda{S,I}(x -> cps(k, first(f.data), x), Base.tail(f.data)...)

cps(k, f::IRTools.Lambda{<:Tuple{typeof(cps),Vararg{Any}}}, args...) = f(k, args...)
cps(k, f::Func, args...) = cps_lambda(f.f, k)(args...)

cps(k, f::Core.IntrinsicFunction, args...) = k(f(args...))
cps(k, ::typeof(Core._apply), f, args...) = Core._apply(cps, (k, f), args...)
cps(k, ::typeof(cond), c, t, f) = c ? cps(k, t) : cps(k, f)
cps(k′, ::typeof(cps), k, args...) = cps(x -> cps(k′, k, x), args...)

for f in :[getfield, typeof, Core.apply_type, typeassert, (===), ifelse,
           Core.sizeof, Core.arrayset, tuple, isdefined, fieldtype, nfields,
           isa, Core.arraysize, repr, print, println, Base.vect, Broadcast.broadcasted,
           Broadcast.materialize, Core.Compiler.return_type, Base.union!, Base.getindex, Base.haskey,
           Base.pop!, Base.setdiff].args
  @eval cps(k, ::typeof($f), args...) = k($f(args...))
end

@dynamo function cps(k, f, args...)
  ir = IR(f, args...)
  ir == nothing && error("No IR for $((f, args...))")
  cpstransform(ir)
end


text:
using CUDA

@inline function CheckBlockBounds(offset, block, grids)
        x::Int32 = offset.x + block.x
        if x < 1 || x > grids
                return true
        end
        x = offset.y + block.y
        if x < 1 || x > grids
                return true
        end
        x = offset.z + block.z
        if x < 1 || x > grids
                return true
        end
        return false
end

@inline function CalcDistance(query, point)
    tempdist::Float32 = 0
    for d = 1:length(query)
        @inbounds tempdist += CUDA.pow(query[d] - point[d], 2)
    end
    tempdist = CUDA.sqrt(tempdist)
    return tempdist
end

function cuGridKnnSimpleViewFunction(devPoints, devQueries,
        Pointsperblock, Queriesperblock,
        IntPointsperblock, IntQueriesperblock,
        devDistances, devNeighbours,
        dimensions, grids, offset)

    if CheckBlockBounds(offset, blockIdx(), grids)
        return nothing
    end
    q_bid::Int32 = ((blockIdx().z - 1)*gridDim().y*gridDim().x
            + (blockIdx().y - 1)*gridDim().x
            + blockIdx().x)
    p_bid::Int32 = ((blockIdx().z + offset.z- 1)*gridDim().y*gridDim().x
            + (blockIdx().y + offset.y - 1)*gridDim().x
            + blockIdx().x + offset.x)
    tid::Int32 = threadIdx().x + (threadIdx().y-1)*blockDim().x
    stride::Int32 = (blockDim().x)*(blockDim().y)

    SharedPoints = @cuDynamicSharedMem(Float32, (dimensions, stride))
    SharedQueries = @cuDynamicSharedMem(Float32, (dimensions, stride), dimensions*stride*sizeof(Float32))

    @inbounds startPoints::Int32 = IntPointsperblock[p_bid]
    @inbounds startQueries::Int32 = IntQueriesperblock[q_bid]
    @inbounds totalPoints::Int32 = Pointsperblock[p_bid]
    @inbounds totalQueries::Int32 = Queriesperblock[q_bid]
    @inbounds Points = @view devPoints[(startPoints+1):(startPoints+totalPoints), :]
    @inbounds Queries = @view devQueries[(startQueries+1):(startQueries+totalQueries), :]
    @inbounds Distances = @view devDistances[(startQueries+1):(startQueries+totalQueries)]
    @inbounds Neighbours = @view devNeighbours[(startQueries+1):(startQueries+totalQueries)]
    @inbounds query = @view SharedQueries[:, tid]
    dist::Float32 = 0
    nb::Int32 = 0
    for q::Int32 in 0:stride:(totalQueries-1)
        if tid + q <= totalQueries
            for d = 1:dimensions
                @inbounds  query[d] = Queries[tid+q, d]
            end
            @inbounds dist = Distances[tid+q]
            @inbounds nb = Neighbours[tid+q]
        end
        for p::Int32 in 0:stride:(totalPoints-1)
            sync_threads()
            if tid + p <= totalPoints
                for d::Int32 in 1:dimensions
                    @inbounds SharedPoints[d, tid] = Points[tid+p, d]
                end
            end
            sync_threads()
            bounds::Int32 = CUDA.min(stride, totalPoints-p)
            for i::Int32 in 1:bounds
                @inbounds point = @view SharedPoints[:, (i+tid-2)%bounds+1]
                tempdist::Float32 = CalcDistance(query, point)
                if tempdist < dist
                    dist = tempdist
                    nb = startPoints + p + (i+tid-2)%bounds+1
                end
            end
        end
        if tid + q <= totalQueries
            @inbounds Distances[tid+q] = dist
            @inbounds Neighbours[tid+q] = nb
        end
    end
    return nothing
end

function cuda_knn_simple_view_function(OrderedPoints, OrderedQueries,
    PointsPerBlock, QueriesPerBlock,
    IntegralPointsPerBlock, IntegralQueriesPerBlock,
    numOfPoints, numOfQueries, numOfGrids, dimensions)

    ## Data transfer
    devPoints = CuArray(OrderedPoints)
    devQueries = CuArray(OrderedQueries)
    devPointsPerBlock = CuArray(PointsPerBlock)
    devQueriesPerBlock = CuArray(QueriesPerBlock)
    devIntegralPointsPerBlock = CuArray(IntegralPointsPerBlock)
    devIntegralQueriesPerBlock = CuArray(IntegralQueriesPerBlock)
    devRes = CUDA.fill(Float32(100), numOfQueries)
    devNeighbours = CuArray(zeros(Int32, numOfQueries))
    ## Config
    dimensions = Int32(dimensions)
    grids = Int32(numOfGrids)
    offset = (x=Int32(0),y=Int32(0),z=Int32(0))
    thread_groups = 2
    shmem=2*thread_groups*32*dimensions*sizeof(Float32)
    ## Running for main block.
    @cuda(blocks=(grids,grids,grids), threads=(32,thread_groups), shmem=shmem,
      cuGridKnnSimpleViewFunction(devPoints, devQueries,
        devPointsPerBlock, devQueriesPerBlock,
        devIntegralPointsPerBlock, devIntegralQueriesPerBlock,
        devRes, devNeighbours, dimensions, grids, offset))
    #Running for neighbour blocks.
    for x::Int32 = -1:1
        for y::Int32 = -1:1
            for z::Int32 = -1:1
                if x == 0 && y == 0 && z == 0
                    continue
                end
                offset = (x=x,y=y,z=z)
                @cuda(blocks=(grids,grids,grids), threads=(32,thread_groups), shmem=shmem,
                  cuGridKnnSimpleViewFunction(devPoints, devQueries,
                    devPointsPerBlock, devQueriesPerBlock,
                    devIntegralPointsPerBlock, devIntegralQueriesPerBlock,
                    devRes, devNeighbours, dimensions, grids, offset))
            end
        end
    end

    return Array(devNeighbours), Array(devRes)

end


text:
#######################################################################################
#                                   Load Packages
#######################################################################################
cd(@__DIR__)
using Pkg
Pkg.activate("../..")
using StatsPlots, ACTRModels, Distributions, Turing
include("Semantic_Model.jl")
Random.seed!(354301)
#######################################################################################
#                                   Generate Data
#######################################################################################
# true blc value
blc = 1.0
# fixed parameters
parms = (noise = true, τ = 0.0, s = .2, mmp = true, δ = 1.0)
stimuli = get_stimuli()
n_reps = 10
data = map(x -> simulate(parms, x, n_reps; blc), stimuli)
#######################################################################################
#                                    Define Model
#######################################################################################
@model model(data, parms) = begin
    blc ~ Normal(1.0, 1.0)
    data ~ Semantic(blc, parms)
end
#######################################################################################
#                                 Estimate Parameters
#######################################################################################
# Settings of the NUTS sampler.
n_samples = 1500
n_adapt = 1500
specs = NUTS(n_adapt, 0.65)
n_chains = 4
chain = sample(model(data, parms), specs, MCMCThreads(), n_samples, n_chains, progress=true)
#######################################################################################
#                                      Plot Chains
#######################################################################################
pyplot()
posteriors = plot(chain, seriestype=:density, grid=false, title="",
    size=(300,175), xaxis=font(8), yaxis=font(8))
plot!(xlabel="blc")
#######################################################################################
#                                  Posterior Predictive
#######################################################################################
font_size = 7
hit_rates(s) = posterior_predictive(x -> hit_rate(parms, s, n_reps; x...), chain, 1000)
preds = map(s -> hit_rates(s), stimuli)
predictive_plot = histogram(preds, xlabel="% Correct" ,ylabel="Probability", xaxis=font(font_size), yaxis=font(font_size),
    grid=false, color=:grey, leg=false, titlefont=font(font_size), xlims=(0,1.1),
    layout=(2,1), ylims=(0,.4), normalize=:probability, size=(300,250), title=["Is a canary a bird?" "Is a canary an animal?"])
savefig(posteriors, "Semantic_Posterior.eps")
savefig(predictive_plot, "Semantic_Postpred.eps")


text:
using GLWindow
using Base.Test

function is_ci()
	get(ENV, "TRAVIS", "") == "true" || get(ENV, "APPVEYOR", "") == "true" || get(ENV, "CI", "") == "true"
end

if !is_ci() # only do test if not CI... this is for automated testing environments which fail for OpenGL stuff, but I'd like to test if at least including works

GLFW.Init()
window = Screen("test", resolution=(500,500))

swapbuffers(window)

pollevents()

end
println("success")


text:
#=
Pkg.add("Plots")
Pkg.add("PyPlot")
Pkg.add("GR")
Pkg.add("UnicodePlots")
Pkg.add("PlotlyJS")

=#
using Plots

unicodeplots();
plot(eq_values)

plot!(
    b,

    label  = "equation of time (calculated)",
    #line=(:black, 0.5, 6, :solid),

    #size=(800, 600),

    #xticks = (1:14:366, datestrings[1:14:366]),
    #yticks = -20:2.5:20,

    ylabel = "Minutes faster or slower than GMT",
    xlabel = "day in year",

    title  = "The Equation of Time",
    xrotation = rad2deg(pi/3),

    #fillrange = 0,
    #fillalpha = 0.25,
    #fillcolor = :lightgoldenrod,

    #background_color = :ivory
    )

    for t in [144,180,216]#[36,72,144,180,216,288]
      nx.newGraph()
      normalise(x) = x/norm(x)

      links = filter(i -> flux[i][t]>0 , 1:len(flux))
      tflux = [log10(flux[i][t]) for i in links]
      weight = 1+normalise(tflux)

      c= 0
      for i in links
          nx.addedge(edges[i][1],edges[i][2],weight[c+=1])
          println(c)
      end

      nx.p()

      a = nx.eigenvector("none")
      b = [try a[i] catch 0.0f0 end for i in nodes]

      a1 = nx.eigenvector()
      b1 = [try a1[i] catch 0.0f0 end for i in nodes]

      difference = [abs(i) for i in (b-b1)]
      difference= [difference[i]/(b[i]+1e-9) for i in 1:len(b)]

    if test
      plot!(
          difference,

          label  = "$t",
          ylabel = "Eigenvector Centrality",
          xlabel = "Species Index",

          title  = "Unicode plots: Eigen-centrality Comparison",
          xrotation = rad2deg(pi/3),
          )

    else
        plot(
            difference,

            label  = "$t",
            ylabel = "Eigenvector Centrality",
            xlabel = "Species Index",

            title  = "Unicode plots: Eigen-centrality Comparison",
            xrotation = rad2deg(pi/3),
            )

    end

      test=true

    end

max_stars_repo_name: wolfiex/RopaJL

max_issues_repo_name: wolfiex/RopaJL

max_forks_repo_name: wolfiex/RopaJL


text:
OS_NAME == :Windows && Pkg.add("FactCheck")

using FactCheck, Compat, Base.Test, Docile

import Compat.String

include("helpers.jl")
VERSION < v"0.4-dev+6619" && include(joinpath("Legacy", "facts.jl"))
include(joinpath("Collector", "facts.jl"))
include(joinpath("Interface", "facts.jl"))
include(joinpath("Runner", "facts.jl"))
include(joinpath("Formats", "facts.jl"))
include(joinpath("Extensions", "facts.jl"))

isinteractive() || FactCheck.exitstatus()


text:
using SymPy

@syms d r_off lambda r_opp

# sipler variant
t_x = lambda*(d - r_opp*r_opp/d) # cos(π/2 + asin(r_opp/d)) = -r_opp/d
t_y = lambda*r_opp*sqrt(1 - (r_opp*r_opp)/(d*d)) #  sin(π/2 + asin(r_opp/d)) = sqrt(1 - (r_opp*r_opp)/(d*d))
toff_x = r_off*(r_off - r_opp)/d # cos(π/2 + asin((r_opp - r_off)/d)) = (r_off - r_opp)/d
toff_y = r_off*sqrt(1 - ((r_opp-r_off)/d)^2) # sin(π/2 + asin((r_opp - r_off)/d)) = sqrt(1 - ((r_opp-r_off)/d)^2)

lambda_sol = solve(t_x - toff_x,lambda)[1]
ty_ins = subs(t_y,lambda,lambda_sol)
l = ty_ins - toff_y
sol = solve(l, r_off)
sol[3]

simplify(subs(sol[3],r_opp,1))

using Plots, LaTeXStrings
gr()

f(d::Real, r::Real) = sqrt((d - r)*(d + r)) + r
dist = linspace(2,3,300)
rad = linspace(0,2,300)
surface(dist,rad,f,xlabel=L"d",ylabel=L"r_{opp}",zlabel="way")


text:
printstyled("  get_seis_channels\n", color=:light_green)
S = randSeisData(20, s=1.0)
S = S[get_seis_channels(S)]
@test get_seis_channels(S, chans=1) == get_seis_channels(S, chans=1:1) == get_seis_channels(S, chans=[1]) == [1]
c1 = collect(1:S.n)
filt_seis_chans!(c1, S)
@test c1 == 1:S.n

C1 = randSeisChannel(s=true)
C2 = randSeisChannel(s=false, c=true)
C2.id = "...YYY"
S = SeisData(C1, C2, S)
chans = collect(1:S.n)
filt_seis_chans!(chans, S)
@test (2 in chans) == false