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

import Base: convert, show

abstract Ty

type TypeName
    name::Symbol
    abs::Bool
    super::Ty    # actually TagT
    # arity
    # representation
    # mutable
    TypeName(name, abs, super) = new(name, abs, super)
    TypeName(name, abs) = new(name, abs)
end

show(io::IO, x::TypeName) = print(io, x.name)

type TagT <: Ty
    name::TypeName
    params
    vararg::Bool
    TagT(n, p, v=false) = new(n, p, v)
end

function show(io::IO, t::TagT)
    print(io, t.name.name)
    isempty(t.params) && return
    print(io, '{')
    l = length(t.params)
    for i=1:l
        show(io, t.params[i])
        if i == l && t.vararg
            print(io, "...")
        elseif i < l
            print(io, ",")
        end
    end
    print(io, '}')
end

type BottomTy <: Ty
end

show(io::IO, ::BottomTy) = print(io, "BottomT")

type UnionT <: Ty
    a
    b
    UnionT() = BottomT
    UnionT(a) = a
    UnionT(a, b) = new(a, b)
    UnionT(a, ts...) = new(a, UnionT(ts...))
end

function show(io::IO, t::UnionT)
    print(io, "UnionT(")
    while true
        show(io, t.a)
        print(io, ",")
        if isa(t.b, UnionT)
            t = t.b
        else
            show(io, t.b)
            break
        end
    end
    print(io, ')')
end

type Var
    name::Symbol
    lb
    ub
    Var(n, lb=BottomT, ub=AnyT) = new(n, lb, ub)
end

function show_var_bounds(io::IO, v::Var)
    if v.lb !== BottomT
        if v.ub === AnyT
            print(io, v.name)
            print(io, ">:")
            show(io, v.lb)
            return
        end
        show(io, v.lb)
        print(io, "<:")
    end
    print(io, v.name)
    if v.ub !== AnyT
        print(io, "<:")
        show(io, v.ub)
    end
end

show(io::IO, v::Var) = print(io, v.name)

type UnionAllT <: Ty
    var::Var
    T
    UnionAllT(v::Var, t) = new(v, t)
    UnionAllT(v::Var, t::Union{Type,Tuple}) = new(v, convert(Ty, t))
end

function show(io::IO, x::UnionAllT)
    print(io, "(@UnionAll ")
    show_var_bounds(io, x.var)
    print(io, " ")
    show(io, x.T)
    print(io, ")")
end

# Any, Bottom, and Tuple

const AnyT = TagT(TypeName(:Any,true), ())
AnyT.name.super = AnyT

const BottomT = BottomTy()

const TupleName = TypeName(:Tuple,false,AnyT)
const TupleT = TagT(TupleName, (AnyT,), true)
tupletype(xs...) = inst(TupleName, xs...)
vatype(xs...) = (t = inst(TupleName, xs...); t.vararg = true; t)

# type application

inst(typename::TypeName, params...) = TagT(typename, params)

inst(t::TagT) = t

inst(t::UnionAllT, param) = subst(t.T, Dict{Any,Any}(t.var => param))
inst(t::UnionAllT, param, rest...) = inst(inst(t,param), rest...)

supertype(t::TagT) = t.name===TupleName ? AnyT : inst(t.name.super, t.params...)

extend(d::Dict, k, v) = (x = copy(d); x[k]=v; x)

subst(t::TagT,    env) = t===AnyT ? t : TagT(t.name, map(x->subst(x,env), t.params), t.vararg)
subst(t::UnionT,  env) = UnionT(subst(t.a,env), subst(t.b,env))
subst(t::Var,     env) = get(env, t, t)
subst(t::UnionAllT, env) = (assert(!haskey(env, t.var));
                            newVar = Var(t.var.name, subst(t.var.lb, env), subst(t.var.ub, env));
                            UnionAllT(newVar, subst(t.T, extend(env, t.var, newVar))))
subst(t, env) = t

# subtype

isequal_type(x, y) = issub(x, y) && issub(y, x)

type Bounds
    # record current lower and upper bounds of a Var
    # right: whether this Var is on the right-hand side of A <: B
    lb
    ub
    right::Bool
end

type UnionState
    depth::Int           # number of union decision points we're inside
    more::Bool           # new union found; need to grow stack
    stack::Vector{Bool}  # stack of decisions
    UnionState() = new(1,0,Bool[])
end

type Env
    vars::Dict{Var,Bounds}
    Lunions::UnionState
    Runions::UnionState
    outer::Bool
    Env() = new(Dict{Var,Bounds}(), UnionState(), UnionState(), true)
end

function issub_env(x, y)
    e = Env()
    ans = forall_exists_issub(x, y, e)
    ans, e.vars
end
issub(x, y) = issub_env(x, y)[1]
issub(x, y, env) = (x === y)
issub(x::Ty, y::Ty, env) = (x === y) || x === BottomT

function forall_exists_issub(x, y, env)
    # for all combinations of elements from Unions on the left, there must
    # exist a combination of elements from Unions on the right that makes
    # issub() true. Unions in invariant position are on both the left and
    # the right in this formula.
    sub = exists_issub(x, y, env)

    if sub && env.Lunions.more
        push!(env.Lunions.stack, false)
        sub = forall_exists_issub(x, y, env)
        if sub
            env.Lunions.stack[end] = true
            sub = forall_exists_issub(x, y, env)
        end
        pop!(env.Lunions.stack)
    end
    return sub
end

function exists_issub(x, y, env)
    env.Lunions.depth = env.Runions.depth = 1
    env.Lunions.more = env.Runions.more = false

    found = issub(x, y, env)

    if env.Lunions.more
        # return up to forall_exists_issub. the recursion must have this shape:
        # ∀₁         ∀₁
        #   ∃₁  =>     ∀₂
        #                ...
        #                ∃₁
        #                  ∃₂
        return true
    end
    if env.Runions.more
        push!(env.Runions.stack, false)
        found = exists_issub(x, y, env)
        if !found
            env.Runions.stack[end] = true
            found = exists_issub(x, y, env)
        end
        pop!(env.Runions.stack)
    end
    return found
end

function issub_union(t, u::UnionT, env, R, state::UnionState)
    env.outer = false
    if state.depth > length(state.stack)
        # indicate that stack needs to grow
        state.more = true
        return true
    end
    ui = state.stack[state.depth]; state.depth += 1
    choice = getfield(u, 1+ui)
    return R ? issub(t, choice, env) : issub(choice, t, env)
end

issub(a::UnionT, b::UnionT, env) = a === b || issub_union(a, b, env, true, env.Runions)
issub(a::UnionT, b::Ty, env) = b===AnyT || issub_union(b, a, env, false, env.Lunions)
issub(a::Ty, b::UnionT, env) =
    a===BottomT || a===b.a || a===b.b || issub_union(a, b, env, true, env.Runions)

# handle vars before unions
issub(a::UnionT, b::Var, env) = var_gt(b, a, env)
issub(a::Var, b::UnionT, env) = var_lt(a, b, env)

function issub(a::TagT, b::TagT, env)
    env.outer = false
    a === b && return true
    b === AnyT && return true
    a === AnyT && return false
    if a.name !== b.name
        return issub(supertype(a), b, env)
    end
    if a.name === TupleName
        va, vb = a.vararg, b.vararg
        la, lb = length(a.params), length(b.params)
        ai = bi = 1
        while ai <= la
            bi > lb && return false
            !issub(a.params[ai], b.params[bi], env) && return false
            ai += 1
            if bi < lb || !vb
                bi += 1
            end
        end
        return (la==lb && va==vb) || (vb && (la >= (va ? lb : lb-1)))
    end
    for i = 1:length(a.params)
        ai, bi = a.params[i], b.params[i]
        # use issub in both directions to test equality
        if !(ai===bi || (issub(ai, bi, env) && issub(bi, ai, env)))
            return false
        end
    end
    return true
end

function join(a,b)
    (a===BottomT || b===AnyT || a === b) && return b
    (b===BottomT || a===AnyT) && return a
    UnionT(a,b)
end

issub(a::Ty, b::Var, env) = var_gt(b, a, env)
issub(a::Var, b::Ty, env) = var_lt(a, b, env)
function issub(a::Var, b::Var, env)
    env.outer = false
    a === b && return true
    aa = env.vars[a]; bb = env.vars[b]
    if aa.right
        # this is a bit odd, but seems necessary to make this case work:
        # (@UnionAll x<:T<:x RefT{RefT{T}}) == RefT{@UnionAll x<:T<:x RefT{T}}
        bb.right && return issub(bb.ub, bb.lb, env)
        return var_lt(a, b, env)
    else
        if !bb.right   # check ∀a,b . a<:b
            # the bounds of left-side variables never change, and can only lead
            # to other left-side variables, so using || here is safe.
            return issub(aa.ub, b, env) || issub(a, bb.lb, env)
        end
        return var_gt(b, a, env)
    end
end

# issub, but taking apart unions before handling vars
issub_ufirst(a::UnionT, b::Var, env) = issub_union(b, a, env, false, env.Lunions)
issub_ufirst(a::Var, b::UnionT, env) = a===b.a || a===b.b || issub_union(a, b, env, true, env.Runions)
issub_ufirst(a, b, env) = issub(a, b, env)

function var_lt(b::Var, a::Union{Ty,Var}, env)
    env.outer = false
    bb = env.vars[b]
    #println("$b($(bb.lb),$(bb.ub)) <: $a")
    !bb.right && return issub_ufirst(bb.ub, a, env)  # check ∀b . b<:a
    !issub_ufirst(bb.lb, a, env) && return false
    # for contravariance we would need to compute a meet here, but
    # because of invariance bb.ub ⊓ a == a here always. however for this
    # to work we need to compute issub(left,right) before issub(right,left),
    # since otherwise the issub(a, bb.ub) check in var_gt becomes vacuous.
    bb.ub = a  # meet(bb.ub, a)
    return true
end

function var_gt(b::Var, a::Union{Ty,Var}, env)
    env.outer = false
    bb = env.vars[b]
    #println("$b($(bb.lb),$(bb.ub)) >: $a")
    !bb.right && return issub_ufirst(a, bb.lb, env)  # check ∀b . b>:a
    !issub_ufirst(a, bb.ub, env) && return false
    bb.lb = join(bb.lb, a)
    return true
end

function rename(t::UnionAllT)
    v = Var(t.var.name, t.var.lb, t.var.ub)
    UnionAllT(v, inst(t,v))
end

function issub_unionall(t::Ty, u::UnionAllT, env, R)
    outer = env.outer
    haskey(env.vars, u.var) && (u = rename(u))
    env.vars[u.var] = Bounds(u.var.lb, u.var.ub, R)
    ans = R ? issub(t, u.T, env) : issub(u.T, t, env)
    !outer && delete!(env.vars, u.var)
    return ans
end

issub(a::UnionAllT, b::UnionAllT, env) = a === b || issub_unionall(a, b, env, true)
issub(a::UnionT, b::UnionAllT, env) = issub_unionall(a, b, env, true)
issub(a::UnionAllT, b::UnionT, env) = issub_unionall(b, a, env, false)
issub(a::Ty, b::UnionAllT, env) = a === BottomT || issub_unionall(a, b, env, true)
issub(a::UnionAllT, b::Ty, env) = b === AnyT || issub_unionall(b, a, env, false)

# convenient syntax

macro UnionAll(var, expr)
    lb = :BottomT
    ub = :AnyT
    if isa(var,Expr) && var.head === :comparison
        if length(var.args) == 3
            v = var.args[1]
            if var.args[2] == :(<:)
                ub = esc(var.args[3])
            elseif var.args[2] == :(>:)
                lb = esc(var.args[3])
            else
                error("invalid bounds in UnionAll")
            end
        elseif length(var.args) == 5
            v = var.args[3]
            if var.args[2] == var.args[4] == :(<:)
                lb = esc(var.args[1])
                ub = esc(var.args[5])
            else
                error("invalid bounds in UnionAll")
            end
        else
            error("invalid bounds in UnionAll")
        end
    elseif isa(var,Expr) && var.head === :(<:)
        v = var.args[1]
        ub = esc(var.args[2])
    elseif isa(var,Expr) && var.head === :(>:)
        v = var.args[1]
        lb = esc(var.args[2])
    elseif !isa(var,Symbol)
        error("invalid variable in UnionAll")
    else
        v = var
    end
    quote
        let $(esc(v)) = Var($(Expr(:quote,v)), $lb, $ub)
            UnionAllT($(esc(v)), $(esc(expr)))
        end
    end
end

# translating from existing julia types

const tndict = ObjectIdDict()

xlate(t) = xlate(t, ObjectIdDict())

xlate(t, env) = t

function xlate(t::Union, env)
    if t === Union{}
        return BottomT
    end
    UnionT(map(x->xlate(x,env), t.types)...)
end

function xlate(t::Tuple, env)
    if length(t) == 0
        return inst(TupleName)
    end
    va = Base.isvarargtype(t[end])
    ts = map(x->(Base.isvarargtype(x) ? xlate(x.parameters[1],env) : xlate(x,env)), t)
    tnew = inst(TupleName, ts...)
    tnew.vararg = va
    tnew
end

function xlate(t::TypeVar, env)
    if haskey(env, t)
        return env[t]
    end
    v = Var(t.name, xlate(t.lb,env), xlate(t.ub,env))
    env[t] = v
    v
end

function xlate(t::DataType, env)
    if t === Any
        return AnyT
    elseif t <: Tuple
        return xlate((t.parameters...,), env)
    elseif !haskey(tndict,t.name)
        para = map(x->xlate(x,env), t.name.primary.parameters)  # adds tvars to env
        sup = xlate(t.name.primary.super, env)
        for i = length(para):-1:1
            sup = UnionAllT(para[i], sup)
        end
        tn = TypeName(t.name.name, t.abstract, sup)
        tndict[t.name] = tn
    else
        tn = tndict[t.name]
    end
    inst(tn, map(x->xlate(x,env), t.parameters)...)
end

convert(::Type{Ty}, t::Union{Type,Tuple}) = xlate(t)
convert(::Type{Ty}, t::TypeVar) = xlate(t)

issub(a::Union{Type,Tuple}, b::Union{Type,Tuple}) = issub(xlate(a), xlate(b))
issub(a::Ty  , b::Union{Type,Tuple})              = issub(a       , xlate(b))
issub(a::Union{Type,Tuple}, b::Ty  )              = issub(xlate(a), b)
issub_env(a::Union{Type,Tuple}, b::Union{Type,Tuple}) = issub_env(xlate(a), xlate(b))
issub_env(a::Ty  , b::Union{Type,Tuple})              = issub_env(a       , xlate(b))
issub_env(a::Union{Type,Tuple}, b::Ty  )              = issub_env(xlate(a), b)

tt(ts...) = Tuple{ts...}
vt(ts...) = Tuple{ts[1:end-1]..., Vararg{ts[end]}}

# tests

AbstractArrayT =
    let AbstractArrayName = TypeName(:AbstractArray, true, @UnionAll T @UnionAll N AnyT)
        @UnionAll T @UnionAll N inst(AbstractArrayName, T, N)
    end

ArrayT =
    let ArrayName = TypeName(:Array, false, @UnionAll T @UnionAll N inst(AbstractArrayT, T, N))
        @UnionAll T @UnionAll N inst(ArrayName, T, N)
    end

PairT = let PairName = TypeName(:Pair, false, @UnionAll A @UnionAll B AnyT)
    @UnionAll A @UnionAll B inst(PairName, A, B)
end

RefT = let RefName = TypeName(:Ref, false, @UnionAll T AnyT)
    @UnionAll T inst(RefName, T)
end

tndict[AbstractArray.name] = AbstractArrayT.T.T.name
tndict[Array.name] = ArrayT.T.T.name
tndict[Pair.name] = PairT.T.T.name

using Base.Test

issub_strict(x,y) = issub(x,y) && !issub(y,x)

# level 1: no varags, union, UnionAll
function test_1()
    @test issub_strict(Int, Integer)
    @test issub_strict(Array{Int,1}, AbstractArray{Int,1})

    @test isequal_type(Int, Int)
    @test isequal_type(Integer, Integer)
    @test isequal_type(Array{Int,1}, Array{Int,1})
    @test isequal_type(AbstractArray{Int,1}, AbstractArray{Int,1})

    @test issub_strict((Int,Int), (Integer,Integer))
    @test issub_strict((Array{Int,1},), (AbstractArray{Int,1},))

    @test isequal_type((Integer,Integer), (Integer,Integer))

    @test !issub((Int,Int), (Int,))
    @test !issub((Int,), (Integer,Integer))
end

# level 2: varargs
function test_2()
    @test issub_strict((Int,Int), vt(Int,))
    @test issub_strict((Int,Int), vt(Int,Int,))
    @test issub_strict((Int,Int), vt(Int,Integer,))
    @test issub_strict((Int,Int), vt(Int,Int,Integer,))
    @test issub_strict(vt(Int,Int), vt(Int,))
    @test issub_strict(vt(Int,Int,Int), vt(Int,))
    @test issub_strict(vt(Int,Int,Int), vt(Integer,Int,))
    @test issub_strict(vt(Int,), vt(Any,))
    @test issub_strict((), vt(Any,))

    @test isequal_type(vt(Int,), vt(Int,))
    @test isequal_type(vt(Integer,), vt(Integer,))

    @test !issub((), vt(Int, Int))
    @test !issub((Int,), vt(Int, Int, Int))

    @test !issub((Int, (Real, Integer)), vt(Int))
end

function test_no_diagonal()
    # from test_3
    @test issub((@UnionAll T @UnionAll S tupletype(T,S)), (@UnionAll T tupletype(T,T)))

    @test isequal_type((@UnionAll T tupletype(T,T)), (@UnionAll T @UnionAll S tupletype(T,S)))

    @test issub(tupletype(inst(ArrayT,Ty(Integer),1), Ty(Int)),
                (@UnionAll T<:Ty(Integer) tupletype(inst(ArrayT,T,1),T)))

    @test issub(Ty((Int,AbstractString,Vector{Any})),
                @UnionAll T tupletype(T, T, inst(ArrayT,T,1)))

    @test isequal_type(Ty(Array{Tuple{Integer,Integer},1}),
                       inst(ArrayT, (@UnionAll T<:Ty(Integer) tupletype(T,T)), 1))

    @test issub(Ty((Float32,Array{Real,1})),
                @UnionAll T<:Ty(Real) @UnionAll S<:inst(AbstractArrayT,T,1) tupletype(T,S))

    @test isequal_type((@UnionAll T tupletype(inst(RefT,T), T)),
                       (@UnionAll T @UnionAll S<:T tupletype(inst(RefT,T),S)))
    @test isequal_type((@UnionAll T tupletype(inst(RefT,T), T)),
                       (@UnionAll T @UnionAll S<:T @UnionAll R<:S tupletype(inst(RefT,T),R)))

    @test  issub((@UnionAll S<:Ty(Int) (@UnionAll R<:Ty(AbstractString) tupletype(S,R,Ty(Vector{Any})))),
                 (@UnionAll T tupletype(T, T, inst(ArrayT,T,1))))
end

# level 3: UnionAll
function test_3()
    @test issub_strict(Ty(Array{Int,1}), @UnionAll T inst(ArrayT, T, 1))
    @test issub_strict((@UnionAll T inst(PairT,T,T)), (@UnionAll T @UnionAll S inst(PairT,T,S)))
    @test issub(inst(PairT,Ty(Int),Ty(Int8)), (@UnionAll T @UnionAll S inst(PairT,T,S)))
    @test issub(inst(PairT,Ty(Int),Ty(Int8)), (@UnionAll S inst(PairT,Ty(Int),S)))

    @test !issub((@UnionAll T<:Ty(Real) T), (@UnionAll T<:Ty(Integer) T))

    # diagonal
    #@test issub_strict((@UnionAll T tupletype(T,T)), (@UnionAll T @UnionAll S tupletype(T,S)))
    @test issub((@UnionAll T tupletype(T,T)), (@UnionAll T @UnionAll S tupletype(T,S)))
    @test isequal_type((@UnionAll T tupletype(T,T)), (@UnionAll R tupletype(R,R)))
    # diagonal
    #@test !issub(tupletype(Ty(Real),Ty(Real)), @UnionAll T<:Ty(Real) tupletype(T,T))

    @test !issub((@UnionAll T<:Ty(Integer) @UnionAll S<:Ty(Number) (T,S)),
                 (@UnionAll T<:Ty(Integer) @UnionAll S<:Ty(Number) (S,T)))

    AUA = inst(ArrayT, (@UnionAll T inst(ArrayT,T,1)), 1)
    UAA = (@UnionAll T inst(ArrayT, inst(ArrayT,T,1), 1))

    @test !issub(AUA, UAA)
    @test !issub(UAA, AUA)
    @test !isequal_type(AUA, UAA)

    @test issub_strict((@UnionAll T Int), (@UnionAll T<:Ty(Integer) Integer))

    @test isequal_type((@UnionAll T @UnionAll S tupletype(T, tupletype(S))),
                       (@UnionAll T tupletype(T, @UnionAll S tupletype(S))))

    @test !issub((@UnionAll T inst(PairT,T,T)), inst(PairT,Ty(Int),Ty(Int8)))
    @test !issub((@UnionAll T inst(PairT,T,T)), inst(PairT,Ty(Int),Ty(Int)))

    @test isequal_type((@UnionAll T tupletype(T)), tupletype(AnyT))
    @test isequal_type((@UnionAll T<:Ty(Real) tupletype(T)), tupletype(Ty(Real)))

    #@test !issub(tupletype(inst(ArrayT,Ty(Integer),1), Ty(Int)),
    #             (@UnionAll T<:Ty(Integer) tupletype(inst(ArrayT,T,1),T)))
    @test  issub(tupletype(inst(ArrayT,Ty(Integer),1), Ty(Int)),
                 (@UnionAll T<:Ty(Integer) @UnionAll S<:T tupletype(inst(ArrayT,T,1),S)))

    @test !issub(tupletype(inst(ArrayT,Ty(Integer),1), Ty(Real)),
                 (@UnionAll T<:Ty(Integer) tupletype(inst(ArrayT,T,1),T)))

    @test !issub(Ty((Int,AbstractString,Vector{Integer})),
                 @UnionAll T tupletype(T, T, inst(ArrayT,T,1)))
    @test !issub(Ty((AbstractString,Int,Vector{Integer})),
                 @UnionAll T tupletype(T, T, inst(ArrayT,T,1)))
    @test !issub(Ty((Int,AbstractString,Vector{Tuple{Integer}})),
                 @UnionAll T tupletype(T,T,inst(ArrayT,tupletype(T),1)))

    #@test !issub(Ty((Int,AbstractString,Vector{Any})),
    #             @UnionAll T tupletype(T, T, inst(ArrayT,T,1)))

    @test isequal_type(Ty(Array{Int,1}), inst(ArrayT, (@UnionAll T<:Ty(Int) T), 1))
    @test isequal_type(Ty(Array{Tuple{Any},1}), inst(ArrayT, (@UnionAll T tupletype(T)), 1))

    @test isequal_type(Ty(Array{Tuple{Int,Int},1}),
                       inst(ArrayT, (@UnionAll T<:Ty(Int) tupletype(T,T)), 1))
    @test !issub(Ty(Array{Tuple{Int,Integer},1}),
                 inst(ArrayT, (@UnionAll T<:Ty(Integer) tupletype(T,T)), 1))

    @test !issub(inst(PairT,Ty(Int),Ty(Int8)), (@UnionAll T inst(PairT,T,T)))

    @test !issub(tupletype(inst(ArrayT,Ty(Int),1), Ty(Integer)),
                 (@UnionAll T<:Ty(Integer) tupletype(inst(ArrayT,T,1),T)))

    @test !issub(tupletype(Ty(Integer), inst(ArrayT,Ty(Int),1)),
                 (@UnionAll T<:Ty(Integer) tupletype(T, inst(ArrayT,T,1))))

    @test !issub(Ty(Array{Array{Int,1},Integer}),
                 (@UnionAll T inst(ArrayT,inst(ArrayT,T,1),T)))

    @test issub(Ty(Array{Array{Int,1},Int}),
                (@UnionAll T inst(ArrayT,inst(ArrayT,T,1),T)))

    @test  issub(Ty((Integer,Int)), @UnionAll T<:Ty(Integer) @UnionAll S<:T tupletype(T,S))
    @test !issub(Ty((Integer,Int)), @UnionAll T<:Ty(Int) @UnionAll S<:T tupletype(T,S))
    @test !issub(Ty((Integer,Int)), @UnionAll T<:Ty(AbstractString) @UnionAll S<:T tupletype(T,S))

    @test issub(Ty((Float32,Array{Float32,1})),
                @UnionAll T<:Ty(Real) @UnionAll S<:inst(AbstractArrayT,T,1) tupletype(T,S))

    @test !issub(Ty((Float32,Array{Float64,1})),
                 @UnionAll T<:Ty(Real) @UnionAll S<:inst(AbstractArrayT,T,1) tupletype(T,S))

    #@test !issub(Ty((Float32,Array{Real,1})),
    #             @UnionAll T<:Ty(Real) @UnionAll S<:inst(AbstractArrayT,T,1) tupletype(T,S))

    @test !issub(Ty((Number,Array{Real,1})),
                 @UnionAll T<:Ty(Real) @UnionAll S<:inst(AbstractArrayT,T,1) tupletype(T,S))

    @test issub((@UnionAll Ty(Int)<:T<:Ty(Integer) T), @UnionAll T<:Ty(Real) T)
    @test issub((@UnionAll Ty(Int)<:T<:Ty(Integer) inst(ArrayT,T,1)),
                (@UnionAll T<:Ty(Real) inst(ArrayT,T,1)))

    @test issub((@UnionAll Ty(Int)<:T<:Ty(Integer) T),
                (@UnionAll Ty(Integer)<:T<:Ty(Real) T))
    @test !issub((@UnionAll Ty(Int)<:T<:Ty(Integer) inst(ArrayT,T,1)),
                 (@UnionAll Ty(Integer)<:T<:Ty(Real) inst(ArrayT,T,1)))

    X = (@UnionAll T<:Ty(Real) @UnionAll S<:inst(AbstractArrayT,T,1) tupletype(T,S))
    Y = (@UnionAll A<:Ty(Real) @UnionAll B<:inst(AbstractArrayT,A,1) tupletype(A,B))
    @test isequal_type(X,Y)
    Z = (@UnionAll A<:Ty(Real) @UnionAll B<:inst(AbstractArrayT,A,1) tupletype(Ty(Real),B))
    @test issub_strict(X,Z)

    @test issub_strict((@UnionAll T @UnionAll S<:T inst(PairT,T,S)),
                       (@UnionAll T @UnionAll S    inst(PairT,T,S)))
    @test issub_strict((@UnionAll T @UnionAll S>:T inst(PairT,T,S)),
                       (@UnionAll T @UnionAll S    inst(PairT,T,S)))

    #@test issub_strict((@UnionAll T tupletype(inst(RefT,T), T)),
    #                   (@UnionAll T @UnionAll S<:T tupletype(inst(RefT,T),S)))
    #@test issub_strict((@UnionAll T tupletype(inst(RefT,T), T)),
    #                   (@UnionAll T @UnionAll S<:T @UnionAll R<:S tupletype(inst(RefT,T),R)))
    @test isequal_type((@UnionAll T tupletype(inst(RefT,T), T)),
                       (@UnionAll T @UnionAll T<:S<:T tupletype(inst(RefT,T),S)))
    @test issub_strict((@UnionAll T tupletype(inst(RefT,T), T)),
                       (@UnionAll T @UnionAll S>:T tupletype(inst(RefT,T),S)))
end

# level 4: Union
function test_4()
    @test isequal_type(UnionT(BottomT,BottomT), BottomT)

    @test issub_strict(Int, Union{Int,AbstractString})
    @test issub_strict(Union{Int,Int8}, Integer)

    @test isequal_type(Union{Int,Int8}, Union{Int,Int8})

    @test isequal_type(UnionT(Ty(Int),Ty(Integer)), Ty(Integer))

    @test isequal_type(tt(Union{Int,Int8},Int16), Union{tt(Int,Int16),tt(Int8,Int16)})

    @test issub_strict((Int,Int8,Int), vt(Union{Int,Int8},))
    @test issub_strict((Int,Int8,Int), vt(Union{Int,Int8,Int16},))

    # nested unions
    @test !issub(UnionT(Ty(Int),inst(RefT,UnionT(Ty(Int),Ty(Int8)))),
                 UnionT(Ty(Int),inst(RefT,UnionT(Ty(Int8),Ty(Int16)))))

    A = Ty(Int);   B = Ty(Int8)
    C = Ty(Int16); D = Ty(Int32)
    @test  issub(UnionT(UnionT(A,UnionT(A,UnionT(B,C))), UnionT(D,BottomT)),
                 UnionT(UnionT(A,B),UnionT(C,UnionT(B,D))))
    @test !issub(UnionT(UnionT(A,UnionT(A,UnionT(B,C))), UnionT(D,BottomT)),
                 UnionT(UnionT(A,B),UnionT(C,UnionT(B,A))))

    @test isequal_type(UnionT(UnionT(A,B,C), UnionT(D)),  UnionT(A,B,C,D))
    @test isequal_type(UnionT(UnionT(A,B,C), UnionT(D)),  UnionT(A,UnionT(B,C),D))
    @test isequal_type(UnionT(UnionT(UnionT(UnionT(A)),B,C), UnionT(D)),
                       UnionT(A,UnionT(B,C),D))

    @test issub_strict(UnionT(UnionT(A,C), UnionT(D)),  UnionT(A,B,C,D))

    @test !issub(UnionT(UnionT(A,B,C), UnionT(D)),  UnionT(A,C,D))

    # obviously these unions can be simplified, but when they aren't there's trouble
    X = UnionT(UnionT(A,B,C),UnionT(A,B,C),UnionT(A,B,C),UnionT(A,B,C),
               UnionT(A,B,C),UnionT(A,B,C),UnionT(A,B,C),UnionT(A,B,C))
    Y = UnionT(UnionT(D,B,C),UnionT(D,B,C),UnionT(D,B,C),UnionT(D,B,C),
               UnionT(D,B,C),UnionT(D,B,C),UnionT(D,B,C),UnionT(A,B,C))
    @test issub_strict(X,Y)
end

# level 5: union and UnionAll
function test_5()
    u = Ty(Union{Int8,Int})

    @test issub(Ty((AbstractString,Array{Int,1})),
                (@UnionAll T UnionT(tupletype(T,inst(ArrayT,T,1)),
                                    tupletype(T,inst(ArrayT,Ty(Int),1)))))

    @test issub(Ty((Union{Vector{Int},Vector{Int8}},)),
                @UnionAll T tupletype(inst(ArrayT,T,1),))

    @test !issub(Ty((Union{Vector{Int},Vector{Int8}},Vector{Int})),
                 @UnionAll T tupletype(inst(ArrayT,T,1), inst(ArrayT,T,1)))

    @test !issub(Ty((Union{Vector{Int},Vector{Int8}},Vector{Int8})),
                 @UnionAll T tupletype(inst(ArrayT,T,1), inst(ArrayT,T,1)))

    @test !issub(Ty(Vector{Int}), @UnionAll T>:u inst(ArrayT,T,1))
    @test  issub(Ty(Vector{Integer}), @UnionAll T>:u inst(ArrayT,T,1))
    @test  issub(Ty(Vector{Union{Int,Int8}}), @UnionAll T>:u inst(ArrayT,T,1))

    @test issub((@UnionAll Ty(Int)<:T<:u inst(ArrayT,T,1)),
                (@UnionAll Ty(Int)<:T<:u inst(ArrayT,T,1)))

    # with varargs
    @test !issub(inst(ArrayT,tupletype(inst(ArrayT,Ty(Int)),inst(ArrayT,Ty(Vector{Int16})),inst(ArrayT,Ty(Vector{Int})),inst(ArrayT,Ty(Int)))),
                 @UnionAll T<:(@UnionAll S vatype(UnionT(inst(ArrayT,S),inst(ArrayT,inst(ArrayT,S,1))))) inst(ArrayT,T))

    @test issub(inst(ArrayT,tupletype(inst(ArrayT,Ty(Int)),inst(ArrayT,Ty(Vector{Int})),inst(ArrayT,Ty(Vector{Int})),inst(ArrayT,Ty(Int)))),
                @UnionAll T<:(@UnionAll S vatype(UnionT(inst(ArrayT,S),inst(ArrayT,inst(ArrayT,S,1))))) inst(ArrayT,T))

    @test !issub(tupletype(inst(ArrayT,Ty(Int)),inst(ArrayT,Ty(Vector{Int16})),inst(ArrayT,Ty(Vector{Int})),inst(ArrayT,Ty(Int))),
                 @UnionAll S vatype(UnionT(inst(ArrayT,S),inst(ArrayT,inst(ArrayT,S,1)))))

    @test issub(tupletype(inst(ArrayT,Ty(Int)),inst(ArrayT,Ty(Vector{Int})),inst(ArrayT,Ty(Vector{Int})),inst(ArrayT,Ty(Int))),
                @UnionAll S vatype(UnionT(inst(ArrayT,S),inst(ArrayT,inst(ArrayT,S,1)))))

    B = @UnionAll S<:u tupletype(S, tupletype(AnyT,AnyT,AnyT), inst(RefT,S))
    # these tests require renaming in issub_unionall
    @test  issub((@UnionAll T<:B tupletype(Ty(Int8), T, inst(RefT,Ty(Int8)))), B)
    @test !issub((@UnionAll T<:B tupletype(Ty(Int8), T, inst(RefT,T))),        B)

    # the `convert(Type{T},T)` pattern, where T is a Union
    # required changing priority of unions and vars
    @test issub(tupletype(inst(ArrayT,u,1),Ty(Int)),
                @UnionAll T tupletype(inst(ArrayT,T,1), T))
    @test issub(tupletype(inst(ArrayT,u,1),Ty(Int)),
                @UnionAll T @UnionAll S<:T tupletype(inst(ArrayT,T,1), S))
end

# tricky type variable lower bounds
function test_6()
    # diagonal
    #@test !issub((@UnionAll S<:Ty(Int) (@UnionAll R<:Ty(AbstractString) tupletype(S,R,Ty(Vector{Any})))),
    #             (@UnionAll T tupletype(T, T, inst(ArrayT,T,1))))
    @test issub((@UnionAll S<:Ty(Int) (@UnionAll R<:Ty(AbstractString) tupletype(S,R,Ty(Vector{Any})))),
                (@UnionAll T tupletype(T, T, inst(ArrayT,T,1))))

    @test !issub((@UnionAll S<:Ty(Int) (@UnionAll R<:Ty(AbstractString) tupletype(S,R,Ty(Vector{Integer})))),
                 (@UnionAll T tupletype(T, T, inst(ArrayT,T,1))))

    t = @UnionAll T tupletype(T,T,inst(RefT,T))
    @test isequal_type(t, rename(t))

    @test !issub((@UnionAll T tupletype(T,Ty(AbstractString),inst(RefT,T))),
                 (@UnionAll T tupletype(T,T,inst(RefT,T))))

    @test !issub((@UnionAll T tupletype(T,inst(RefT,T),Ty(AbstractString))),
                 (@UnionAll T tupletype(T,inst(RefT,T),T)))

    i = Ty(Int); ai = Ty(Integer)
    @test isequal_type((@UnionAll i<:T<:i   inst(RefT,T)), inst(RefT,i))
    @test isequal_type((@UnionAll ai<:T<:ai inst(RefT,T)), inst(RefT,ai))

    # Pair{T,S} <: Pair{T,T} can be true with certain bounds
    @test issub_strict((@UnionAll i<:T<:i @UnionAll i<:S<:i inst(PairT,T,S)),
                       @UnionAll T inst(PairT,T,T))

    @test issub_strict(tupletype(i, inst(RefT,i)),
                       (@UnionAll T @UnionAll S<:T tupletype(S,inst(RefT,T))))

    @test !issub(tupletype(Ty(Real), inst(RefT,i)),
                 (@UnionAll T @UnionAll S<:T tupletype(S,inst(RefT,T))))

    # S >: T
    @test issub_strict(tupletype(Ty(Real), inst(RefT,i)),
                       (@UnionAll T @UnionAll S>:T tupletype(S,inst(RefT,T))))

    @test !issub(tupletype(inst(RefT,i), inst(RefT,ai)),
                 (@UnionAll T @UnionAll S>:T tupletype(inst(RefT,S),inst(RefT,T))))

    @test issub_strict(tupletype(inst(RefT,Ty(Real)), inst(RefT,ai)),
                       (@UnionAll T @UnionAll S>:T tupletype(inst(RefT,S),inst(RefT,T))))

    @test issub_strict(tupletype(Ty(Real), inst(RefT,tupletype(i))),
                       (@UnionAll T @UnionAll S>:T tupletype(S,inst(RefT,tupletype(T)))))

    @test !issub(tupletype(inst(RefT,tupletype(i)), inst(RefT,tupletype(ai))),
                 (@UnionAll T @UnionAll S>:T tupletype(inst(RefT,tupletype(S)),inst(RefT,tupletype(T)))))

    @test issub_strict(tupletype(inst(RefT,tupletype(Ty(Real))), inst(RefT,tupletype(ai))),
                       (@UnionAll T @UnionAll S>:T tupletype(inst(RefT,tupletype(S)),inst(RefT,tupletype(T)))))

    # (@UnionAll x<:T<:x Q{T}) == Q{x}
    @test isequal_type(inst(RefT,inst(RefT,i)),
                       inst(RefT,@UnionAll i<:T<:i inst(RefT,T)))
    @test isequal_type((@UnionAll i<:T<:i inst(RefT,inst(RefT,T))),
                       inst(RefT,@UnionAll i<:T<:i inst(RefT,T)))
    @test !issub((@UnionAll i<:T<:i inst(RefT,inst(RefT,T))),
                 inst(RefT,@UnionAll T<:i inst(RefT,T)))

    u = Ty(Union{Int8,Int64})
    A = inst(RefT,BottomT)
    B = @UnionAll S<:u inst(RefT,S)
    @test issub(inst(RefT,B), @UnionAll A<:T<:B inst(RefT,T))

    C = @UnionAll S<:u S
    @test issub(inst(RefT,C), @UnionAll u<:T<:u inst(RefT,T))

    BB = @UnionAll S<:BottomT S
    @test issub(inst(RefT,B), @UnionAll BB<:U<:B inst(RefT,U))
end

# tests that don't pass yet
function test_failing()
    # TODO: S <: Array{T} cases
end

function test_all()
    test_1()
    test_2()
    test_no_diagonal()
    test_3()
    test_4()
    test_5()
    test_6()
    test_failing()
end

const menagerie =
    Any[BottomT, AnyT, Ty(Int), Ty(Int8), Ty(Integer), Ty(Real),
        Ty(Array{Int,1}), Ty(AbstractArray{Int,1}),
        Ty(vt(Int,Integer,)), Ty(vt(Integer,Int,)), Ty(()),
        Ty(Union{Int,Int8}),
        (@UnionAll T inst(ArrayT, T, 1)),
        (@UnionAll T inst(PairT,T,T)),
        (@UnionAll T @UnionAll S inst(PairT,T,S)),
        inst(PairT,Ty(Int),Ty(Int8)),
        (@UnionAll S inst(PairT,Ty(Int),S)),
        (@UnionAll T tupletype(T,T)),
        (@UnionAll T<:Ty(Integer) tupletype(T,T)),
        (@UnionAll T @UnionAll S tupletype(T,S)),
        (@UnionAll T<:Ty(Integer) @UnionAll S<:Ty(Number) (T,S)),
        (@UnionAll T<:Ty(Integer) @UnionAll S<:Ty(Number) (S,T)),
        inst(ArrayT, (@UnionAll T inst(ArrayT,T,1)), 1),
        (@UnionAll T inst(ArrayT, inst(ArrayT,T,1), 1)),
        inst(ArrayT, (@UnionAll T<:Ty(Int) T), 1),
        (@UnionAll T<:Ty(Real) @UnionAll S<:inst(AbstractArrayT,T,1) tupletype(T,S)),
        UnionT(Ty(Int),inst(RefT,UnionT(Ty(Int),Ty(Int8)))),
        (@UnionAll T UnionT(tupletype(T,inst(ArrayT,T,1)),
                            tupletype(T,inst(ArrayT,Ty(Int),1)))),
        ]

let new = Any[]
    # add variants of each type
    for T in menagerie
        push!(new, inst(RefT, T))
        push!(new, tupletype(T))
        push!(new, tupletype(T,T))
        push!(new, vatype(T))
        push!(new, @UnionAll S<:T S)
        push!(new, @UnionAll S<:T inst(RefT,S))
    end
    append!(menagerie, new)
end

function test_properties()
    x→y = !x || y
    ¬T = @UnionAll X>:T inst(RefT,X)

    for T in menagerie
        # top and bottom identities
        @test issub(BottomT, T)
        @test issub(T, AnyT)
        @test issub(T, BottomT) → isequal_type(T, BottomT)
        @test issub(AnyT, T) → isequal_type(T, AnyT)

        # unionall identity
        @test isequal_type(T, @UnionAll S<:T S)
        @test isequal_type(inst(RefT,T), @UnionAll T<:U<:T inst(RefT,U))

        # equality under renaming
        if isa(T, UnionAllT)
            @test isequal_type(T, rename(T))
        end

        # inequality under wrapping
        @test !isequal_type(T, inst(RefT,T))

        for S in menagerie
            issubTS = issub(T, S)
            # transitivity
            if issubTS
                for R in menagerie
                    if issub(S, R)
                        if !issub(T, R)
                            @show T
                            @show S
                            @show R
                        end
                        @test issub(T, R)  # issub(S, R) → issub(T, R)
                        @test issub(inst(RefT,S), @UnionAll T<:U<:R inst(RefT,U))
                    end
                end
            end

            # union subsumption
            @test isequal_type(T, UnionT(T,S)) → issub(S, T)

            # invariance
            @test isequal_type(T, S) == isequal_type(inst(RefT,T), inst(RefT,S))

            # covariance
            @test issubTS == issub(tupletype(T), tupletype(S))
            @test issubTS == issub(vatype(T), vatype(S))
            @test issubTS == issub(tupletype(T), vatype(S))

            # contravariance
            @test issubTS == issub(¬S, ¬T)
        end
    end
end

# function non_terminating_F()
#     # undecidable F_<: instance
#     ¬T = @ForAll α<:T α
#     θ = @ForAll α ¬(@ForAll β<:α ¬β)
#     a₀ = Var(:a₀, BottomT, θ)
#     issub(a₀, (@ForAll a₁<:a₀ ¬a₁))
# end

# attempt to implement non-terminating example from
# "On the Decidability of Subtyping with Bounded Existential Types"
function non_terminating_2()
    C = let CName = TypeName(:C, false, @UnionAll T AnyT)
        @UnionAll T inst(CName, T)
    end
    D = let DName = TypeName(:D, false, @UnionAll T AnyT)
        @UnionAll T inst(DName, T)
    end
    ¬T = @UnionAll X>:T inst(D,X)
    U = AnyT
    X = Var(:X, BottomT, ¬U)
    e = Env()
    e.vars[X] = Bounds(BottomT, ¬U, false)
    issub(X, ¬inst(C,X), e)
end

#=
Replacing type S with a union over var T(ᶜ)<:S can affect the subtype
relation in all possible ways:

    (C{Real},C{Real}) < @U T<:Real (C{T},C{T})
    (Real,Real)       = @U T<:Real (T,T)

    (C{Real},C{Real}) ! @U Tᶜ<:Real (C{T},C{T})  (*)  # neither <= nor >=
    (Real,Real)       > @U Tᶜ<:Real (T,T)

    (C{Real},) < @U T<:Real (C{T},)
    (Real,)    = @U T<:Real (T,)
    (C{Real},) ! @U Tᶜ<:Real (C{T},)  (*)
    (Real,)    = @U Tᶜ<:Real (T,)

cases (*) are probably not needed in practice.

transform all covariant occurrences of non-concrete types to
vars: (Real,Real) => @U T<:Real @U S<:Real (T,S)
and then treat all vars as invariant.
this can maybe even be done on the fly without actually making the
UnionAll types up-front.
in fact this may only need to be done on the left.
which implies right away that the algorithm works as-is if the left
side is a concrete type.

possible explanation:

1. match all tvars with invariant rule
2. imagine left side is concrete
3. therefore all (right-) vars will have to equal concrete types
   (this completes the base case)
4. inductive hypothesis: tvars match concrete types
5. "build up" an abstract type on the left over the initial concrete type
6. handle all abstract types T in covariant position on the left by treating
   them as variables V<:T instead
7. by the inductive hypothesis tvars continue to match concrete types

it seems you need to plug in the isconcrete predicate in step 6.

=#


text:
#=====================================================#
#== MATHEMATICS AND ARITHMETIC TEMPLATES FOR HUMANS ==#
#=====================================================#
module MATH

using LinearAlgebra

#===========================================#
#== first use case: function manipulation ==#
#===========================================#
#== general function sum ==#
macro ∑(bounds::UnitRange{Int64}, op::Function)
  return :( general_driver(+, $(esc(bounds)) , $(esc(op)) ) )
end

#== general function product ==#
macro ∏(bounds::UnitRange{Int64}, op::Function)
  return :( general_driver(*, $(esc(bounds)) , $(esc(op)) ) )
end

#== common driver ==#
function general_driver(funct::Function, bounds::UnitRange{Int64}, op::Function)
  #GC.enable(false)
  return reduce(funct, map(op, collect(bounds)))
  #GC.enable(true)
  #GC.gc()
end

#================================================#
#== second use case: array/matrix manipulation ==#
#================================================#
#== summations involving single arrays ==#
macro ∑(array)
  return :( sum_over_one( $(esc(array)) ) )
end

@inline function sum_over_one(array::T) where {T<:AbstractArray}
  return sum(array)
end

@inline function sum_over_one(range::UnitRange{T}) where {T<:Number}
  return sum(collect(range))
end

#== summations involving two arrays ==#
macro ∑(array_1, array_2)
  return :( sum_over_two( $(esc(array_1)), $(esc(array_2)) ) )
end

@inline function sum_over_two(array_1::T,
  array_2::T) where {T<:AbstractArray}

  return LinearAlgebra.dot(array_1,array_2)
end

@inline function sum_over_two(array_1::Array{T,1},
  array_2::Array{T,1}) where {T<:Number}

  return LinearAlgebra.dot(array_1,array_2)
end

@inline function sum_over_two(array_1::Array{T,2},
  array_2::Array{T,2}) where {T<:Number}

  return LinearAlgebra.dot(array_1,array_2)
end

#== quadratic transformation summation ==#
macro ∑(array_1,array_2,array_center)
  return :( quadratic_transform( $(esc(array_1)), $(esc(array_2)) ,
    $(esc(array_center)) ) )
end

@inline function quadratic_transform(array_1::Array{T},
  array_2::Array{T}, array_center::Array{T,2}) where {T<:Number}

  return transpose(array_2)*array_center*array_1
end

#== products involving single arrays ==#
macro ∏(array)
  return :( prod_over_one( $(esc(array)) ) )
end

@inline function prod_over_one(array::T) where {T<:AbstractArray}
  return prod(array)
end

@inline function prod_over_one(range::UnitRange{T}) where {T<:Number}
  return prod(collect(range))
end

export @∑
export @∏

end

ext: jl


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
    Arrow.Map

An `ArrowVector` where each element is a "map" of some kind, like a `Dict`.
"""
struct Map{T, O, A} <: ArrowVector{T}
    validity::ValidityBitmap
    offsets::Offsets{O}
    data::A
    ℓ::Int
    metadata::Union{Nothing, Dict{String, String}}
end

Base.size(l::Map) = (l.ℓ,)

@propagate_inbounds function Base.getindex(l::Map{T}, i::Integer) where {T}
    @boundscheck checkbounds(l, i)
    @inbounds lo, hi = l.offsets[i]
    if Base.nonmissingtype(T) !== T
        return l.validity[i] ? ArrowTypes.fromarrow(T, Dict(x.key => x.value for x in view(l.data, lo:hi))) : missing
    else
        return ArrowTypes.fromarrow(T, Dict(x.key => x.value for x in view(l.data, lo:hi)))
    end
end

keyvalues(KT, ::Missing) = missing
keyvalues(KT, x::AbstractDict) = [KT(k, v) for (k, v) in pairs(x)]

keyvaluetypes(::Type{NamedTuple{(:key, :value), Tuple{K, V}}}) where {K, V} = (K, V)

arrowvector(::MapKind, x::Map, i, nl, fi, de, ded, meta; kw...) = x

function arrowvector(::MapKind, x, i, nl, fi, de, ded, meta; largelists::Bool=false, kw...)
    len = length(x)
    validity = ValidityBitmap(x)
    ET = eltype(x)
    DT = Base.nonmissingtype(ET)
    KT = KeyValue{keytype(DT), valtype(DT)}
    VT = Vector{KT}
    T = DT !== ET ? Union{Missing, VT} : VT
    flat = ToList(T[keyvalues(KT, y) for y in x]; largelists=largelists)
    offsets = Offsets(UInt8[], flat.inds)
    data = arrowvector(flat, i, nl + 1, fi, de, ded, nothing; lareglists=largelists, kw...)
    K, V = keyvaluetypes(eltype(data))
    return Map{withmissing(ET, Dict{K, V}), eltype(flat.inds), typeof(data)}(validity, offsets, data, len, meta)
end

function compress(Z::Meta.CompressionType, comp, x::A) where {A <: Map}
    len = length(x)
    nc = nullcount(x)
    validity = compress(Z, comp, x.validity)
    offsets = compress(Z, comp, x.offsets.offsets)
    buffers = [validity, offsets]
    children = Compressed[]
    push!(children, compress(Z, comp, x.data))
    return Compressed{Z, A}(x, buffers, len, nc, children)
end

function makenodesbuffers!(col::Union{Map{T, O, A}, List{T, O, A}}, fieldnodes, fieldbuffers, bufferoffset, alignment) where {T, O, A}
    len = length(col)
    nc = nullcount(col)
    push!(fieldnodes, FieldNode(len, nc))
    @debug 1 "made field node: nodeidx = $(length(fieldnodes)), col = $(typeof(col)), len = $(fieldnodes[end].length), nc = $(fieldnodes[end].null_count)"
    # validity bitmap
    blen = nc == 0 ? 0 : bitpackedbytes(len, alignment)
    push!(fieldbuffers, Buffer(bufferoffset, blen))
    @debug 1 "made field buffer: bufferidx = $(length(fieldbuffers)), offset = $(fieldbuffers[end].offset), len = $(fieldbuffers[end].length), padded = $(padding(fieldbuffers[end].length, alignment))"
    # adjust buffer offset, make array buffer
    bufferoffset += blen
    blen = sizeof(O) * (len + 1)
    push!(fieldbuffers, Buffer(bufferoffset, blen))
    @debug 1 "made field buffer: bufferidx = $(length(fieldbuffers)), offset = $(fieldbuffers[end].offset), len = $(fieldbuffers[end].length), padded = $(padding(fieldbuffers[end].length, alignment))"
    bufferoffset += padding(blen, alignment)
    if eltype(A) == UInt8
        blen = length(col.data)
        push!(fieldbuffers, Buffer(bufferoffset, blen))
        @debug 1 "made field buffer: bufferidx = $(length(fieldbuffers)), offset = $(fieldbuffers[end].offset), len = $(fieldbuffers[end].length), padded = $(padding(fieldbuffers[end].length, alignment))"
        bufferoffset += padding(blen, alignment)
    else
        bufferoffset = makenodesbuffers!(col.data, fieldnodes, fieldbuffers, bufferoffset, alignment)
    end
    return bufferoffset
end

function writebuffer(io, col::Union{Map{T, O, A}, List{T, O, A}}, alignment) where {T, O, A}
    @debug 1 "writebuffer: col = $(typeof(col))"
    @debug 2 col
    writebitmap(io, col, alignment)
    # write offsets
    n = writearray(io, O, col.offsets.offsets)
    @debug 1 "writing array: col = $(typeof(col.offsets.offsets)), n = $n, padded = $(padding(n, alignment))"
    writezeros(io, paddinglength(n, alignment))
    # write values array
    if eltype(A) == UInt8
        n = writearray(io, UInt8, col.data)
        @debug 1 "writing array: col = $(typeof(col.data)), n = $n, padded = $(padding(n, alignment))"
        writezeros(io, paddinglength(n, alignment))
    else
        writebuffer(io, col.data, alignment)
    end
    return
end


text:
#=
HackerRank - Binomial distribution
https://www.hackerrank.com/challenges/s10-binomial-distribution-1/problem
https://www.hackerrank.com/challenges/s10-binomial-distribution-2/problem
=#
using Printf

function comb(n::Int, r::Int)
    r == 0 ? 1.0 : reduce(*, n-r+1:n, init=1) / reduce(*, 1:r, init=1)
end

function binomial(p::Float64, n::Int64, r::Int64)
    comb(n, r) * p ^ r * (1 - p) ^ (n - r)
end

# Binomial Problem 1
# What is the probability of having at least three boys out of six kids.
let ratio = 1.09    # Ratio of a boy or an event
    local p = ratio / (1.0 + ratio)
    # at least 3 out of 6
    local n = 6
    local c = 3
    local result = mapreduce(c -> binomial(p, n, c), +, c:n)
    @printf("%.3f", result)
end

## Binomial Problem 2
## What are the probabilities of:
##  No more than two rejected, and
##  At least two rejected.
let defect_ratio = 0.12, n = 10
    # No more than two rejected
    local result1 = mapreduce(x -> binomial(defect_ratio, n, x), +, 0:2)
    @printf("%.3f\n", result1)
    local result2 = 1.0 - mapreduce(x -> binomial(defect_ratio, n, x), +, 0:1)
    @printf("%.3f\n", result2)
end


text:
##########################################################################
# Helper functions for set partitions.

# A block is just a set of integers.
immutable Block
    elements::IntSet
end
Block() = Block(IntSet())

# Inherit part of the interface from IntSet
Base.show(io::IO, b::Block) = showcompact(io, b.elements)
Base.length(b::Block) = length(b.elements)
Base.push!(b::Block, e::Int) = push!(b.elements, e)

# The Partition type describes a partition of integers into blocks.
immutable Partition
    blocks::Set{Block}
end
Partition() = Partition(Set{Block}())

# Inherit part of the interface from Set
Base.show(io::IO, p::Partition) = showcompact(io, p.blocks)
Base.length(p::Partition) = length(p.blocks)

# Create partition from restricted growth string
function Partition(a::Vector; skip_check::Bool=false)
    I = length(a)
    a_max = maximum(a)
    num_blocks = 1 + a_max

    # Consistency checks
    if !skip_check
        all(sort(unique(a)) .== 0:a_max) || Lumberjack.error("Restricted growth string may not skip values.", @compat Dict( :a => a ))
    end

    # Build blocks and partition
    blocks = [ Block() for n = 1:num_blocks ]
    for l = 1:I
        push!(blocks[a[l] + 1], l)
    end

    return Partition(Set(blocks))
end

# Create partition from assignment matrix
Partition(A::Matrix; skip_check::Bool=false) =
    Partition(restricted_growth_string(A); skip_check=skip_check)

# Convert from partition to restricted growth string
restricted_growth_string(p::Partition) = restricted_growth_string(logical_matrix(p))

# Convert from partition to assignment matrix
function logical_matrix(p::Partition)
    I = maximum([ maximum(block.elements) for block in p.blocks ])
    A = eye(Int, I, I)

    for block in p.blocks
        elements = collect(block.elements)
        A[elements, elements] = 1
    end

    return A
end

# Convert from assignment matrix to restricted growth string
function restricted_growth_string(A::Matrix)
    I = size(A, 1)

    # Consistency checks
    A' == A || Lumberjack.error("Tried to create partition from asymmetric assignment matrix.", @compat Dict(:A => A))
    diag(A) == ones(I) || Lumberjack.error("Elements must belong to their own block.", @compat Dict(:A => A))

    # Build restricted growth string
    Ac = copy(A)
    a = zeros(Int, I)
    block_no = 0
    for i = 1:I
        # Find all integers in the block containing integer i.
        elements = find(Ac[i,:] .== 1)

        if length(elements) > 0
            # Consistency check
            for j in elements[2:end]
                find(Ac[j,:] .== 1) == elements || Lumberjack.error("Incorrect structure in assignment matrix.", @compat Dict(:A => Ac))
            end

            # Add all elements to block and increment number of blocks.
            a[elements] = block_no
            block_no += 1

            # Do not investigate these elements more
            Ac[elements,elements] = 0
        end
    end

    return a
end

# Convert from restricted growth string to assignment matrix
function logical_matrix(a::Vector)
    I = length(a)
    a_max = maximum(a)
    num_blocks = 1 + a_max

    # Consistency checks
    all(sort(unique(a)) .== 0:a_max) || Lumberjack.error("Restricted growth string may not skip values.", @compat Dict(:a => a))

    # Create matrix iteratively
    A = eye(Int, I, I)
    for n = 0:(num_blocks-1)
        elements = find(a .== n)
        A[elements, elements] = 1
    end

    return A
end

# PartitionIterator gives all partitions of the set 1:n. This iterator is
# based on Algorithm H from TAoCP 7.2.1.5. Note that a is zero-based.
immutable PartitionIterator
    n::Int

    PartitionIterator(n) = n >= 2 ? new(n) : error("only applicable for n >= 2.")
end

type PartitionIteratorState
    a::Vector{Int}
    b::Vector{Int}
    m::Int
    terminate::Bool
end
PartitionIteratorState(n) = PartitionIteratorState(zeros(n), ones(n - 1), 1, false)

Base.start(p::PartitionIterator) = PartitionIteratorState(p.n) # H1
Base.done(p::PartitionIterator, s) = s.terminate
function Base.next(p::PartitionIterator, s)
    # H2
    current_partition = Partition(s.a, skip_check=true) # By construction, a is a valid restricted growth string.

    if s.a[p.n] < s.m
        # H3
        s.a[p.n] += 1
    else
        # H4
        j = p.n - 1
        while s.a[j] == s.b[j]
            j -= 1
        end

        # H5
        if j == 1
            s.terminate = true
        else
            s.a[j] += 1

            # H6
            s.m = s.b[j] + (s.a[j] == s.b[j] ? 1 : 0)
            j += 1
            while j < p.n
                s.a[j] = 0
                s.b[j] = s.m
                j += 1
            end
            s.a[p.n] = 0
        end
    end

    # H2 cont'd.
    return current_partition, s
end

##########################################################################
# Stirling and Bell numbers

# Stirling2NumberCoefficientIterator gives all coefficients that when summed
# result in a Stirling number of the second kind. Note that the coefficients
# may be real numbers, but when summed are integers. The coefficients are
# thus calculated as floating point numbers.
immutable Stirling2NumberCoefficientIterator
    n::BigInt
    k::BigInt
    k_fac::BigInt
end
Stirling2NumberCoefficientIterator(n, k) = Stirling2NumberCoefficientIterator(n, k, factorial(k))
Base.start(s::Stirling2NumberCoefficientIterator) = BigInt(0) # state: j
Base.done(s::Stirling2NumberCoefficientIterator, j) = (j > s.k)
Base.next(s::Stirling2NumberCoefficientIterator, j) = ((iseven(s.k - j) ? 1. : -1.)*(binomial(s.k, j)*j^s.n/s.k_fac), j + 1)

# Stirling2NumberIterator gives all Stirling numbers of the second kind
# up to n. The Stirling numbers are integers, but the Stirling number
# coefficients are not. We therefore round the sum of the coefficients,
# in order to end up with an integer.
immutable Stirling2NumberIterator
    n::BigInt
end
Base.start(s::Stirling2NumberIterator) = BigInt(0) # state: k
Base.done(s::Stirling2NumberIterator, k) = (k > s.n)
Base.next(s::Stirling2NumberIterator, k) = (BigInt(round(sum(Stirling2NumberCoefficientIterator(s.n, k)))), k + 1)

# The Bell number can we calculated as the sum of all
# Stirling numbers of the second kind.
bell_number(n) = sum(Stirling2NumberIterator(n))


text:
# N.B. This file is maintained in TensorFlowBuilder source
export dtypesym

const DTYPE_SYMS = Symbol[:float32,
                          :float64,
                          :bfloat16,
                          :complex64,
                          :int8,
                          :uint8,
                          :int32,
                          :int64,
                          :bool,
                          :string,
                          :qint8,
                          :quint8,
                          :qint16,:quint16,:qint32,]
dtypesym(s::Symbol) = symbol("DT_$(uppercase(string(s)))")
dtypesym(o::PyObject) = dtypesym(symbol(o.x[:name]))

const DTYPES = Dict{PyObject,AbstractString}()

for dtype in DTYPE_SYMS
    sym = dtypesym(dtype)
    symsym = string(sym)
    if dtype in names(tf)
        @eval begin
            const $sym = Dtype(tf.$(dtype))
            export $sym
            DTYPES[$sym.x] = $symsym
        end
    end
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %%
using Plots
default(size=(400, 250))

Base.:^(f::Function, n::Integer) = n == 0 ? identity : foldr(∘, (f for _ in 1:n))

# %%
f(x) = sin(x)
plot(f, 0, 2; label="") |> display

plot([n^(1/2)*(f^n)(1.0) for n in 0:10:200]; label="")
hline!([√3]; label="")

# %%
f(x) = x - x^3/10
plot(f, 0, 2; label="") |> display

plot([n^(1/2)*(f^n)(1.0) for n in 0:10:200]; label="")

# %%
f(x) = x - x^2/4
plot(f, 0, 2; label="") |> display

plot([n*(f^n)(1.0) for n in 0:25:500]; label="")

# %%
f(x) = x - x^(3/2)/2
plot(f, 0, 2; label="") |> display

plot([n^2*(f^n)(1.0) for n in 0:25:500]; label="")

# %%

max_stars_repo_path: 0023/sin sin ... sin x.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0023/sin sin ... sin x.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0023/sin sin ... sin x.jl

max_forks_repo_name: genkuroki/public


text:
using Roots

"""
Return y'(x) using forward numerical differentiation.
"""
function num_differentiate(xs,ys)
    ders=zeros(length(xs)-1)
    for i in 1:length(xs)-1
        ders[i]=(ys[i+1]-ys[i])/(xs[i+1]-xs[i])
    end
    return ders
end

"""
Function where root x is relative extension of
worm-like chain in l>P regime with
    F = force
    P = persistence length

Bouchiat C, Wang MD, Allemand J, Strick T, Block SM, Croquette V.
Estimating the persistence length of a worm-like chain molecule from force-extension measurements.
Biophys J. 1999;76:409-413.
doi:10.1016/s0006-3495(99)77207-3
"""
function g_long(x, F, P, kT)
    a_two=-0.5164228
    a_three=-2.737418
    a_four=16.07497
    a_five=-38.87607
    a_six=39.49944
    a_seven=-14.17718

    return F*P/kT -(1-x)^(-2)/4. + 1. /4. -x - a_two*x^2-a_three*x^3-a_four*x^4-a_five*x^5-a_six*x^6-a_seven*x^7
end

"""
Function where root x is relative extension of
worm-like chain in other regime with l<P
    F = force
    P = persistence length
    kT = k_B T

Broedersz C.P., MacKintosh F.C.
Modeling semiflexible polymer networks
Rev. Mod. Phys. 2014;86:995-1036
doi:10.1103/RevModPhys.86.995
"""
function g_short(x, F, P, l, kT)
    ϕ=F*l^2/(pi^2*P*kT)
    return ϕ-9/pi^2*((1-x)^(-2)-1-x/3)
end

"""
Finds the relative extension of a chain with
    F = force
    P = persistence length
    kT = k_B T
    Marko_Siggia : if true, use Marko-Siggia, else
"""
function wlc_ext(L,P,F,kT=1, Marko_Siggia=true)
    if Marko_Siggia
        return L*find_zero(x->g_long(x,F,P,kT),(0,1),Bisection())
    else
        return L*find_zero(x->g_short(x,F,P,L,kT),(0,1),Bisection())
    end

end


text:
module SauterSchwab3D

using LinearAlgebra
using CompScienceMeshes
using StaticArrays
using FastGaussQuadrature
using GrundmannMoeller
using ShunnHamQuadrature

export sauterschwab_parameterized
export SauterSchwab3DStrategy,Singularity
export PositiveDistance6D, CommonVertex6D, CommonEdge6D, CommonFace6D, CommonVolume6D
export PositiveDistance6D_S, CommonVertex6D_S, CommonEdge6D_S, CommonFace6D_S, CommonVolume6D_S
export PositiveDistance5D, CommonVertex5D, CommonEdge5D, CommonFace5D
export PositiveDistance5D_S, CommonVertex5D_S, CommonEdge5D_S, CommonFace5D_S
export PositiveDistance4D, CommonVertex4D, CommonEdge4D, CommonFace4D
export PositiveDistance4D_S, CommonVertex4D_S, CommonEdge4D_S, CommonFace4D_S
export reorder,reverse_face_order, reverse_edge_order
#=
export CommonEdgeS, CommonEdgeS_a, CommonVertexS
export CommonVertexGM, Orig, Orig2, New, OrigEdge, OrigVertex, NewEdge, GM, GMEdge, GMVertex
=#
function sauterschwab_parameterized end

include("pulled_back_integrals.jl")
include("sauterschwabintegral.jl")
include("parametric_kernel_generator.jl")
include("reorder_vertices.jl")

end # module


text:
using ImageSmooth
using Test, ReferenceTests, TestImages
using FileIO
using ImageTransformations, ImageQualityIndexes
using ImageBase
using OffsetArrays

include("testutils.jl")

@testset "ImageSmooth.jl" begin
    include("algorithms/l0_smooth.jl")
end


text:
module StateVectorGates

    export SV_Gate, SV_OneQubitGate, SV_TwoQubitGate
    abstract type SV_Gate end
    abstract type SV_OneQubitGate <: SV_Gate end
    abstract type SV_TwoQubitGate <: SV_Gate end

    struct X <: SV_OneQubitGate
        _target::Int64
    end

    struct Y <: SV_OneQubitGate
        _target::Int64
    end

    struct Z <: SV_OneQubitGate
        _target::Int64
    end

    struct H <: SV_OneQubitGate
        _target::Int64
    end

    struct T <: SV_OneQubitGate
        _target::Int64
    end

    struct S <: SV_OneQubitGate
        _target::Int64
    end

    struct RX <: SV_OneQubitGate
        _target::Int64
        _theta::Float64
    end

    struct RY <: SV_OneQubitGate
        _target::Int64
        _theta::Float64
    end

    struct RZ <: SV_OneQubitGate
        _target::Int64
        _theta::Float64
    end

    struct U3 <: SV_OneQubitGate
        _target::Int64
        _theta::Float64
        _phi::Float64
        _lambd::Float64
    end

    struct CZ <: SV_TwoQubitGate
        _control::Int64
        _target::Int64
    end

    struct CX <: SV_TwoQubitGate
        _control::Int64
        _target::Int64
    end

    struct CP <: SV_TwoQubitGate
        _control::Int64
        _target::Int64
    end

    struct CRX <: SV_TwoQubitGate
        _control::Int64
        _target::Int64
        _theta::Float64
    end

    struct CRY <: SV_TwoQubitGate
        _control::Int64
        _target::Int64
        _theta::Float64
    end

    struct CRZ <: SV_TwoQubitGate
        _control::Int64
        _target::Int64
        _theta::Float64
    end

    mutable struct M <: SV_OneQubitGate
        _target::Int64
        _cregidx::Int64
        _result::Int64
    end

    # export all SV_Gate and TwoQubitGate
    for n in names(@__MODULE__; all=true)
        if Base.isidentifier(n) && n ∉ (Symbol(@__MODULE__), :eval, :include)
            if @eval typeof($n) <: DataType
                if @eval $n <: SV_Gate
                    @eval export $n
                end
            end
        end
    end
end

max_stars_repo_name: terasakisatoshi/blueqat-julia


text:
module scratchpad

#Use this to print
#=
using Weave
codeName = "TimeSeriesA6"
weave(Pkg.dir("$(pwd())\\$(codeName).jl"),
  informat="script",
  out_path = "$(pwd())\\output\\$(codeName)_Appendix.html",
  doctype = "md2html")
=#

#NOTE: To prepare the data, conform missing data to "". That means eliminating
#missing codes "N/A", "missing", and "--"

if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, CSV, StatsBase, GZip, ForwardDiff, Query,
  Gadfly, NLopt, HCubature, StaticArrays, Measures, Formatting

importall CTIO, CTReg, CTStat, CTVAR, CTNumerical

#############parameters
const VOL_PERIODS = 1

#path constants
const DATA_PATH = pwd() * "\\data" #path where we store the data
const OUTPUT_PATH = pwd() * "\\output"
const DATE_FORMAT = "m/d/yyyy"
const CHINA_DS_NAME = "ChinaDS1.1"
const CHINA_DS_NAME_L = "$(CHINA_DS_NAME)_L"
const CHINA_DS_NAME_W = "$(CHINA_DS_NAME)_W"
const R_TEST_PATH = "C:\\Users\\Clinton\\Dropbox\\Projects\\InterconnectRTester\\trainTester\\trainTester"
const ACCEPTABLE_RANGE = #use this to check if values are in acceptable ranges
    Dict(:value=>(-10., 10.), :marketCap=>(.1*10.^9.,Inf))
const DROP_OUTSIDE_RANGE = true
const MIN_SZSECOM_MARKETCAP = 6.*10.^9.
const KEEP_MISSING_MARKETCAP = false

const MIN_SAMPLES_IN_INDUSTRY = 5
#const MIN_SAMPLES_INTERACTED = 50 #obsolete
const MIN_COMPANY_QUARTERS = 1 #threshold for dropping companies
const SET_MISSING_IF_BELOW = true
const KEEP_LARGEST_WEIGHT_IF_DUP = true # keep the largest weight if multiple shareclasses

const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
const HEADER_NAME = "header.tex"
const DECIMALS = 2
const NUM_FORMAT = "{:.$(DECIMALS)f}"
#const PlotContainer = Union{Plot,Gadfly.Compose.Context,Vector{Plot}} #holds graphs
const OVERRIDE_STAR_STRINGS = ["\\ostar","\\dstar","\\tstar"]
const COLORS = ["DodgerBlue", "OrangeRed", "Green",
  "Brown", "Black", "BlueViolet"]

const PlotContainer = Union{Plot,Gadfly.Compose.Context,Vector{Plot}} #holds graphs

const DATE_SYMBOLS =  Symbol("12/31/2010"), Symbol("3/31/2011"), Symbol("6/30/2011"),
  Symbol("9/30/2011"), Symbol("12/31/2011"), Symbol("3/31/2012"), Symbol("6/30/2012"),
  Symbol("9/30/2012"), Symbol("12/31/2012"), Symbol("3/31/2013"), Symbol("6/30/2013"),
  Symbol("9/30/2013"), Symbol("12/31/2013"), Symbol("3/31/2014"), Symbol("6/30/2014"),
  Symbol("9/30/2014"), Symbol("12/31/2014"), Symbol("3/30/2015"), Symbol("6/30/2015"),
  Symbol("9/30/2015"), Symbol("12/31/2015"), Symbol("3/30/2016"), Symbol("6/30/2016"),
  Symbol("9/30/2016"), Symbol("12/31/2016"), Symbol("3/30/2017"), Symbol("6/30/2017"),
  Symbol("9/30/2017"), Symbol("12/31/2017"), Symbol("3/31/2018")

const FOCAL_DATES = Date("12/31/2010",  DATE_FORMAT):Dates.Month(3):Date("3/31/2018",  DATE_FORMAT)
const OUTCOMES = ["NORMALIZED_ACCRUALS_CF_METHOD"]

const MFloat64 = Union{Missing, Float64}
const MInt = Union{Missing, Int}
const MBool = Union{Missing, Bool}
const MSymbol = Union{Missing, Symbol}

n2s(x::T where T<:Real) = num2Str(x, 3, Ints=true)

function script(N::Int = 10^4)
X::Matrix{Float64}= [ones(N) rand(N) rand(N).*10000.]
Y::Vector{Float64} = X[:,2].*10 + rand(N).*5

println("First regression results:")
lm1::CTLM = CTLM(X[:,1:2],Y)
println("Coef: $(lm1.β)")
println("SE: $(diag(getWhiteΣ!(lm1)).^0.5)")

println("Noisy regression results:")
lm2::CTLM = CTLM(X,Y)
println("Coef: $(lm2.β)")
println("SE: $(diag(getWhiteΣ!(lm2)).^0.5)")

end

script()

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# Play each machine once.
# Then always play machine j that maximizes UCB value.
# Like sqrt(log(sum(counts).^2)) / sqrt(counts[j])

type UCB <: BanditAlgorithm
  counts::Vector{Int64}
  values::Vector{Float64}
end

function UCB(n_arms::Int64)
  UCB(zeros(Int64, n_arms), zeros(n_arms))
end

function initialize(algo::UCB, n_arms::Int64)
  algo.counts = zeros(Int64, n_arms)
  algo.values = zeros(n_arms)
end

function initialize(algo::UCB, arms::Vector)
  initialize(algo, length(arms))
end

function select_arm(algo::UCB)
  if any(algo.counts .== 0)
    return find(algo.counts .== 0)[1]
  else
    ucb_values = algo.values + sqrt((2 * log(sum(algo.counts))) ./ algo.counts)
    return ind_max(ucb_values)
  end
end

function update(algo::UCB, chosen_arm::Int64, reward::Real)
  algo.counts[chosen_arm] += 1
  n = algo.counts[chosen_arm]

  value = algo.values[chosen_arm]
  algo.values[chosen_arm] = ((n - 1) / n) * value + (1 / n) * reward
end

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


text:
using HDF5, JLD;

Qs = Array(Dict{Array{Uint8,1}, Dict{Uint8, Float64}}, 10);

JLD.save("qs.jld", "Qs", Qs);

Qs2 = JLD.load("qs.jld", "Qs");

max_stars_repo_name: tudor-berariu/phd

max_issues_repo_name: tudor-berariu/phd

max_forks_repo_name: tudor-berariu/phd


text:
module GaussianMarkovRandomFields

using LDLFactorizations
import LDLFactorizations: LDLFactorization
using LinearAlgebra, SparseArrays
using Distributions
using Random
using Memoize # For recursive marginal variance implementation

export GMRF, cholesky_ldl, logdet_ldl, prec

function cholesky_ldl(A)
    F = ldl(A)
    L = cholesky(F)
    return L, F.P
end

function LinearAlgebra.cholesky(F::LDLFactorizations.LDLFactorization)
    # https://en.wikipedia.org/wiki/Cholesky_decomposition#LDL_decomposition
	# Need next two lines for type-stability, the `getproperty` method for LDLFactorization
	# objects is not type-stable
	d = getfield(F, :d)
	L = SparseMatrixCSC(F.n, F.n, F.Lp, F.Li, F.Lx)
    all(d .> 0) || PosDefException(0)
    L = ((L + I) * Diagonal(sqrt.(d)))#[invperm(F.P), invperm(F.P)]
    return L
end

function logdet_from_chol(L)
    return 2 * sum(log(L[i]) for i in diagind(L))
end

function logdet_ldl(M::AbstractMatrix)
    L, P = cholesky_ldl(M)
    return logdet_from_chol(L)
end

struct GMRF{Tv<:AbstractVector{<:Real},
        Tm<:AbstractMatrix{<:Real},
        Tl<:AbstractMatrix{<:Real}} <: AbstractMvNormal
    μ::Tv
    Q::Tm
    L::Tl
end

function GMRF(μ, Q)#::AbstractVector, Q::AbstractMatrix)
    n = LinearAlgebra.checksquare(Q)
    length(μ) == n || DimensionMismatch("The dimensions of μ and Q are inconsistent.")
    L, P = cholesky_ldl(Q)
    return GMRF(μ, Q, L)
end

function GMRF(Q::Tv) where {Tv <: AbstractMatrix{<:Real}}
    n = LinearAlgebra.checksquare(Q)
    return GMRF(zeros(eltype(Q), n), Q)
end

function GMRF(μ::AbstractVector, F::LDLFactorization)
	P = F.P
	L = cholesky(F)
	Q = (L * L')[invperm(P), invperm(P)]
    return GMRF(μ, Q, F.L)
end
GMRF(F::LDLFactorization) = GMRF(zeros(size(F.L, 1)), F)

function Distributions._logpdf(d::GMRF, x::AbstractVector{T}) where T
    k = length(d.μ)
    ld = logdet_from_chol(d.L)
    x0 = x .- d.μ
    return -0.5 * (-ld + dot(x0, d.Q, x0) + k*log(2pi))
end

Base.length(d::GMRF) = length(mean(d))
Base.eltype(d::GMRF) = eltype(mean(d))
prec(d::GMRF) = d.Q
Distributions.mean(d::GMRF) = d.μ

# Uses the recursive method of Rue (2005)
# Finds the "future" non-zero indices
function ℐ(L::AbstractSparseMatrix, i::Integer)
	Lcol = L[:, i]
	ind, val = findnz(Lcol)
	ind[ind .> i]
end
# Recursively calculate the marginal variances. Should be memoized to avoid
# duplicated calculations and (probably more importantly) stack overflows in
# large precision matrices.
@memoize function cov_element(L::AbstractSparseMatrix, i::Integer, j::Integer)
	# Reverse indices to calculate upper triangle entry of Σ, otherwise get zero below
	if (i > j)
		 (j, i) = (i, j)
	end
	Σij = (i == j) ? 1 / L[i, i]^2 : zero(L[i, i])
	for k in ℐ(L, i)
        # `cov_element` will be zero if k > j; indices are switched if
        # necessary; this is where recursion occurs.
		Σij -= 1 / L[i, i] * L[k, i] * cov_element(L, k, j)
	end
	Σij
end
function Distributions.var(d::GMRF)
    n = length(d)
    v = Vector{Float64}(undef, n)
    # Iterate *backwards*
    for idx in reverse(eachindex(v))
        v[idx] = cov_element(d.L, idx, idx)
    end
    return v
end

Distributions.cov(d::GMRF) = error( "ERROR: The covariance matrix of a GMRF `g` can often be " *
    "dense and can cause the computer to run out of memory. If you are sure you have "*
    "enough memory, you can invert the precision matrix with `inv(Matrix(prec(g)))`")

function Distributions._rand!(rng::Random.AbstractRNG, d::GMRF, x::AbstractArray)
    z = randn(rng, size(x))
    x .= mean(d) .+ d.L \ z
end

end # module


text:
# This file was generated, do not modify it. # hide
raw_data = urldownload("https://github.com/tlienart/DataScienceTutorialsData.jl/blob/master/data/wri_global_power_plant_db_be_022020.csv?raw=true")
data = DataFrame(raw_data);


text:
using LowRankModels

# test losses in losses.jl
srand(1);

losses = [
QuadLoss(),
QuadLoss(10),
L1Loss(),
L1Loss(5.2),
HuberLoss(),
HuberLoss(4),
HuberLoss(3.1, crossover=3.2),
PeriodicLoss(2*pi),
PeriodicLoss(2*pi, 4),
PoissonLoss(20),
PoissonLoss(22,4.1),
OrdinalHingeLoss(1,10),
OrdinalHingeLoss(2,7,5),
LogisticLoss(),
LogisticLoss(0.2),
WeightedHingeLoss(),
WeightedHingeLoss(11),
WeightedHingeLoss(1.5, case_weight_ratio=4.3),
MultinomialLoss(4),
MultinomialLoss(6, .5),
# OrdisticLoss(5),
MultinomialOrdinalLoss(3)
] #tests what should be successful constructions

# TODO: do some bad constructions and test that they fail with catches
bad_losses = [
:(QuadLoss(10,RealDomain)),
:(HuberLoss(3.1, 3.2)),
:(PeriodicLoss(scale=2*pi)),
:(PeriodicLoss(2*pi, scale=4)),
:(PeriodicLoss())
]
for expression in bad_losses
	try
		eval(expression);
		println("test FAILED for $expression")
	catch
		println("test PASSED for $expression (failed to construct)")
	end
end

m,n,k = 1000, length(losses), 5;
d = embedding_dim(losses)
X_real, Y_real = 2*randn(m,k), 2*randn(k,d);
XY_real = X_real*Y_real;

# tests default imputations and implicit domains
# we can visually inspect the differences between A and A_real to make sure imputation is right
A = impute(losses, XY_real)

regscale = 1
yregs = Array(Regularizer, length(losses))
for i=1:length(losses)
	if typeof(losses[i]) == MultinomialOrdinalLoss ||
		typeof(losses[i]) == OrdisticLoss
		yregs[i] = OrdinalReg(QuadReg(regscale))
	else
		yregs[i] = QuadReg(regscale)
	end
end

# tests all the M-estimators with scale=false, offset=false
glrm = GLRM(A, losses, QuadReg(regscale), yregs, 5, scale=false, offset=false);

# interestingly adding an offset to a model with multidimensional ordinal data causes a segfault
# but let's test the offset for everything but ordinals
# oops we still get a segfault...
# tamecols = [typeof(losses[i]) !== MultinomialOrdinalLoss &&
# 			typeof(losses[i]) !== OrdisticLoss
# 			for i=1:length(losses)]
# glrm = GLRM(A[:, tamecols],
# 	losses[tamecols],
# 	QuadReg(regscale),
# 	yregs[tamecols],
# 	5, scale=false, offset=true)

# tests eval and grad
@time X,Y,ch = fit!(glrm);

# tests initialization
init_svd!(glrm)


text:
function set_mp_prop!(g::MetaDiGraph{T, U}, mp::MainPathResult, sym::Symbol) where T <: Integer where U <: Real
        [set_prop!(g, e, sym, false) for e in edges(g)]
        [set_prop!(g, e, sym, true) for e in mp.edges]
        vmp = unique([src.(mp.edges); dst.(mp.edges)])
        [set_prop!(g, v, sym, false) for v in vertices(g)]
        [set_prop!(g, v, sym, true) for v in vmp]
end

# function add_weights(g::SimpleDiGraph{T}, weights::Vector{U}) where T <: Integer where U <: Real
#     from = [src(e) for e in edges(g)]
#     to = [dst(e) for e in edges(g)]
#     SimpleWeightedDiGraph(from, to, weights)
# end


text:
abstract type AbstractParticleSwarm <: MLJTuning.TuningStrategy end

struct ParticleSwarmState{T, R, P, I}
    ranges::R
    parameters::P
    indices::I
    X::Matrix{T}
    V::Matrix{T}
    pbest_X::Matrix{T}
    gbest_X::Matrix{T}
    pbest::Vector{T}
    gbest::Vector{T}
end


text:
using Random
using Distributions
using Statistics
using RoboLib.Geom: Pose2D
#TODO(cx) extend or pass fn?
import mushr_pf.MotionModels: reset!
import mushr_pf.MotionModels: predict!
import mushr_pf.SensorModels: reweight!
using LinearAlgebra
using StatsBase: Weights, sample!

# TODO(cxs): port what works out back to upstream
struct BufferedWeightedParticleBelief{T,W}
    particles::Vector{T}
    _particle_memory::Vector{T}

    weights::Vector{W}
    _weight_memory::Vector{W}

end

@inline resetweights!(b::BufferedWeightedParticleBelief{T,W}) where {T,W} = fill!(b.weights, 1/length(b.weights))

@inline function BufferedWeightedParticleBelief(particles::AbstractVector{T}, weights::AbstractVector{W}) where {T,W}
    BufferedWeightedParticleBelief{T,W}(particles, similar(particles), weights, similar(weights))
end
@inline function BufferedWeightedParticleBelief{T,W}(n::Integer) where {T,W}
    BufferedWeightedParticleBelief(Vector{T}(undef, n), Vector{W}(undef, n))
end
@inline BufferedWeightedParticleBelief{T}(n::Integer) where {T} = BufferedWeightedParticleBelief{T, Float64}(n)
@inline function BufferedWeightedParticleBelief{W}(particles::AbstractVector) where {W}
    weights = fill(one(W)/length(particles), length(particles))
    BufferedWeightedParticleBelief(particles, weights)
end
@inline BufferedWeightedParticleBelief(particles::AbstractVector) = BufferedWeightedParticleBelief{Float64}(particles)

Distributions.mean(b::BufferedWeightedParticleBelief) = sum((b.weights .* b.particles)) ./ sum(b.weights)
Distributions.mode(b::BufferedWeightedParticleBelief) = b.particles[argmax(b.weights)]
Distributions.mean(b::BufferedWeightedParticleBelief{<:Pose2D}) = mean(b.particles, b.weights)

struct StatefulParticleFilter{PC,P,R,S,RNG<:AbstractRNG}
    belief::PC
    predictmodel::P
    reweightmodel::R
    resample!::S
    rng::RNG
end

#@inline function StatefulParticleFilter(pc::BufferedWeightedParticleBelief, predictmodel, reweightmodel)
#    StatefulParticleFilter(pc, predictmodel, reweightmodel, Random.GLOBAL_RNG)
#end
#@inline function StatefulParticleFilter(n::Integer, predictmodel, reweightmodel)
#    StatefulParticleFilter(BufferedWeightedParticleBelief
#end

#todo: handle no reweight
function resample_lowvar!(b::BufferedWeightedParticleBelief, rng::AbstractRNG)
    error("fix me")
    ws = sum(b.weights)
    n = length(b.weights)
    pm = b._particle_memory
    wm = b._weight_memory

    dr = 1/n*ws
    r = rand(rng) * dr
    c = b.weights[1]
    U = r
    i = 1
    @inbounds for m in 1:n
        while U > c
            i += 1
            c += b.weights[i]
        end
        U += dr
        pm[m] = b.particles[i]
        wm[m] = b.weights[i]
    end
    copy!(b.particles, pm)
    #NOTE:(cxs): were weights not reset before?)
    #before: copy!(b.weights, wm)
    resetweights!(b)
end

function resample_naive!(b::BufferedWeightedParticleBelief, rng::AbstractRNG)
    weights = Weights(b.weights)
    @inbounds sample!(b.particles, weights, b.particles)
    resetweights!(b)
end

@inline function update_act!(up::StatefulParticleFilter, a, dt)
    predict!(up.belief.particles, a, up.rng, dt, up.predictmodel)
    up
end

@inline function update_obs!(up::StatefulParticleFilter, o, dt)
    b = up.belief
    reweight!(b.weights, b.particles, o, up.rng, dt, up.reweightmodel)
    up
end

@inline function reset!(pf::StatefulParticleFilter, state)
    reset!(pf.belief.particles, state, pf.rng, pf.predictmodel)
    resetweights!(pf.belief)
    pf
end

@inline resample!(up::StatefulParticleFilter) = up.resample!(up.belief, up.rng)

#function update!(up::StatefulParticleFilter, a, o, dt)
#    b = up.belief
#    pm = b._particle_memory
#    p = b.particles
#    wm = b._weight_memory
#    w = b.weights
#    resize!(pm, length(b.particles))
#    resize!(wm, length(b.particles))
#    predictmodel(p, up.predict_model, a, up.rng, dt)
#    reweightmodel(w, p, up.reweight_model, o, up.rng, dt)
#    # TODO(cxs): inplace resampler
#    resample!(b, up.rng)
#    return up
#end


text:
using GrowthMaps, GeoData, ModelParameters, Unitful, Test, Plots
using Unitful: °C, K, hr, d, mol, cal
using GrowthMaps: rate, condition, conditionalrate

#= This file just tests that the fitting routines run, not
that they work, which is difficult with Interact.jl. =#

p = Param(0.3; units=nothing, bounds=(0.0, 1.0))
ΔH_A = Param(2e4; units=u"cal*mol^-1", bounds=(2e3, 2e5))
ΔH_L = Param(-1e5; units=u"cal*mol^-1", bounds=(1e-6, 1e-4))
ΔH_H = Param(3e5; units=u"cal*mol^-1", bounds=(3e4, 3e6))
T_halfL = Param(250.0; units=u"K", bounds=(200, 300))
T_halfH = Param(300.0; units=u"K", bounds=(200, 300))
T_ref = K(25.0°C)
growth = Layer(:temp, K, SchoolfieldIntrinsicGrowth(p, ΔH_A, ΔH_L, T_halfL, ΔH_H, T_halfH, T_ref))
model = Model(growth)

obs = [(10.0°C, 0.0), (20.0°C, 0.2), (30.0°C, 0.1), (40.0°C, 0.0)]

@testset "Auto fit" begin
    fitted = fit!(model, obs)
end

@testset "Manual fit interface" begin
    data = (temp=0.0:1.0:45.0,)

    withunits(model, :bounds)
    interface = manualfit!(model, data; obs=obs);
end

@testset "Map fit interface" begin
    ## Set up series data
    tempdata = GeoArray.(
        [rand(0.0:40.0, 10, 10) for i in 1:5], Ref((X, Y));
        name=:stress, missingval=-99.0
    )
    stacks = [GeoStack((temp=tempdata[i],)) for i in 1:length(tempdata)]
    timedim = (Ti((1:1:5)hr; mode=Sampled(; span=Regular(1hr))),)
    modelkwargs = (series=GeoSeries(stacks, timedim), tspan=1hr:1hr:5hr)
    interface = mapfit!(model, modelkwargs; occurrence=[(1, 2), (9, 10)]);
end

# To test interfaces manually
# using Blink
# w = Blink.Window()
# body!(w, interface)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "sais"
version = v"2.4.1"

# Collection of sources required to complete build
sources = [
    "https://sites.google.com/site/yuta256/sais-2.4.1.zip" =>
    "467b7b0b6ec025535c25e72174d3cc7e29795643e19a3f8a18af9ff28eca034a",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/sais-2.4.1/
if [[ $nbits == 64 ]]; then
  # modify the build to put sais and sais64 into the same library
  build64=ON
  sed -ri \
    -e "s/sais sais.c/sais sais.c sais64.c/" \
    -e "s/sais PROPERTIES/sais PROPERTIES\n  COMPILE_FLAGS \"-DSAIS_BUILD_64BIT\"/" \
    lib/CMakeLists.txt
else
  build64=OFF
fi
mkdir build && cd build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=$CMAKE_TARGET_TOOLCHAIN -DCMAKE_BUILD_TYPE=Release -DBUILD_SHARED_LIBS=ON -DBUILD_EXAMPLES=OFF -DBUILD_SAIS64=$build64 ..
make -j$nproc
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libsais", :libsais)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module StanMCMCChain

# package code goes here

using Reexport

@reexport using CmdStan, MCMCChain, Plots, Statistics, JLD
import CmdStan: convert_a3d

include("utilities/convert_a3d.jl")

export
  convert_a3d

end # module


text:
#Annotators and their dependencies.
#See https://stanfordnlp.github.io/CoreNLP/dependencies.html

#Value is an unordered list of *all* annotators required by the Key.
ANNOTATOR_DEPENDENCIES = Dict(
    "tokenize" => String[],
    "quote" => String[],
    "regexner" => String[],
    "cleanxml" => ["tokenize"],
    "ssplit" => ["tokenize"],
    "pos" => ["tokenize", "ssplit"],
    "lemma" => ["tokenize", "ssplit", "pos"],
    "ner" => ["tokenize", "ssplit", "pos", "lemma"],
    "parse" => ["tokenize", "ssplit", "pos"],
    "depparse" => ["tokenize", "ssplit", "pos"],
    "sentiment" => ["tokenize", "ssplit", "pos", "parse"],
    "natlog" => ["tokenize", "ssplit", "pos", "lemma", "depparse"],
    "dcoref" => ["tokenize", "ssplit", "pos", "lemma", "ner", "parse"],
    "relation" => ["tokenize", "ssplit", "pos", "lemma", "ner", "parse", "depparse"],
    "mention" => ["tokenize", "ssplit", "pos", "parse", "lemma", "ner", "depparse"],
    "coref" => ["tokenize", "ssplit", "pos", "parse", "lemma", "ner", "depparse", "mention"],
    "openie" => ["tokenize", "ssplit", "pos", "lemma", "depparse", "natlog"]
)

#Each annotator is to the right of its dependencies.
ORDERED_ANNOTATOR_LIST = ["tokenize", "regexner", "quote", "ssplit",
        "cleanxml", "pos", "lemma", "parse", "depparse", "sentiment", "ner",
        "natlog", "dcoref", "relation", "mention", "coref", "openie"]


text:
module TimeToContact

using Plots
using Images
using FastConv
using ImageView
import VideoIO

include("./utils.jl");

end # module


text:
function dissociation_calculation(df::DataFrame)
    result = DataFrame(delta_t = Float64[], dissociation_rate = Float64[])
    tmp_df = tmp_df = df[df.frame_num.==2, :]
    for i = 1:maximum(df.lifetime)
        delta_t = i * (1 / 45)
        rate = nrow(tmp_df[tmp_df.lifetime.==i, :])
        push!(result, [delta_t, rate])
    end
    result[:, :rate] ./= maximum(result[:, :dissociation_rate])
end

# function fit_dissociation(
#     df::DataFrame
#     )
#     model(x, p)
# end

# function fit_msd(
#     df::DataFrame;
#     max_time::Int64 = 10,
#     loc_error::Float64 = 0.03,
#     p0::AbstractVector = [1.0, 1.0]
#     )
#     @. model(x, p) = 4 * p[1] * x^p[2] + 4 * 0.03^2
#     fit = curve_fit(model, df.delta_t[1:max_time], df.msd[1:max_time], p0)
#     fit.param
# end

max_stars_repo_name: SosUts/spt

max_issues_repo_name: SosUts/spt

max_forks_repo_name: SosUts/spt


text:
function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
    TS = MOI.OPTIMIZE_NOT_CALLED
    ts = status(optimizer)
    if ts == :iteration_limit
        TS = MOI.ITERATION_LIMIT
    elseif ts == :time_limit
        TS = MOI.TIME_LIMIT
    elseif ts == :solution_limit
        TS = MOI.SOLUTION_LIMIT
    end
    return TS
end

function MOI.get(optimizer::Optimizer, ::MOI.PrimalStatus)
    return has_solution(optimizer) ? MOI.FEASIBLE_POINT : MOI.NO_SOLUTION
end

MOI.get(::Optimizer, ::MOI.DualStatus) = MOI.NO_SOLUTION

function MOI.get(optimizer::Optimizer, ::MOI.ObjectiveValue)
    return is_sat(optimizer) ? 0 : best_value(optimizer)
end

MOI.get(optimizer::Optimizer, ::MOI.ObjectiveBound) = _best_bound(optimizer)
MOI.get(optimizer::Optimizer, ::MOI.ResultCount) = 1

function MOI.get(optimizer::Optimizer, ::MOI.VariablePrimal, vi::MOI.VariableIndex)
    if has_solution(optimizer)
        return best_values(optimizer)[vi.value]
    else
        return get_value(optimizer, vi.value)
    end
end

MOI.get(optimizer::Optimizer, ::MOI.SolveTime) = time_info(optimizer)[:total_run]

function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    return has_solution(optimizer) ? "Satisfying solution" : "No solutions"
end


text:
using GeometricFlux
using Flux
using Zygote
using CUDA

Flux.@functor GeometricFlux.FeaturedGraph

include("cuda/cuda.jl")
include("representation/representation.jl")
include("nn_structures/nn_structures.jl")
include("utils.jl")


text:
export create_curl

# Wrapper to create the discrete curl by default
create_curl(isfwd::AbsVecBool,  # isfwd[w] = true|false: ∂w is forward|backward difference
            N::AbsVecInteger,  # size of grid
            ∆l⁻¹::NTuple{K,Number}=ntuple(k->1.0,length(isfwd)),  # ∆l⁻¹[w]: inverse of uniform distance between grid planes in x-direction
            isbloch::AbsVecBool=fill(true,K),  # boundary conditions in x, y, z
            e⁻ⁱᵏᴸ::AbsVecNumber=ones(K);  # Bloch phase factor in x, y, z
            cmp_shp::AbsVecInteger=1:K,
            cmp_out::AbsVecInteger=1:K,
            cmp_in::AbsVecInteger=1:K,
            order_cmpfirst::Bool=true  # true to use Cartesian-component-major ordering for more tightly banded matrix
            ) where {K} =
    # I should not cast e⁻ⁱᵏᴸ into a complex vector, because then the entire curl matrix
    # becomes a complex matrix.  Sometimes I want to keep it real (e.g., when no PML and
    # Bloch phase factors are used).  In fact, this is the reason why I accept e⁻ⁱᵏᴸ instead
    # of constructing it from k and L as exp.(-im .* k .* L), which is always complex even
    # if k = 0.
    #
    # I should not cast ∆l⁻¹ to a vector of any specific type (e.g., Float, ComplexF), either,
    # because sometimes I would want to even create an integral curl operator.
    create_curl(isfwd, fill.(∆l⁻¹,(N...,)), isbloch, e⁻ⁱᵏᴸ; cmp_shp, cmp_out, cmp_in, order_cmpfirst)

# Wrapper to convert AbstractVector's to SVec's
create_curl(isfwd::AbsVecBool,  # isfwd[w] = true|false: ∂w is forward|backward difference
            ∆l⁻¹::NTuple{K,AbsVecNumber},  # ∆l⁻¹[w]: inverse of uniform distance between grid planes in x-direction
            isbloch::AbsVecBool=fill(true,K),  # boundary conditions in x, y, z
            e⁻ⁱᵏᴸ::AbsVecNumber=ones(K);  # Bloch phase factor in x, y, z
            cmp_shp::AbsVecInteger=1:K,
            cmp_out::AbsVecInteger=1:K,
            cmp_in::AbsVecInteger=1:K,
            order_cmpfirst::Bool=true  # true to use Cartesian-component-major ordering for more tightly banded matrix
            ) where {K} =
    (Kout = length(cmp_out); Kin = length(cmp_in);
     create_curl(SVec{K}(isfwd), ∆l⁻¹, SVec{K}(isbloch), SVec{K}(e⁻ⁱᵏᴸ);
                 cmp_shp=SInt{K}(cmp_shp), cmp_out=SInt{Kout}(cmp_out), cmp_in=SInt{Kin}(cmp_in), order_cmpfirst))

# ∇×F = x̂ (∂y Fz - ∂z Fy) + ŷ (∂z Fx - ∂x Fz) + ẑ (∂x Fy - ∂y Fx)
#
# In matrix form,
# ⎡ 0  -∂z  ∂y⎤
# ⎢ ∂z  0  -∂x⎥
# ⎣-∂y  ∂x  0 ⎦
#
# CURL_BLK takes only the signs of the blocks:
const CURL_BLK = S²Int{3,9}(0,  1, -1,  # 1st column (not row)
                           -1,  0,  1,  # 2nd column (not row)
                            1, -1,  0)  # 3rd column (not row)

function create_curl(isfwd::SBool{K},  # isfwd[w] = true|false: ∂w is forward|backward difference
                     ∆l⁻¹::NTuple{K,AbsVecNumber},  # ∆l⁻¹[w]: inverse of distances between grid planes in w-direction
                     isbloch::SBool{K},  # isblock[w]: boundary conditions in w-direction
                     e⁻ⁱᵏᴸ::SNumber{K};  # e⁻ⁱᵏᴸ[w]: Bloch phase factor in w-direction
                     cmp_shp::SInt{K}=SVec(ntuple(identity, Val(K))),
                     cmp_out::SInt{Kout}=SVec(ntuple(identity, Val(K))),
                     cmp_in::SInt{Kin}=SVec(ntuple(identity, Val(K))),
                     order_cmpfirst::Bool=true  # true to use Cartesian-component-major ordering for more tightly banded matrix
                     ) where {K,Kout,Kin}  # {space dimension, output field dimension, input field dimension}
    T = promote_type(eltype.(∆l⁻¹)..., eltype(e⁻ⁱᵏᴸ))  # eltype(eltype(∆l⁻¹)) can be Any if ∆l⁻¹ is inhomogeneous
    N = SVec(length.(∆l⁻¹))
    M = prod(N)

    num_blk = sum(abs, CURL_BLK[cmp_out,cmp_in])  # no allocations

    # Below, create_∂info() is called num_blk times in the double for loops, and each call
    # returns I, J, V of length 2M, so we preallocate Itot, Jtot, Vtot of length num_blk⋅2M.
    Itot = VecInt(undef, num_blk * 2M)
    Jtot = VecInt(undef, num_blk * 2M)
    Vtot = Vector{T}(undef, num_blk * 2M)

    ind_blk = 0  # index of matrix block
    for ind_nv = 1:Kout
        nv = cmp_out[ind_nv]  # Cartesian compotent of output field
        istr, ioff = order_cmpfirst ? (Kout, ind_nv-Kout) : (1, M*(ind_nv-1))  # (row stride, row offset)
        for ind_nu = 1:Kin
            nu = cmp_in[ind_nu]  # Cantesian component of input field

            parity = CURL_BLK[nv,nu]
            iszero(parity) && continue  # if zero block, skip to next nu

            jstr, joff = order_cmpfirst ? (Kin, ind_nu-Kin) : (1, M*(ind_nu-1))  # (column stride, column offset)

            nw = 6 - nv - nu  # direction of differentiation; 6 = 1 + 2 + 3
            is_nw = cmp_shp.==nw
            sum(is_nw)==1 || @error "cmp_shp = $cmp_shp does not have one and only one nw = $nw"
            ind_nw = findfirst(is_nw)
            I, J, V = create_∂info(ind_nw, isfwd[ind_nw], N, ∆l⁻¹[ind_nw], isbloch[ind_nw], e⁻ⁱᵏᴸ[ind_nw])

            @. I = istr * I + ioff
            @. J = jstr * J + joff
            V .*= parity

            # For some reason, using .= below is slower because it uses 1 allocatiotn.  On
            # the other hand, using = does not use allocation and therefore faster.
            indₛ, indₑ = ind_blk*2M + 1, (ind_blk+1)*2M
            Itot[indₛ:indₑ] = I
            Jtot[indₛ:indₑ] = J
            Vtot[indₛ:indₑ] = V
            ind_blk += 1
        end
    end

    return dropzeros!(sparse(Itot, Jtot, Vtot, Kout*M, Kin*M))
end


text:
function convert_to_image(var::AbstractVector, grd::AbstractGrid)
    x = Matrix{Float64}(undef, reverse(cells(grd))...)
    fill!(x, NaN)
    xind, yind =  indices(grd, 1), indices(grd,2) #since matrices are drawn from upper left corner
    [x[yind[i], xind[i]] = val for (i, val) in enumerate(var)]
    x
end

RecipesBase.@recipe function f(var::AbstractVector, grd::AbstractGrid)
    seriestype := :heatmap
    aspect_ratio --> :equal
    grid --> false
    xrange(grd), yrange(grd), convert_to_image(var, grd)
end

# RecipesBase.@recipe function f(sit::SiteFields) # not sure what SiteFields are
#     ones(nsites(sit)), sit
# end

RecipesBase.@recipe function f(var::AbstractVector, pnt::AbstractPoints)
    seriestype := :scatter
    aspect_ratio --> :equal
    grid --> false
    marker_z := var
    legend --> false
    colorbar --> true
    cd = coordinates(pnt)
    cd[:,1], cd[:,2]
end

RecipesBase.@recipe function f(asm::AbstractAssemblage; showempty = false)
    var = richness(asm)
    if !showempty
        var = [Float64(v) for v in var]
        (var[var.==0] .= NaN)
    end
    var, getcoords(places(asm))
end

RecipesBase.@recipe function f(var::AbstractVector, asm::AbstractAssemblage)
    var, getcoords(places(asm))
end

RecipesBase.@recipe function f(g::Function, asm::AbstractAssemblage)
    g(asm), getcoords(places(asm))
end


text:
import Base: show

show(io::IO, ::MIME"text/html", p::MExpr{:Graphics}) =
  print(io, ExportString(p, "SVG"))

show(io::IO, ::MIME"image/png", p::MExpr{:Graphics3D}) =
  print(io, @math ExportString($p, "PNG"))


text:
export make_gui, add_callbacks, add_spikes, add_video, add_ttl_cov, add_labels, add_times, add_spike_ts

function make_gui(mypath)

    window = glscreen()

    description = """
    Usage of the play slider widget
    """
    #Create screens

    editarea, viewarea = y_partition(window.area, 20)
    edit_screen = Screen(window, area = editarea)
    viewscreen = Screen(window, name = :viewscreen,area = viewarea)

    dataarea, imgarea = x_partition(viewarea,50)
    imgscreen = Screen(viewscreen, area = imgarea)
    datascreen = Screen(viewscreen, area = dataarea)

    slider_area, time_area = x_partition(editarea,50)
    slider_screen = Screen(edit_screen, area=slider_area)
    time_screen = Screen(edit_screen, area=time_area)

    vid_path = string(mypath,filter(x->contains(x,".mp4"),readdir(mypath))[1])

    max_time=10000

    #Y data is the covariate of interest
    y_data = rand(1f0:100f0,max_time,3)

    detect=DetectNeg();
    cluster=ClusterTemplate();
    align=AlignMin();
    feature=FeatureTime();
    reduce=ReductionNone();
    thres=ThresholdMeanN();
    num_channels=1;

    s=create_multi(detect,cluster,align,feature,reduce,thres,1,48,Float32);
    buf=Spike[Spike() for i=1:10000,j=1:1];
    nums=zeros(Int64,1)
    s[1].thres=-1.0

    tif = Array{ColorTypes.Gray{FixedPointNumbers.Normed{UInt8,8}},3}(0,0,0)
    whiskers=[Whisker() for i=1:1, j=1:1]

    analysis_gui(window,edit_screen,viewscreen,imgscreen,datascreen,time_screen, slider_screen,
    0.2f0,
    Point2f0[Point2f0(i*5,0.0)  for i=1:500,j=1:3],
    vid_path,mypath,y_data,max_time,ones(Float32,3),zeros(Int64,max_time),0,100,
    1,zeros(Int64,0),s[1],buf,nums,false,[Point2f0(0.0,0.0) for i=1:50,j=1:4],
    Reactive.Signal(1),[zeros(Int64,0) for i=0:0],zeros(Int64,length(-3000:300:3000)-1),1:10:10000,tif,whiskers,false,[0])
end

function add_spikes(gui,channel_num,data_type="Intan")

    if data_type == "Intan"
        spike_path = string(gui.folder_path,"v.bin")
        spikes = Intan.parse_v(spike_path)
    else

        spike_path = get_spike_path(gui.folder_path,channel_num)
        io_spike = open(spike_path,"r");

        spikes = SampleArray(Float32,io_spike);
        spikes = Array(spikes)
        close(io_spike)
    end

    gui.y_data=zeros(Float32,length(spikes),3)
    gui.y_data[:,1]=spikes

    nothing
end

function add_times(gui,channel_num,data_type="Intan")

    if data_type == "Intan"
        spike_path = string(gui.folder_path,"v.bin")
        gui.start_time=1
        gui.max_time=size(Intan.parse_v(spike_path),1)
    else
        spike_path = get_spike_path(gui.folder_path,channel_num)
        io_spike = open(spike_path,"r")
        times=TimeArray(Int64,io_spike);

        gui.start_time=times[1]
        gui.max_time=length(times)
        close(io_spike)
    end

    slider_step = round(Int64,gui.max_time / ((gui.max_time-100) / 1000 * 10))
    gui.slider_values = 30000:slider_step:gui.max_time

    nothing
end

function add_video(gui,channel_num,data_type="Intan")

    if data_type == "Intan"
        myttl = Intan.parse_ttl(string(gui.folder_path,"ttl.bin"))
        my_vid_times=myttl[channel_num][find(diff(myttl[channel_num]).>1)]
        my_vid_array = zeros(Int64,gui.max_time)
        for i=2:length(my_vid_times)
            my_vid_array[my_vid_times[i-1]:my_vid_times[i]]=i
        end
        my_vid_array[my_vid_times[end]:end]=length(my_vid_times)

        gui.video_ts=my_vid_array
    else

        event_path = string(gui.folder_path, "all_channels.events")
        io_event = open(event_path,"r");

        xx=parse_ttl(io_event,channel_num,gui.start_time,gui.max_time)
        close(io_event)
        gui.video_ts=xx[3]

    end

    #If i seek to the last frame, mpv crashes, so set to second to last frame
    gui.video_ts[gui.video_ts.==gui.video_ts[end]]=gui.video_ts[end]-1;

    stdout, stdin, process = mpv_open(gui.vid_path)
    sleep(1.0)
    run(pause_cmd)

    nothing
end

#If there are not any events, should probably state that in a message and not error
function add_ttl_cov(gui,channel_num,cov_num,data_type="Intan")

    if data_type == "Intan"
        myttl = Intan.parse_ttl(string(gui.folder_path,"ttl.bin"))
        gui.y_data[myttl[channel_num],cov_num] = 1
        push!(gui.event_ts,myttl[channel_num][find(diff(myttl[channel_num]).>1)])
    else

        event_path = string(gui.folder_path, "all_channels.events")
        io_event = open(event_path,"r");

        xx=parse_ttl(io_event,channel_num,gui.start_time,gui.max_time)

        gui.y_data[:,cov_num] = xx[2]

        push!(gui.event_ts,xx[1]-gui.start_time)

        close(io_event)
    end

    nothing
end

function add_spike_ts(gui,filename)

    gui.pre_sorted=true

    myspikes=matread(string(gui.folder_path,"v_sorted.mat"))

    times=myspikes["Channel01"][:,1];

    gui.pre_sorted_ts=round(Int64,times.*30000)

    nothing
end

function get_spike_path(folder_path,channel_num)

    adc_channels=filter(x->contains(x,".continuous"),readdir(folder_path))

    spike_path=string(folder_path,adc_channels[channel_num])

    spike_path
end

function add_callbacks(gui)

    #slider is used to start and stop animation, as well as drag to certain point
    #Play slider renders at 30 fps
    #Our DAQ records at 30000, so we can roughly control our speed of playback by how far
    #advance in time with each frame in that 30 fps
    #For instance, if we play just 1 data point per frame, that would be 30000 / 30 = 1/1000 of real time
    #Lets slow down to 1/10 of normal

    iconsize = 8mm
    play_viz, slider_value = play_slider(
        gui.sliderscreen, iconsize, gui.slider_values
    )
    #Each slider value change sends a signal to covariate plotter
    #100 points are plotted on the x axis.
    #The time scale dictates how far in the future we can visualize
    #For instance, 100 / 30000 = is only 3 ms
    #We want to see the previous points as time advances, so these should always be some multiple of the slowdown
    time_scale=1;
    my_animation = map(slider_value) do t

        video_frame_rate = 25 #We should calculate this
        #Because the video is not continuous, at T we should
        #find the total frame count at that index and use that instead

            frame_num = gui.video_ts[round(Int,t)]
            gui.t = round(Int,t)
            if gui.current_frame != frame_num
                gui.current_frame = frame_num
                run(myseek(frame_num / video_frame_rate))
            end

            plot_lines(gui,t,increment=2,p_per_frame=100)

        gui.cov1
    end

    my_spikes = map(slider_value) do t

        if gui.show_spikes
            if gui.pre_sorted ==false
                sort_spikes(gui,t)
            else
                get_sorted_ts(gui,t)
            end
            if length(gui.event_ts)>0
                event_triggered(gui,t)
            end
            plot_spikes(gui)
        end

        gui.spikes
    end

    my_colors = map(my_spikes) do t

        waveform_color=[RGBA(0f0, 0f0, 1f0,1f0) for i=1:(size(t,2)-3)*size(t,1)]

        thres_color = [RGBA(0f0, 0f0, 0f0,1f0) for i=1:size(t,1)*3]

        [waveform_color; thres_color]
    end

    my_time = map(slider_value) do t

        minutes = div(t, 30000 * 60)
        seconds = round((t- minutes*60*30000)/30000,2)

        mystring=string("              ", minutes,"m ",seconds,"s")
    end

    gamma_slider, gamma_slider_s = labeled_slider(0.0f0:.1f0:1.0f0,gui.sliderscreen)

    change_gamma = map(gamma_slider_s) do t
       gui.gamma = t
    end

    controls = Pair[
        "play" => play_viz,
        "gamma" => gamma_slider,
    ]

    gui.p_slider=slider_value

    _view(visualize(my_animation, :lines,thickness=2f0), gui.datascreen)
    _view(visualize(my_spikes, :lines, thickness=1f0,color = my_colors), gui.imgscreen)
    _view(visualize(
        controls,
        text_scale = 4mm,
        width = 8*iconsize
    ), gui.sliderscreen, camera = :fixed_pixel)
    _view(visualize(my_time, color=RGBA(0f0,0f0,0f0)),gui.timescreen)

end

function add_labels(points,labels,myscreen)

    mypoints = GeometryTypes.Point2f0[points[i] for i=1:length(points)]

    (mytext,mywidths)=GLVisualize.annotated_text(mypoints,labels)

    _view(mytext,myscreen)
end


text:
abstract type AbstractSynthesizeMethod end

"""
    Raw <: AbstractSynthesizeMethod

Pixel values inside SuperPixel is not changed.
"""
struct Raw <: AbstractSynthesizeMethod end
"""
    Mean <: AbstractSynthesizeMethod

Pixel values inside SuperPixel is averaged.
"""
struct Average <: AbstractSynthesizeMethod end

"""
    synthesize([CT], img, [method=Raw()])

Collect pixels in `img` and returns an image of type `Array{CT}`.

Argument `method` controls how pixels are collected, supported methods
are [`Raw`](@ref) and [`Average`](@ref). The default value is `Raw`.

!!! note
    Overlaps between superpixels are averaged. Empty pixels are filled
    by zeros.
"""
synthesize(img::AbstractArray{<:SuperPixel}, method = Raw()) =
    synthesize(color_type(img), img, method)
synthesize(::Type{T}, img::AbstractArray{<:SuperPixel}, method = Raw()) where T =
    synthesize(T, img, method)

function synthesize(::Type{CT},
                    img::AbstractArray{<:SuperPixel},
                    ::Raw) where CT <: Union{AbstractRGB, AbstractGray}

    out = zeros(CT, imsize(img))
    count = zeros(Int, imsize(img)) # overlap count
    for SP in img
        isempty(SP) && continue
        R = position(SP)
        out[R] .+= color(SP)
        count[R] .+= 1
    end
    return out./count
end

function synthesize(::Type{CT},
                    img::AbstractArray{<:SuperPixel},
                    ::Average) where CT <: Union{AbstractRGB, AbstractGray}
    averaged_img = map(img) do SP
        SuperPixel(fill(mean(SP.color), size(SP.color)),
                   position(SP))
    end
    return synthesize(CT, averaged_img, Raw())
end