text:
isdense(::Type{<:DenseArray}) = true

# """
# ShortVector{T} simply wraps a Vector{T}, but uses a different hash function that is faster for short vectors to support using it as the keys of a Dict.
# This hash function scales O(N) with length of the vectors, so it is slow for long vectors.
# """
# struct ShortVector{T} <: DenseVector{T}
#     data::Vector{T}
# end
# Base.@propagate_inbounds Base.getindex(x::ShortVector, I...) = x.data[I...]
# Base.@propagate_inbounds Base.setindex!(x::ShortVector, v, I...) = x.data[I...] = v
# ShortVector{T}(::UndefInitializer, N::Integer) where {T} = ShortVector{T}(Vector{T}(undef, N))
# @inbounds Base.length(x::ShortVector) = length(x.data)
# @inbounds Base.size(x::ShortVector) = size(x.data)
# @inbounds Base.strides(x::ShortVector) = strides(x.data)
# @inbounds Base.push!(x::ShortVector, v) = push!(x.data, v)
# @inbounds Base.append!(x::ShortVector, v) = append!(x.data, v)
# function Base.hash(x::ShortVector, h::UInt)
#     @inbounds for n ∈ eachindex(x)
#         h = hash(x[n], h)
#     end
#     h
# end
# function Base.isequal(a::ShortVector{T}, b::ShortVector{T}) where {T}
#     length(a) == length(b) || return false
#     @inbounds for i ∈ 1:length(a)
#         a[i] === b[i] || return false
#     end
#     true
# end
# Base.convert(::Type{Vector}, sv::ShortVector) = sv.data
# Base.convert(::Type{Vector{T}}, sv::ShortVector{T}) where {T} = sv.data

struct Loop
    itersymbol::Symbol
    rangehint::Int
    rangesym::Symbol
    hintexact::Bool # if true, rangesym ignored and rangehint used for final lowering
end
function Loop(itersymbol::Symbol, rangehint::Int)
    Loop( itersymbol, rangehint, Symbol("##UNDEFINED##"), true )
end
function Loop(itersymbol::Symbol, rangesym::Symbol, rangehint::Int = 1_024)
    Loop( itersymbol, rangehint, rangesym, false )
end

# load/compute/store × isunroled × istiled × pre/post loop × Loop number
struct LoopOrder <: AbstractArray{Vector{Operation},5}
    oporder::Vector{Vector{Operation}}
    loopnames::Vector{Symbol}
    bestorder::Vector{Symbol}
end
function LoopOrder(N::Int)
    LoopOrder( [ Operation[] for i ∈ 1:24N ], Vector{Symbol}(undef, N), Vector{Symbol}(undef, N) )
end
LoopOrder() = LoopOrder(Vector{Operation}[],Symbol[],Symbol[])
Base.empty!(lo::LoopOrder) = foreach(empty!, lo.oporder)
function Base.resize!(lo::LoopOrder, N::Int)
    Nold = length(lo.loopnames)
    resize!(lo.oporder, 32N)
    for n ∈ 32Nold+1:32N
        lo.oporder[n] = Operation[]
    end
    resize!(lo.loopnames, N)
    resize!(lo.bestorder, N)
    lo
end
Base.size(lo::LoopOrder) = (4,2,2,2,length(lo.loopnames))
Base.@propagate_inbounds Base.getindex(lo::LoopOrder, i::Int) = lo.oporder[i]
Base.@propagate_inbounds Base.getindex(lo::LoopOrder, i...) = lo.oporder[LinearIndices(size(lo))[i...]]

# Must make it easy to iterate
# outer_reductions is a vector of indixes (within operation vectors) of the reduction operation, eg the vmuladd op in a dot product
struct LoopSet
    loops::Dict{Symbol,Loop} # sym === loops[sym].itersymbol
    opdict::Dict{Symbol,Operation}
    operations::Vector{Operation} # Split them to make it easier to iterate over just a subset
    outer_reductions::Vector{Int} # IDs of reduction operations that need to be reduced at end.
    loop_order::LoopOrder
    # stridesets::Dict{ShortVector{Symbol},ShortVector{Symbol}}
    preamble::Expr # TODO: add preamble to lowering
    includedarrays::Vector{Tuple{Symbol,Int}}
    syms_aliasing_refs::Vector{Symbol} # O(N) search is faster at small sizes
    refs_aliasing_syms::Vector{ArrayReference}
    cost_vec::Matrix{Float64}
    reg_pres::Matrix{Int}
    # sym_to_ref_aliases::Dict{Symbol,ArrayReference}
    # ref_to_sym_aliases::Dict{ArrayReference,Symbol}
end

function cost_vec_buf(ls::LoopSet)
    cv = @view(ls.cost_vec[:,2])
    @inbounds for i ∈ 1:4
        cv[i] = 0.0
    end
    cv
end
function reg_pres_buf(ls::LoopSet)
    ps = @view(ls.reg_pres[:,2])
    @inbounds for i ∈ 1:4
        ps[i] = 0
    end
    ps
end
function save_tilecost!(ls::LoopSet)
    @inbounds for i ∈ 1:4
        ls.cost_vec[i,1] = ls.cost_vec[i,2]
        ls.reg_pres[i,1] = ls.reg_pres[i,2]
    end
end

# function op_to_ref(ls::LoopSet, op::Operation)
    # s = op.variable
    # id = findfirst(ls.syms_aliasing_regs)
    # @assert id !== nothing
    # ls.refs_aliasing_syms[id]
# end
pushpreamble!(ls::LoopSet, ex) = push!(ls.preamble.args, ex)

function includesarray(ls::LoopSet, array::Symbol)
    for (a,i) ∈ ls.includedarrays
        a === array && return i
    end
    -1
end
function LoopSet()
    LoopSet(
        Dict{Symbol,Loop}(),
        Dict{Symbol,Operation}(),
        Operation[],
        Int[],
        LoopOrder(),
        Expr(:block,),
        Tuple{Symbol,Int}[],
        Symbol[],
        ArrayReference[],
        Matrix{Float64}(undef, 4, 2),
        Matrix{Int}(undef, 4, 2)
    )
end
num_loops(ls::LoopSet) = length(ls.loops)
function oporder(ls::LoopSet)
    N = length(ls.loop_order.loopnames)
    reshape(ls.loop_order.oporder, (4,2,2,2,N))
end
names(ls::LoopSet) = ls.loop_order.loopnames
isstaticloop(ls::LoopSet, s::Symbol) = ls.loops[s].hintexact
looprangehint(ls::LoopSet, s::Symbol) = ls.loops[s].rangehint
looprangesym(ls::LoopSet, s::Symbol) = ls.loops[s].rangesym
# itersyms(ls::LoopSet) = keys(ls.loops)
function getop(ls::LoopSet, var::Symbol, elementbytes::Int = 8)
    get!(ls.opdict, var) do
        # might add constant
        op = add_constant!(ls, var, elementbytes)
        pushpreamble!(ls, Expr(:(=), mangledvar(op), var))
        op
    end
end
function getop(ls::LoopSet, var::Symbol, deps, elementbytes::Int = 8)
    get!(ls.opdict, var) do
        # might add constant
        op = add_constant!(ls, var, deps, gensym(:constant), elementbytes)
        pushpreamble!(ls, Expr(:(=), mangledvar(op), var))
        op
    end
end
getop(ls::LoopSet, i::Int) = ls.operations[i + 1]

@inline extract_val(::Val{N}) where {N} = N
function vec_looprange(ls::LoopSet, s::Symbol, isunrolled::Bool, W::Symbol, U::Int, loop = ls.loops[s])
    incr = if isunrolled
        Expr(:call, lv(:valmuladd), W, U, -1)
    else
        Expr(:call, lv(:valsub), W, 1)
    end
    if loop.hintexact
        Expr(:call, :<, s, Expr(:call, :-, loop.rangehint, incr))
    else
        Expr(:call, :<, s, Expr(:call, :-, loop.rangesym, incr))
    end
end
function looprange(ls::LoopSet, s::Symbol, incr::Int = 1, mangledname::Symbol = s, loop = ls.loops[s])
    incr -= 1
    if iszero(incr)
        Expr(:call, :<, mangledname, loop.hintexact ? loop.rangehint : loop.rangesym)
    else
        Expr(:call, :<, mangledname, loop.hintexact ? loop.rangehint - incr : Expr(:call, :-, loop.rangesym, incr))
    end
end
# function looprange(ls::LoopSet, s::Symbol, incr::Expr, mangledname::Symbol = s, loop = ls.loops[s])
#     increxpr = Expr(:call, :-, incr, 1)
#     increxpr = if loop.hintexact
#         Expr(:call, :-, loop.rangehint, increxpr)
#     else
#         Expr(:call, :-, loop.rangesym, increxpr)
#     end
#     Expr(:call, :<, mangledname, increxpr)
# end

function Base.length(ls::LoopSet, is::Symbol)
    ls.loops[is].rangehint
end

function Operation(
    ls::LoopSet, variable, elementbytes, instruction,
    node_type, dependencies, reduced_deps, parents, ref = NOTAREFERENCE
)
    Operation(
        length(operations(ls)), variable, elementbytes, instruction,
        node_type, dependencies, reduced_deps, parents, ref
    )
end

# load_operations(ls::LoopSet) = ls.loadops
# compute_operations(ls::LoopSet) = ls.computeops
# store_operations(ls::LoopSet) = ls.storeops
# function operations(ls::LoopSet)
    # Base.Iterators.flatten((
        # load_operations(ls),
        # compute_operations(ls),
        # store_operations(ls)
    # ))
# end
operations(ls::LoopSet) = ls.operations
function pushop!(ls::LoopSet, op::Operation, var::Symbol = name(op))
    push!(ls.operations, op)
    ls.opdict[var] = op
    op
end
function add_block!(ls::LoopSet, ex::Expr, elementbytes::Int = 8)
    for x ∈ ex.args
        x isa Expr || continue # be that general?
        push!(ls, x, elementbytes)
    end
end
function register_single_loop!(ls::LoopSet, looprange::Expr)
    itersym = (looprange.args[1])::Symbol
    r = (looprange.args[2])::Expr
    @assert r.head === :call
    f = first(r.args)
    loop::Loop = if f === :(:)
        lower = r.args[2]
        upper = r.args[3]
        lii::Bool = lower isa Integer
        uii::Bool = upper isa Integer
        if lii & uii
            Loop(itersym, 1 + convert(Int,upper) - convert(Int,lower))
        else
            N = gensym(Symbol(:loop, itersym))
            ex = if lii
                if lower == 1
                    pushpreamble!(ls, Expr(:(=), N, upper))
                else
                    pushpreamble!(ls, Expr(:(=), N, Expr(:call, :-, upper, lower - 1)))
                end
            else
                ex = if uii
                    Expr(:call, :-, upper + 1, lower)
                else
                    Expr(:call, :-, Expr(:call, :+, upper, 1), lower)
                end
                pushpreamble!(ls, Expr(:(=), N, ex))
            end
            Loop(itersym, N)
        end
    elseif f === :eachindex
        N = gensym(Symbol(:loop, itersym))
        pushpreamble!(ls, Expr(:(=), N, Expr(:call, :length, r.args[2])))
        Loop(itersym, N)
    else
        throw("Unrecognized loop range type: $r.")
    end
    ls.loops[itersym] = loop
    nothing
end
function register_loop!(ls::LoopSet, looprange::Expr)
    if looprange.head === :block # multiple loops
        for lr ∈ looprange.args
            register_single_loop!(ls, lr)
        end
    else
        @assert looprange.head === :(=)
        register_single_loop!(ls, looprange)
    end
end
function add_loop!(ls::LoopSet, q::Expr, elementbytes::Int = 8)
    register_loop!(ls, q.args[1])
    body = q.args[2]
    if body.head === :block
        add_block!(ls, body, elementbytes)
    else
        push!(ls, q, elementbytes)
    end
end
function add_loop!(ls::LoopSet, loop::Loop)
    ls.loops[loop.itersym] = loop
end
function add_vptr!(ls::LoopSet, indexed::Symbol, id::Int, ptr::Symbol = Symbol("##vptr##_", indexed))
    if includesarray(ls, indexed) < 0
        push!(ls.includedarrays, (indexed, id))
        pushpreamble!(ls, Expr(:(=), ptr, Expr(:call, lv(:stridedpointer), indexed)))
    end
    nothing
end
function intersection(depsplus, ls)
    deps = Symbol[]
    for dep ∈ depsplus
        dep ∈ ls && push!(deps, dep)
    end
    deps
end
function loopdependencies(ref::ArrayReference, ls::LoopSet)
    deps = loopdependencies(ref)
    loopset = keys(ls.loops)
    for dep ∈ deps
        dep ∈ loopset || return intersection(deps, loopset)
    end
    deps
end
function add_load!(
    ls::LoopSet, var::Symbol, ref::ArrayReference, elementbytes::Int = 8
)
    if ref.loaded[] == true
        op = getop(ls, var, elementbytes)
        @assert var === op.variable
        return op
    end
    id = findfirst(r -> r == ref, ls.refs_aliasing_syms)
    if id === nothing
        push!(ls.syms_aliasing_refs, var)
        push!(ls.refs_aliasing_syms, ref)
    else
        opp = getop(ls, ls.syms_aliasing_refs[id], elementbytes)
        return isstore(opp) ? getop(ls, first(parents(opp))) : opp
    end
    ref.loaded[] = true
    # ls.sym_to_ref_aliases[ var ] = ref
    # ls.ref_to_sym_aliases[ ref ] = var
    op = Operation(
        length(operations(ls)), var, elementbytes,
        :getindex, memload, loopdependencies(ref, ls),
        NODEPENDENCY, NOPARENTS, ref
    )
    add_vptr!(ls, ref.array, identifier(op), ref.ptr)
    pushop!(ls, op, var)
end
function add_load_ref!(ls::LoopSet, var::Symbol, ex::Expr, elementbytes::Int = 8)
    ref = ref_from_ref(ex)
    add_load!(ls, var, ref, elementbytes)
end
function add_load_getindex!(ls::LoopSet, var::Symbol, ex::Expr, elementbytes::Int = 8)
    ref = ref_from_getindex(ex)
    add_load!(ls, var, ref, elementbytes)
end
function instruction(x)
    x isa Symbol ? x : last(x.args).value
end

function addsetv!(s::AbstractVector{T}, v::T) where {T}
    for sᵢ ∈ s
        sᵢ === v && return nothing
    end
    push!(s, v)
    nothing
end
function mergesetv!(s1::AbstractVector{T}, s2::AbstractVector{T}) where {T}
    for s ∈ s2
        addsetv!(s1, s)
    end
    nothing
end
function mergesetdiffv!(
    s1::AbstractVector{T},
    s2::AbstractVector{T},
    s3::AbstractVector{T}
) where {T}
    for s ∈ s2
        s ∉ s3 && addsetv!(s1, s)
    end
    nothing
end
function setdiffv!(s3::AbstractVector{T}, s1::AbstractVector{T}, s2::AbstractVector{T}) where {T}
    for s ∈ s1
        (s ∈ s2) || (s ∉ s3 && push!(s3, s))
    end
end
# This version has no dependencies, and thus will not be lowered
### if it is a literal, that literal is either var"##ZERO#Float##", var"##ONE#Float##", or has to have been assigned to var in the preamble.
# if it is a literal, that literal has to have been assigned to var in the preamble.
function add_constant!(ls::LoopSet, var::Symbol, elementbytes::Int = 8)
    pushop!(ls, Operation(length(operations(ls)), var, elementbytes, LOOPCONSTANT, constant, NODEPENDENCY, Symbol[], NOPARENTS), var)
end
function add_constant!(ls::LoopSet, var, elementbytes::Int = 8)
    sym = gensym(:temp)
    op = Operation(length(operations(ls)), sym, elementbytes, LOOPCONSTANT, constant, NODEPENDENCY, Symbol[], NOPARENTS)
    pushpreamble!(ls, Expr(:(=), mangledvar(op), var))
    pushop!(ls, op, sym)
end
# This version has loop dependencies. var gets assigned to sym when lowering.
function add_constant!(ls::LoopSet, var::Symbol, deps::Vector{Symbol}, sym::Symbol = gensym(:constant), elementbytes::Int = 8)
    # length(deps) == 0 && push!(ls.preamble.args, Expr(:(=), sym, var))
    pushop!(ls, Operation(length(operations(ls)), sym, elementbytes, var, constant, deps, NODEPENDENCY, NOPARENTS), sym)
end
function add_constant!(ls::LoopSet, var, deps::Vector{Symbol}, sym::Symbol = gensym(:constant), elementbytes::Int = 8)
    sym2 = gensym(:temp)
    pushpreamble!(ls, Expr(:(=), sym2, var))
    pushop!(ls, Operation(length(operations(ls)), sym, elementbytes, sym2, constant, deps, NODEPENDENCY, NOPARENTS), sym)
end
function pushparent!(parents::Vector{Operation}, deps::Vector{Symbol}, reduceddeps::Vector{Symbol}, parent::Operation)
    push!(parents, parent)
    mergesetdiffv!(deps, loopdependencies(parent), reduceddependencies(parent))
    if !(isload(parent) || isconstant(parent))
        mergesetv!(reduceddeps, reduceddependencies(parent))
    end
    nothing
end
function maybe_cse_load!(ls::LoopSet, expr::Expr, elementbytes::Int = 8)
    if expr.head === :ref
        offset = 0
        # array = first(expr.args)::Symbol
        # args = @view expr.args[2:end]
        # ref = ref_from_ref(expr)
    elseif expr.head === :call && first(expr.args) === :getindex
        offset = 1
        # array = (expr.args[2])::Symbol
        # args = @view expr.args[3:end]
        # ref = ref_from_getindex(expr)
    else
        return add_operation!(ls, gensym(:temporary), expr, elementbytes)
    end
    ref = ArrayReference(
        expr.args[1+offset],
        @view(expr.args[2+offset:end]),
        Ref(false)
    )::ArrayReference
    id = findfirst(r -> r == ref, ls.refs_aliasing_syms)
    if id === nothing
        add_load!( ls, gensym(:temporary), ref, elementbytes )
    else
        getop(ls, ls.syms_aliasing_refs[id], elementbytes)
    end
end
function add_parent!(
    parents::Vector{Operation}, deps::Vector{Symbol}, reduceddeps::Vector{Symbol}, ls::LoopSet, var, elementbytes::Int = 8
)
    parent = if var isa Symbol
        getop(ls, var, elementbytes)
    elseif var isa Expr #CSE candidate
        maybe_cse_load!(ls, var, elementbytes)
    else # assumed constant
        add_constant!(ls, var, elementbytes)
    end
    pushparent!(parents, deps, reduceddeps, parent)
end
function add_reduction!(
    parents::Vector{Operation}, deps::Vector{Symbol}, reduceddeps::Vector{Symbol}, ls::LoopSet, var::Symbol, elementbytes::Int = 8
)
    get!(ls.opdict, var) do
        add_constant!(ls, var, elementbytes)
    end
    # pushparent!(parents, deps, reduceddeps, parent)
end
function add_reduction_update_parent!(
    parents::Vector{Operation}, deps::Vector{Symbol}, reduceddeps::Vector{Symbol}, ls::LoopSet,
    var::Symbol, instr::Symbol, elementbytes::Int = 8
)
    parent = getop(ls, var, elementbytes)
    setdiffv!(reduceddeps, deps, loopdependencies(parent))
    pushparent!(parents, deps, reduceddeps, parent) # deps and reduced deps will not be disjoint
    op = Operation(length(operations(ls)), var, elementbytes, instr, compute, deps, reduceddeps, parents)
    parent.instruction === LOOPCONSTANT && push!(ls.outer_reductions, identifier(op))
    pushop!(ls, op, var) # note this overwrites the entry in the operations dict, but not the vector
end
function add_compute!(ls::LoopSet, var::Symbol, ex::Expr, elementbytes::Int = 8, ref = nothing)
    @assert ex.head === :call
    instr = instruction(first(ex.args))::Symbol
    args = @view(ex.args[2:end])
    parents = Operation[]
    deps = Symbol[]
    reduceddeps = Symbol[]
    # op = Operation( length(operations(ls)), var, elementbytes, instr, compute )
    reduction = false
    for arg ∈ args
        if var === arg
            reduction = true
            add_reduction!(parents, deps, reduceddeps, ls, arg, elementbytes)
        elseif ref == arg
            reduction = true
            add_load!(ls, var, ref, elementbytes)
        else
            add_parent!(parents, deps, reduceddeps, ls, arg, elementbytes)
        end
    end
    if reduction # arg[reduction] is the reduction
        add_reduction_update_parent!(parents, deps, reduceddeps, ls, var, instr, elementbytes)
    else
        op = Operation(length(operations(ls)), var, elementbytes, instr, compute, deps, reduceddeps, parents)
        pushop!(ls, op, var)
    end
end
function add_store!(
    ls::LoopSet, var::Symbol, ref::ArrayReference, elementbytes::Int = 8
)
    ldref = loopdependencies(ref)
    parent = getop(ls, var, ldref, elementbytes)
    pvar = parent.variable
    if pvar ∉ ls.syms_aliasing_refs
        push!(ls.syms_aliasing_refs, pvar)
        push!(ls.refs_aliasing_syms, ref)
    end
    op = Operation( length(operations(ls)), ref.array, elementbytes, :setindex!, memstore, ldref, reduceddependencies(parent), [parent], ref )
    add_vptr!(ls, ref.array, identifier(op), ref.ptr)
    pushop!(ls, op, ref.array)
end
function add_store_ref!(ls::LoopSet, var::Symbol, ex::Expr, elementbytes::Int = 8)
    ref = ref_from_ref(ex)::ArrayReference
    add_store!(ls, var, ref, elementbytes)
end
function add_store_setindex!(ls::LoopSet, ex::Expr, elementbytes::Int = 8)
    ref = ref_from_setindex(ex)::ArrayReference
    add_store!(ls, (ex.args[2])::Symbol, ref, elementbytes)
end
# add operation assigns X to var
function add_operation!(
    ls::LoopSet, LHS::Symbol, RHS::Expr, elementbytes::Int = 8
)
    if RHS.head === :ref
        add_load_ref!(ls, LHS, RHS, elementbytes)
    elseif RHS.head === :call
        f = first(RHS.args)
        if f === :getindex
            add_load_getindex!(ls, LHS, RHS, elementbytes)
        elseif f === :zero || f === :one
            c = gensym(:constant)
            pushpreamble!(ls, Expr(:(=), c, RHS))
            add_constant!(ls, c, [keys(ls.loops)...], LHS, elementbytes)
        else
            add_compute!(ls, LHS, RHS, elementbytes)
        end
    else
        throw("Expression not recognized:\n$x")
    end
end
function add_operation!(
    ls::LoopSet, LHS_sym::Symbol, RHS::Expr, LHS_ref::ArrayReference, elementbytes::Int = 8
)
    if RHS.head === :ref# || (RHS.head === :call && first(RHS.args) === :getindex)
        add_load!(ls, LHS_sym, LHS_ref, elementbytes)
    elseif RHS.head === :call
        f = first(RHS.args)
        if f === :getindex
            add_load!(ls, LHS_sym, LHS_ref, elementbytes)
        elseif f === :zero || f === :one
            c = gensym(:constant)
            pushpreamble!(ls, Expr(:(=), c, RHS))
            add_constant!(ls, c, [keys(ls.loops)...], LHS_sym, elementbytes)
        else
            add_compute!(ls, LHS_sym, RHS, elementbytes, LHS_ref)
        end
    else
        throw("Expression not recognized:\n$x")
    end
end
function Base.push!(ls::LoopSet, ex::Expr, elementbytes::Int = 8)
    if ex.head === :call
        finex = first(ex.args)::Symbol
        if finex === :setindex!
            add_store_setindex!(ls, ex, elementbytes)
        else
            throw("Function $finex not recognized.")
        end
    elseif ex.head === :(=)
        LHS = ex.args[1]
        RHS = ex.args[2]
        if LHS isa Symbol
            if RHS isa Expr
                add_operation!(ls, LHS, RHS, elementbytes)
            else
                add_constant!(ls, RHS, [keys(ls.loops)...], LHS, elementbytes)
            end
        elseif LHS isa Expr
            @assert LHS.head === :ref
            local lrhs::Symbol
            if RHS isa Symbol
                lrhs = RHS
            elseif RHS isa Expr
                # need to check of LHS appears in RHS
                # assign RHS to lrhs
                ref = ArrayReference(LHS)
                id = findfirst(r -> r == ref, ls.refs_aliasing_syms)
                lrhs = if id === nothing
                    gensym(:RHS)
                else
                    ls.syms_aliasing_refs[id]
                end
                add_operation!(ls, lrhs, RHS, ref, elementbytes)
            end
            add_store_ref!(ls, lrhs, LHS, elementbytes)
        else
            throw("LHS not understood:\n$LHS")
        end
    elseif ex.head === :block
        add_block!(ls, ex)
    elseif ex.head === :for
        add_loop!(ls, ex)
    else
        throw("Don't know how to handle expression:\n$ex")
    end
end

function fillorder!(ls::LoopSet, order::Vector{Symbol}, loopistiled::Bool)
    lo = ls.loop_order
    ro = lo.loopnames # reverse order; will have same order as lo
    nloops = length(order)
    if loopistiled
        tiled    = order[1]
        unrolled = order[2]
    else
        tiled = Symbol("##UNDEFINED##")
        unrolled = first(order)
    end
    ops = operations(ls)
    nops = length(ops)
    included_vars = fill(false, nops)
    place_after_loop = fill(true, nops)
    # to go inside out, we just have to include all those not-yet included depending on the current sym
    empty!(lo)
    for _n ∈ 1:nloops
        n = 1 + nloops - _n
        ro[_n] = loopsym = order[n]
        #loopsym = order[n]
        for (id,op) ∈ enumerate(ops)
            included_vars[id] && continue
            loopsym ∈ loopdependencies(op) || continue
            included_vars[id] = true
            isunrolled = (unrolled ∈ loopdependencies(op)) + 1
            istiled = (loopistiled ? (tiled ∈ loopdependencies(op)) : false) + 1
            optype = Int(op.node_type) + 1
            after_loop = place_after_loop[id] + 1
            push!(lo[optype,isunrolled,istiled,after_loop,_n], ops[id])
            set_upstream_family!(place_after_loop, op, false) # parents that have already been included are not moved, so no need to check included_vars to filter
        end
    end
end

# function depends_on_assigned(op::Operation, assigned::Vector{Bool})
#     for p ∈ parents(op)
#         p === op && continue # don't fall into recursive loop when we have updates, eg a = a + b
#         assigned[identifier(op)] && return true
#         depends_on_assigned(p, assigned) && return true
#     end
#     false
# end
# ind gets increased across tiles / unroll, so we need steps.
# function replace_ind_in_offset!(offset::Vector, op::Operation, ind::Int, t)
#     t == 0 && return nothing
#     var = op.variable
#     siter = op.symbolic_metadata[ind]
#     striden = op.numerical_metadata[ind]
#     strides = Symbol(:stride_, var)
#     offset[ind] = if tstriden == -1
#         Expr(:call, :*, Expr(:call, :+, strides, t), siter)
#     else
#         Expr(:call, :*, striden + t, siter)
#     end
#     nothing
# end


text:
module TestExtras
    using Base.Test
    using DataArrays
    using StatsBase
    using Compat

    ##########
    ## countmap
    ##########

    d = @data [NA,3,3]
    w = weights([1.1,2.2,3.3])
    cm = Dict{(@compat Union{Int, NAtype}), Int}([(NA, 1), (3, 2)])
    cmw = Dict{(@compat Union{Int, NAtype}), Real}([(NA, 1.1), (3, 5.5)])
    @assert isequal(countmap(d), cm)
    @assert isequal(countmap(d, w), cmw)

    ##########
    ## cut
    ##########

    @assert isequal(cut([2, 3, 5], [1, 3, 6]), PooledDataArray(["(1,3]", "(1,3]", "(3,6]"]))
    @assert isequal(cut([2, 3, 5], [3, 6]), PooledDataArray(["[2,3]", "[2,3]", "(3,6]"]))
    @assert isequal(cut([2, 3, 5, 6], [3, 6]), PooledDataArray(["[2,3]", "[2,3]", "(3,6]", "(3,6]"]))
    @assert isequal(cut([1, 2, 4], [1, 3, 6]), PooledDataArray(["[1,3]", "[1,3]", "(3,6]"]))
    @assert isequal(cut([1, 2, 4], [3, 6]), PooledDataArray(["[1,3]", "[1,3]", "(3,6]"]))
    @assert isequal(cut([1, 2, 4], [3]), PooledDataArray(["[1,3]", "[1,3]", "(3,4]"]))
    @assert isequal(cut([1, 5, 7], [3, 6]), PooledDataArray(["[1,3]", "(3,6]", "(6,7]"]))

    ages = [20, 22, 25, 27, 21, 23, 37, 31, 61, 45, 41, 32]
    bins = [18, 25, 35, 60, 100]
    cats = cut(ages, bins)
    pdv = PooledDataArray(["(18,25]", "(18,25]", "(18,25]",
                           "(25,35]", "(18,25]", "(18,25]",
                           "(35,60]", "(25,35]", "(60,100]",
                           "(35,60]", "(35,60]", "(25,35]"])
    @assert isequal(cats, pdv)

    ##########
    ## repeat
    ##########

    @test isequal(repeat(@data [3.0, 2.0, NA]; inner = 2, outer = 1),
                  @data [3.0, 3.0, 2.0, 2.0, NA, NA])
    @test isequal(repeat(@pdata ["a", "b", NA]; inner = 2, outer = 1),
                  @pdata ["a", "a", "b", "b", NA, NA])
    @test isequal(repeat(@data [1 2; 3 NA]; inner = [1, 2], outer = [2, 1]),
                  @data [1 1 2 2; 3 3 NA NA; 1 1 2 2; 3 3 NA NA])
    @test isequal(repeat(@pdata [:a :b NA]; inner = [2,1], outer = [1,3]),
                  @pdata [:a :b NA :a :b NA :a :b NA;
                          :a :b NA :a :b NA :a :b NA])
end

ext: jl


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

JLLWrappers.@generate_main_file_header("mold")
JLLWrappers.@generate_main_file("mold", UUID("0903f2b2-81ff-5be9-abf7-91be1d67eee8"))
end  # module mold_jll


text:
macro import_basic_huge()
    :(
        using ComplexVisual:
            CV_Error, cv_error,
            CV_TranslateByOffset, CV_MultiplyByFactor,
            CV_CyclicValue, cv_set_value!,
            cv_create_angle_cross_test, cv_half,
            CV_AttachType, cv_north, cv_south, cv_east, cv_west
     )
end

# {{{ Error Handling
struct CV_Error <: Exception
    msg :: String;
end
cv_error(messages::String...) = throw(CV_Error(join(messages)))

cv_error(something...) = cv_error(
    (@sprintf("%s", thing) for thing in something)...)
# }}}

# {{{ helper functions for showing concrete datatypes with their fields

"""
return fieldnames of a Datatype. Put the `parent_...` fields at the end.
"""
function cv_sorted_fieldnames(t::DataType)
    fnames = fieldnames(t)
    return tuple(
        filter(x -> !contains(string(x), "parent_"), fnames)...,
        filter(x ->  contains(string(x), "parent_"), fnames)...)
end

function cv_show_value_replacements(value)
    if value isa Cairo.CairoSurfaceImage  || value isa Cairo.CairoContext
        value = value.ptr
    end
    if value isa CV_LineSegments
        value = string(length(value)) * " CV_LineSegments"
    end
    return value
end

"""
implementation for `show(io, obj)` for ComplexVisual-objects with
`typeof(obj) isa DataType`
"""
function cv_show_impl(io::IO, obj)
    t = typeof(obj)::DataType
    print(io, string(t.name.name), '(')
    first = true
    fnames = cv_sorted_fieldnames(t)
    for name_sym in fnames
        value = getfield(obj, name_sym)
        vtype = typeof(value)
        value = cv_show_value_replacements(value)
        if !first
            print(io, ", ")
        end
        print(io, string(name_sym), ": ")
        if contains(string(vtype.name.name), '#')
            show(io, MIME("text/plain"), value)
        else
            show(io, value)
        end
        first = false
    end
    print(io, ')')
    return nothing
end

"""
implementation for `show(io, mime, obj)` for ComplexVisual-objects with
`typeof(obj) isa DataType` and where `mime` is "text/plain".
"""
function cv_show_impl(io::IO, m::MIME{Symbol("text/plain")}, obj)
    t = typeof(obj)
    outer_indent = (get(io, :cv_indent, "")::AbstractString)
    indent = outer_indent * "  "
    iio = IOContext(io, :cv_indent => indent)
    println(io, string(t.name.name), '(')
    fnames = cv_sorted_fieldnames(t)
    for name_sym in fnames
        print(io, indent, string(name_sym), ": ")
        show(iio, m, cv_show_value_replacements(getfield(obj, name_sym)))
        println(io)
    end
    print(io, outer_indent, ')')
    return nothing
end
# }}}

# {{{ Dynamic helpers
"""
A mutable struct with one (mutable) `value` (i.e. degree of freedom).
The struct is callable. When called, the input is translated by the
given `value`.
"""
mutable struct CV_TranslateByOffset{N<:Number} # {{{
    value    :: N

    function CV_TranslateByOffset(::Type{T}) where {T<:Number}
        return new{T}(zero(T))
    end
end
function (tbp::CV_TranslateByOffset{N})(z::N) where {N}
    return z + tbp.value
end # }}}

"""
A mutable struct with one `factor` as degree of freedom.
The struct is callable. When called, the input is multiplied by
the given `factor.

If the factor is a complex number then this can be used to code
a scaling and rotation operator.
"""
mutable struct CV_MultiplyByFactor{N<:Number} # {{{
    factor   :: N
    function CV_MultiplyByFactor(::Type{T}) where {T<:Number}
        return new{T}(one(T))
    end
end
function (mbf::CV_MultiplyByFactor{N})(z::N) where {N}
    return mbf.factor * z
end
# }}}

"""
A mutable struct with an integer `value` in the interval [1, maxvalue].
The struct is callable. When called, the value is increased by 1 or
set to 1 if maxvalue was reached.
"""
mutable struct CV_CyclicValue{maxV} # {{{
    value :: Int

    function CV_CyclicValue(maxV::Int)
        maxV < 1 && cv_error("CV_CyclicValue: maxV must be >= 1")
        return new{maxV}(1)
    end
end
function (cv::CV_CyclicValue{maxV})() where {maxV}
    cv.value = (cv.value == maxV) ? 1 : cv.value + 1
    return nothing
end

function cv_set_value!(cv::CV_CyclicValue{maxV}, new_value::Int) where {maxV}
    if !(1 ≤ new_value ≤ maxV)
        cv_error("cv_set_value: new_value must be in [1, ", string(maxV),
            "]; but found: ", string(new_value))
    end
    cv.value = new_value
    return nothing
end

# }}}

# }}}

# {{{ Helper functions
"""
creates a test, which checks if a line segment [z,w] intersects
the line segemnt [rmin*exp(iϕ), rmax*exp(iϕ)].
"""
function cv_create_angle_cross_test(ϕ::Real, rmin::Real, rmax::Real; δ=100eps())
    emϕ, rmin, rmax = exp(-Float64(ϕ)*1im), Float64(rmin), Float64(rmax)

    return (z, w) -> begin
        hz = z*emϕ
        rhz, ihz = real(hz), imag(hz)
        z == w &&  return (rmin ≤ rhz ≤ rmax) && (abs(ihz) ≤ δ)

        hw = w*emϕ
        rhw, ihw = real(hw), imag(hw)
        if (abs(ihz) ≤ δ) && (abs(ihw) ≤ δ)
            # both are very close to the cross-boundary
            return (rmin ≤ rhz ≤ rmax) || (rmin ≤ rhw ≤ rmax)
        else
            sign(ihz)*sign(ihw) != -1 && return false
            r = (rhz - rhw + rhw*ihz - rhz*ihw)/(ihz - ihw)
            return rmin ≤ r ≤ rmax
        end
    end
end

cv_half(x::N) where {N<:Integer}        =  x ÷ N(2)
cv_half(x::N) where {N<:AbstractFloat}  =  x/2

# }}}

# {{{ Directions/Locations
const cv_north, CV_northT = Val(:north), Val{:north}
const cv_south, CV_southT = Val(:south), Val{:south}
const cv_east, CV_eastT = Val(:east), Val{:east}
const cv_west, CV_westT = Val(:west), Val{:west}
const CV_AttachType = Union{CV_northT, CV_southT, CV_eastT, CV_westT}
# }}}

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
# Sim.jl

export simGPMTD!, simGPMTD_full!, forecast_sim!;

function simGPMTD!(n::Int, nburn::Int, intercept::InterceptNormal,
    mixcomps::Vector{MixComponentNormal},
    λ::Vector{T}) where T <: Real

    L = length(mixcomps)
    y = zeros(Float64, n + nburn + L)
    ζvec = zeros(Int, n + nburn + L)
    y[1:L] = sqrt(intercept.σ2 / 100.0) .* randn(L)
    xx = [ Matrix{Float64}(undef, 0, 1) for ℓ = 1:L ]
    nζ = zeros(Int, L+1)

    for t = (L + 1):(L + nburn + n)

        ζ = StatsBase.sample(0:L, Weights(λ))

        if ζ == 0

            y[t] = sqrt(intercept.σ2)*randn() + intercept.μ

        else

            if nζ[ζ+1] == 0

                f = sqrt(mixcomps[ζ].κ * mixcomps[ζ].σ2) * randn()
                mixcomps[ζ].fx[1] = deepcopy(f)
                y[t] = sqrt(mixcomps[ζ].σ2)*randn() + mixcomps[ζ].μ + f

            else

                mixcomps[ζ].D = expanD(mixcomps[ζ].D, y[t-ζ], xx[ζ])
                mixcomps[ζ].Cor = corrMat(mixcomps[ζ].D, mixcomps[ζ].corParams)

                Cov = PDMat_adj((mixcomps[ζ].κ * mixcomps[ζ].σ2) .* mixcomps[ζ].Cor)
                ncov = size(Cov.mat,1)
                f = rfullcond_fstar!(mixcomps[ζ].fx, Cov.mat[2:ncov, 2:ncov],
                    Cov.mat[1,1], Cov.mat[1,2:ncov], mixcomps[ζ].rng)
                pushfirst!(mixcomps[ζ].fx, f)

                y[t] = sqrt(mixcomps[ζ].σ2)*randn() + mixcomps[ζ].μ + f

            end

            # pushfirst!(xx[ζ][:,1], y[t-ζ])
            xx[ζ] = reshape(vcat(y[t-ζ], xx[ζ][:,1]), size(xx[ζ],1)+1, 1)
            nζ[ζ+1] += 1
            ζvec[t] = deepcopy(ζ)

        end
    end

    yout = y[(L + nburn + 1):(L + nburn + n)]
    ζout = ζvec[(L + nburn + 1):(L + nburn + n)]
    fx = [ deepcopy(mixcomps[ℓ].fx) for ℓ = 1:L ]

    return (yout, xx, fx, ζout)
end

function simGPMTD_full!(n::Int, nburn::Int, intercept::InterceptNormal,
    mixcomps::Vector{MixComponentNormal},
    λ::Vector{T}) where T <: Real

    L = length(mixcomps)
    y = zeros(Float64, n + nburn + L)
    X = zeros(Float64, n + nburn, L)
    ζvec = zeros(Int, n + nburn + L)
    y[1:L] = sqrt(intercept.σ2 / 100.0) .* randn(L)

    for ℓ = 1:L
        mixcomps[ℓ].fx[1] = sqrt(mixcomps[ℓ].κ * mixcomps[ℓ].σ2) * randn()
        X[1,ℓ] = deepcopy( y[L+1-ℓ] )
    end

    ζ = StatsBase.sample(0:L, Weights(λ))
    y[L+1] = sqrt(mixcomps[ζ].σ2)*randn() + mixcomps[ζ].μ + mixcomps[ζ].fx[1]
    ζvec[L+1] = deepcopy(ζ)

    for t = (L + 2):(L + nburn + n)
        ii = t - L

        for ℓ = 1:L
            X[ii,ℓ] = deepcopy( y[t-ℓ] )
            mixcomps[ℓ].D = expanD(mixcomps[ℓ].D, X[ii,ℓ], X[(ii-1):-1:1,[ℓ]])
            mixcomps[ℓ].Cor = corrMat(mixcomps[ℓ].D, mixcomps[ℓ].corParams)

            Cov = PDMat_adj((mixcomps[ℓ].κ * mixcomps[ℓ].σ2) .* mixcomps[ℓ].Cor)
            ncov = size(Cov.mat,1)
            f = rfullcond_fstar!(mixcomps[ℓ].fx, Cov.mat[2:ncov, 2:ncov],
                                 Cov.mat[1,1], Cov.mat[1,2:ncov], mixcomps[ℓ].rng)
            pushfirst!(mixcomps[ℓ].fx, f)
        end

        ζ = StatsBase.sample(0:L, Weights(λ))
        ζvec[t] = deepcopy(ζ)

        if ζ == 0
            y[t] = sqrt(intercept.σ2)*randn() + intercept.μ
        else
            y[t] = sqrt(mixcomps[ζ].σ2)*randn() + mixcomps[ζ].μ + mixcomps[ζ].fx[1]
        end

    end

    for ℓ = 1:L
        mixcomps[ℓ].fx = deepcopy(reverse(mixcomps[ℓ].fx)[(nburn+1):(nburn+n)])
    end

    yout = y[(L + nburn + 1):(L + nburn + n)]
    ζout = ζvec[(L + nburn + 1):(L + nburn + n)]
    Xout = X[(nburn+1):(nburn+n),:]

    return (yout, Xout, ζout)
end

## forecasting

function draw_f_next!(k::Int, DD::Matrix{T}, # XXk::T, XX1tok::Matrix{T},
    cc::T, mxcomp::Dict{Symbol, Any}, mod_mxcomp::MixComponentGPMTD, tol::T) where T <: Real
    ## for easy parallelization; each of these arguments is an element of a vector of length L

    # DD = expanD(DD, XXk, XX1tok, front=false)
    Cov = PDMat_adj(covMat(DD, cc, mxcomp[:corParams]), tol) # redundant, but not expensive
    ncov = size(Cov.mat, 1)

    f = rfullcond_fstar!(mxcomp[:fx], Cov.mat[1:(ncov-1), 1:(ncov-1)],
                         Cov.mat[ncov,ncov], Cov.mat[ncov,1:(ncov-1)], mod_mxcomp.rng, tol) # redundant and expensive, but necessary until we have a numerically stable shortcut

    push!(mxcomp[:fx], f)

    return mxcomp
end

function forecast_sim!(model::Model_GPMTD, K::Int, postsamp::Dict{Symbol, Any}; tol=1.0e-3, n_procs=1)
    ## modifies each rng of model.state.mixcomps

    mxcomp = deepcopy(postsamp[:mixcomps])
    yy = vcat( deepcopy(model.y), zeros(K) )
    XX = vcat( deepcopy(model.X), zeros(K, model.L) )

    λw = Weights(exp.(postsamp[:lλ]))
    ζ_new = [ StatsBase.sample(0:model.L, λw) for k = 1:K ]

    cc = [ mxcomp[j][:κ] * mxcomp[j][:σ2] for j = 1:model.L ]
    DD = [ deepcopy(model.D[j]) for j = 1:model.L ]

    ## appears to be totally unstable
      # C_inv = [ inv(PDMat_adj(covMat(model.D[j], cc[j], mxcomp[j][:corParams]), 1.0e-09).chol) for j = 1:model.L ] # inv(PDMat(C).chol) goes straight to inverse of C

    for k = 1:K
        XX[model.n + k, :] = deepcopy(yy[(model.n - 1 + k):-1:(model.n - model.L + k)])

        ## part of the unstable shortcut
        # dd = [ pairDistMat( XX[[model.n + k], [j]], XX[1:(model.n + k - 1), [j]] ) for j = 1:model.L ]
        # cc_vec = [ vec(covMat(dd[j], cc[j], mxcomp[j][:corParams])) for j = 1:model.L ]
        # sd_f = [ cc[j] - cc_vec[j]'C_inv[j]*cc_vec[j] for j = 1:model.L ]

        # XXk = [ deepcopy(XX[model.n+k, j]) for j = 1:model.L ]
        # XX1tok = [ deepcopy(XX[1:(model.n + k - 1), [j]]) for j = 1:model.L ]
        # DD = [ expanD(DD[j], XXk[j], XX1tok[j], front=false) for j = 1:model.L ]

        DD = [ expanD(DD[j], XX[model.n+k, j], XX[1:(model.n + k - 1), [j]], front=false) for j = 1:model.L ]

        if n_procs == 1
            for j = 1:model.L
                # DD[j] = expanD(DD[j], XXk[j], XX1tok[j], front=false)
                # Cov = PDMat_adj(covMat(DD[j], cc[j], mxcomp[j][:corParams]), tol)
                # ncov = size(Cov.mat, 1)
                #
                # f = rfullcond_fstar!(mxcomp[j][:fx], Cov.mat[1:(ncov-1), 1:(ncov-1)],
                #                      Cov.mat[ncov,ncov], Cov.mat[ncov,1:(ncov-1)], model.state.mixcomps[j].rng, tol)
                # push!(mxcomp[j][:fx], f)

                ## instead of above
                mxcomp[j] = draw_f_next!(k, DD[j], cc[j], mxcomp[j], model.state.mixcomps[j], tol) # apparently can't modify DD if I pass DD[j]?
            end
        elseif n_procs > 1 # parallelization does not solve the issue of large n which can really slow down forecasting, but hopefully speed things up modestly
            mxcomp = pmap(draw_f_next!, fill(k, model.L), DD, cc, mxcomp, model.state.mixcomps, fill(tol, model.L)) # pmap will update in place on each worker
        end

        if ζ_new[k] == 0
            yy[model.n + k] = sqrt(postsamp[:intercept].σ2) * randn() + postsamp[:intercept].μ
        else
            yy[model.n + k] = sqrt(mxcomp[ζ_new[k]][:σ2]) * randn() + mxcomp[ζ_new[k]][:μ] + mxcomp[ζ_new[k]][:fx][model.n + k]
        end
    end

    return yy[(model.n + 1):(model.n + K)]
end

function forecast_sim!(model::Model_GPMTD, K::Int, sims::Vector{Any}; tol=1.0e-3, n_procs=1)
    nsim = length(sims)
    Yforec = permutedims( hcat( [ forecast_sim!(model, K, sims[i], tol=tol, n_procs=n_procs) for i = 1:nsim ]... ) )
    return Yforec
end


text:
# ----------------------------------------
# 1. Frequency abstract types
# ----------------------------------------

"""
    Frequency

`Frequency` is an abstract type.

`Monthly`, `Quarterly`, `Yearly`, `Unit` abstract types are subtypes of `Frequency`.

### Examples
```julia-repl
julia> Quarterly <: Frequency
true
```
"""
abstract type Frequency end

"""
    Monthly

See also: [`Frequency`](@ref)
```
"""
abstract type Monthly <: Frequency end

"""
    Quarterly

See also: [`Frequency`](@ref)
```
"""
abstract type Quarterly <: Frequency end

"""
    Yearly

See also: [`Frequency`](@ref)
```
"""
abstract type Yearly <: Frequency end

"""
    Unit

See also: [`Frequency`](@ref)
```
"""
abstract type Unit <: Frequency end

"""
    ppy(::Frequency)

Returns the number of periods per year for a given `Frequency`, `MIT`, and `TSeries`

### Examples
```julia-repl
julia> ppy(Quarterly)                   # Frequency
4

julia> ppy(mm(2020, 1))                 # MIT
12

julia> ppy(TSeries(yy(2020), ones(3)))   # TSeries
1
```
"""
ppy(::Type{T}) where T <: Monthly   = 12
ppy(::Type{T}) where T <: Quarterly = 4
ppy(::Type{T}) where T <: Yearly    = 1
ppy(::Type{T}) where T <: Unit      = 1

# ----------------------------------------
# 2. MIT{T} - main workhorse of date declaration in TimeSeriesEcon.jl
# ----------------------------------------
"""
    MIT{Frequency}(x::Int64)

MIT is a primitive type (Signed 64) representing a discrete date.

__Note:__ Please use `yy`, `qq`, 'mm', `ii` to instantiate `MIT`s.
`MIT`s are mainly used internally for the package development.

### Examples
```julia-repl
julia> MIT{Yearly}(2000)
2000Y
julia> MIT{Monthly}(2020*12 + 3)
2020M4
```
"""
primitive type MIT{T <: Frequency} <: Signed 64 end

function MIT{T}(x::Int64) where T <: Frequency
    reinterpret(MIT{T}, x)
end

MIT{F}(y::Integer, p::Integer) where F <: Frequency = 1 <= p <= ppy(F) ?
                                                      reinterpret(MIT{F}, y * ppy(F) + p - 1) :
                                                      throw(ArgumentError("$F period must be in 1…$(ppy(F))"))

export frequencyof
"""
    frequencyof(::MIT)
    frequencyof(::Type{MIT})

Return the Frequency type of the given MIT instance of type.
"""
frequencyof(::Any) = nothing
frequencyof(::MIT{T}) where T <: Frequency = T
frequencyof(::Type{MIT{T}}) where T <: Frequency = T

# converting to Int is just the internal representation
Int64(x::MIT{T}) where T <: Frequency = reinterpret(Int64, x)
Base.promote_rule(::Type{MIT{S}}, ::Type{T}) where S <: Frequency where T <: Integer = Int64

# converting to Float has whole part equal to the year and fractional part
# equal to the fraction of the year for the period
(T::Type{<:AbstractFloat})(x::MIT{F}) where F <: Frequency = year(x) + (period(x) - 1) / ppy(F)
Base.promote_rule(::Type{MIT{F}}, ::Type{T}) where F <: Frequency where T <: AbstractFloat = T

Base.one(::MIT) = Int(1)
Base.one(::Type{<:MIT}) = Int(1)

"""
    year(x::MIT)

Return `::Int64` representing year for a given `MIT`

__Note:__ an internal method for now.

### Examples
```
julia> year(mm(2020, 1))
2020
```
"""
function year(x::MIT{T}) where T <: Frequency
    div(reinterpret(Int64, x), ppy(T))
end

"""
    period(x::MIT)

Return `::Int64` representing period for a given `MIT`

### Examples
```
julia> period(mm(2020, 1))
1
```
"""
function period(x::MIT{T}) where T <: Frequency
    1 + rem(reinterpret(Int64, x), ppy(T))
end

function Base.show(io::IO, x::MIT)
    F = first("$(frequencyof(x))")
    Y = year(x)
    Y = get(io, :compact, true) ? string(year(x)) : lpad(Y, 4)
    np = ppy(x)
    if np == 1
        print(io, Y, F)
    else
        n = length(string(np))
        print(io, Y, F, rpad(period(x), n))
    end
end

Base.string(m::MIT{T}) where T <: Frequency = repr(m)

# ----------------------------------------
# 2.1 MIT{T} constructors
#   - ii, mm, qq, yy
# ----------------------------------------

"""
    ii(x::Int64)

Represents an _Integer_ date and return `MIT{Unit}` type instance

### Examples
```julia-repl
julia> ii(123)
ii(123)

julia> ii(123) + 5
ii(128)

julia> typeof(ii(123))
MIT{Unit}
```
"""
ii(x::Int64) = MIT{Unit}(x, 1)

"""
    mm(y::Int64, p::Int64)

Represents a `Monthly` date and returns `MIT{Monthly}` type instance.

### Examples
```julia-repl
julia> mm(2020, 1)
2020M1

julia> mm(2020, 1) + 5
2020M6
```
"""
mm(y, p) = MIT{Monthly}(y, p)

"""
    qq(y::Int64, p::Int64)

Represents a `Quarterly` date and returns `MIT{Quarterly}` type instance.

### Examples
```julia-repl
julia> qq(2020, 1)
2020Q1

julia> qq(2020, 1) + 5
2021Q2
```
"""
qq(y, p) = MIT{Quarterly}(y, p)

"""
    yy(y::Int64)

Represents a `Yearly` date and returns `MIT{Yearly}` type instance.

### Examples
```julia-repl
julia> yy(2020)
2020Y

julia> yy(2020) + 5
2025Y
```
"""
yy(y, p=1) = MIT{Yearly}(y, p)

# ----------------------------------------
# 2.2 MIT{T} operations
# ----------------------------------------

Base.isequal(x::MIT, y::MIT) = frequencyof(x) == frequencyof(y) && isequal(Int(x), Int(y))
Base.:(==)(x::MIT, y::MIT) = frequencyof(x) == frequencyof(y) && Int(x) == Int(y)

Base.:(<)(x::MIT, y::MIT) = frequencyof(x) == frequencyof(y) ? Int(x) < Int(y) : throw(ArgumentError("Cannot compare MIT of different frequency."))
Base.:(<=)(x::MIT, y::MIT) = (x < y) || (x == y)

Base.:-(a::MIT, b::MIT) = frequencyof(a) == frequencyof(b) ? Int(a)-Int(b) : throw(ArgumentError("Cannot subtract MIT of different frequencies."))
Base.:-(m::MIT, x::Integer) = reinterpret(typeof(m), Int(m) - x)
Base.:-(::Integer, ::MIT) = throw(ArgumentError("Cannot subtract MIT from an Integer."))

Base.:+(m::MIT, x::Integer) = reinterpret(typeof(m), Int(m) + x)
Base.:+(x::Integer, m::MIT) = reinterpret(typeof(m), Int(m) + x)
Base.:+(::MIT, ::MIT) = throw(ArgumentError("Cannot add two MIT."))

# ----------------------------------------
# 2.2 MIT{T} vector and dict support
# ----------------------------------------

# added so MIT can be used as dictionary keys
Base.hash(x::MIT{T}) where T <: Frequency = begin
    reinterpret(Int64, x) |> hash
end

# added for sorting Vector{MIT{T}} where T <: Frequency
Base.sub_with_overflow(x::MIT{T}, y::MIT{T}) where T <: Frequency = begin
    Base.checked_ssub_int(reinterpret(Int64, x), reinterpret(Int64, y))
end

# ----------------------------------------
# 4 Convenience constants
# ----------------------------------------

struct _FPConst{F <: Frequency,P}
    _FPConst{F,P}() where F <: Frequency where P = 0 < P <= ppy(F) ? new{F,P}() : throw(ArgumentError("Invalid period $P for frequency $F."))
end
_FPConst(F::Type{<:Frequency}, p::Integer=1) = _FPConst{F,p}()
Base.:(*)(y::Integer, fc::_FPConst{F,P}) where F <: Frequency where P = MIT{F}(y, P)

export Y, U, Q1, Q2, Q3, Q4, M1, M2, M3, M4, M5, M6, M7, M8, M9, M10, M11, M12
global const Y = _FPConst(Yearly)
global const U = _FPConst(Unit)
global const Q1 = _FPConst(Quarterly, 1)
global const Q2 = _FPConst(Quarterly, 2)
global const Q3 = _FPConst(Quarterly, 3)
global const Q4 = _FPConst(Quarterly, 4)
global const M1 = _FPConst(Monthly, 1)
global const M2 = _FPConst(Monthly, 2)
global const M3 = _FPConst(Monthly, 3)
global const M4 = _FPConst(Monthly, 4)
global const M5 = _FPConst(Monthly, 5)
global const M6 = _FPConst(Monthly, 6)
global const M7 = _FPConst(Monthly, 7)
global const M8 = _FPConst(Monthly, 8)
global const M9 = _FPConst(Monthly, 9)
global const M10 = _FPConst(Monthly, 10)
global const M11 = _FPConst(Monthly, 11)
global const M12 = _FPConst(Monthly, 12)


text:
using CairoMakie, ColorSchemes

export plot_neutron_falcon_stats, plot_mass_stats, plot_other_stats, plot_rocket_specs

noto_sans = assetpath("fonts", "NotoSans-Regular.ttf")
noto_sans_bold = assetpath("fonts", "NotoSans-Bold.ttf")

function plot_neutron_falcon_stats()
    colors = Makie.wong_colors()

    fig = Figure(resolution = (1000, 700), font = noto_sans)
    grid = fig[1:3, 1:3] = GridLayout()

    heightaxis = Axis(grid[1:3, 1],
                      xticks = (1:3, ["Falcon 9", "Falcon 9 v1.0", "Neutron"]),
                      ylabel = "meters")

    # position of each bar along x-axis
    # repeated x values means several bars are placed in the same position
    xs =    [1, 1, 1, 2, 2, 2, 3, 3, 3]

    # indicate order for y-values placed on the same x-axis location.
    group = [1, 2, 3, 1, 2, 3, 1, 2, 3]

    # y-axis value for each bar
    heights = [70, 3.7, 5.2,  47.8, 3.7, 5.2, 40, 7, 5]
    barplot!(
        heightaxis,
        xs, heights,
        dodge = group, # cluster bars per property
        color = colors[group])

    labels = ["Height", "Diameter", "Fairing Diameter"]
    elements = [PolyElement(polycolor = colors[i]) for i in 1:length(labels)]
    Legend(grid[1, 2], elements, labels)
    fig
end

function plot_mass_stats()
    colors = Makie.wong_colors()

    fig = Figure(resolution = (1000, 700), font = noto_sans)
    grid = fig[1:3, 1:3] = GridLayout()

    massaxis = Axis(grid[1:3, 1],
                      xticks = (1:3, ["Falcon 9", "Falcon 9 v1.0", "Neutron"]),
                      ylabel = "tonne")

    # position of each bar along x-axis
    # the max value is 3, which suggests the x-axis goes from 1 - 3
    # the same x-value is repeated twice, which means two y-values get
    # placed on the same location at every x-value
    xs =    [1, 1, 2, 2, 3, 3]

    # assigning each bar to a group.
    group = [1, 2, 1, 2, 1, 2]

    # y-axis value for each bar
    masses = [549, 22.8, 333, 9, 480, 15]
    barplot!(
        massaxis,
        xs, masses,
        dodge = group, # cluster bars per property
        color = colors[group])

    labels = ["Mass", "Payload (LEO)"]
    elements = [PolyElement(polycolor = colors[i]) for i in 1:length(labels)]
    Legend(grid[1, 2], elements, labels)
    fig
end

function plot_other_stats()
    fig = Figure(resolution = (1000, 700), font = noto_sans)
    grid = fig[1:3, 1:3] = GridLayout()

    # massaxis = Axis(grid[1:3, 1],
    #                   xticks = (1:3, ["Falcon 9", "Falcon 9 v1.0", "Neutron"]),
    #                   ylabel = "Mass (tonne)")
    #
    # payloadaxis = Axis(grid[1:3, 1],
    #                 xticks = (1:3, ["Falcon 9", "Falcon 9 v1.0", "Neutron"]),
    #                 ylabel = "Payload (tonne)")
    #
    # masses = [549, 333, 480]
    #
    # barplot!(
    #     massaxis,
    #     xs, masses,
    #     dodge = group, # cluster bars per property
    #     color = colors[group])
end

function plot_rocket_specs()
    # Picked from catalog: https://juliagraphics.github.io/ColorSchemes.jl/stable/catalogue/
    colors = ColorSchemes.Spectral_4.colors

    df = load_rocket_specs()
    fig = Figure(resolution = (1024, 1024), font = noto_sans)

    cols = [:height, :diameter, :fairing, :mass, :payload, :thrust, :Isp]
    units = ["meter", "meter", "meter", "tonne", "tonne", "kilo Newton", "Secons"]

    # column position in grid for plot for a given index
    colpos = [1, 1, 1,  # meters
              2, 2,     # tonnes
              3, 3]
    rowpos = [1, 2, 3,
              1, 2,
              1, 2]

    rows = nrow(df)
    for (i, col) in enumerate(cols)
        ax = Axis(fig[colpos[i], rowpos[i]],
                  ylabel = units[i],    # label along y-axis
                  title = string(col),
                  xticks = 0:0          # num of tick marks and their spacing
             )
        # hidedecorations!(ax)
        # hidespines!(ax)
        barplot!(ax, 1:rows, df[:, col],
            color=colors[1:rows],
            # color_over_bar=:white,
            # label_size=12,
            # flip_labels_at=0,
            # font=noto_sans_bold,
            # bar_labels=:y
        )
    end
    for i in 1:3
        colsize!(fig.layout, i, Relative(1/3))
    end
    elements = polyelement.(colors[1:rows])
    Legend(fig[2, 3], elements, df[:, :rocket])

    fig
end


text:
module ImmutableArrays

export ImmutableArray, freeze, thaw

struct ImmutableArray{T,N} <: DenseArray{T,N}
    data::Array{T,N}
end
freeze(A::Array) = ImmutableArray(copy(A))
thaw(A::ImmutableArray) = copy(A.data)

Base.size(A::ImmutableArray) = size(A.data)
Base.@propagate_inbounds Base.getindex(A::ImmutableArray, I::Integer...) = getindex(A.data, I...)
Base.@propagate_inbounds Base.getindex(A::ImmutableArray, I...) = freeze(getindex(A.data, I...))
Base.strides(A::ImmutableArray) = strides(A.data)
Base.pointer(A::ImmutableArray) = pointer(A.data)
Base.pointer(A::ImmutableArray, i::Integer) = pointer(A.data, i)
Base.unsafe_convert(::Type{Ptr{T}}, A::ImmutableArray) where {T} = Base.unsafe_convert(Ptr{T}, A.data)
Base.reshape(A::ImmutableArray, args...) = freeze(reshape(A.data, args...)) # This could share data if possible
Base.reshape(A::ImmutableArray, arg::Tuple{Vararg{Int64,N}} where N) = freeze(reshape(A.data, arg))


text:
mutable struct HexGrid

    node_size::Float64
    grid::Array{Array{Float64}}

    function HexGrid(res::Float64; dilat::Float64 = 1.0, theta::Float64 = 0.0,
                     bor_x::Array{Float64,1} = [0.0,1.0], bor_y::Array{Float64,1} = [0.0,1.0],
                     node_size::Float64 = 1.0, scale::Bool = true)
        grid = ref_grid(res)
        if dilat != 1.0
            dilat_grid!(grid, dilat)
        end
        if theta != 0
            rotate_grid!(grid, theta)
        end
        if scale
            scale_grid!(grid,bor_x[1], bor_x[2], bor_y[1], bor_y[2])
        end
        new(node_size,grid)
    end

end

function ref_grid(hex_size::Float64)
    w = sqrt(3)*hex_size
    h = 2*hex_size
    hi,lo = 2,-1
    center_x = [i*w/2+lo for i in 0:div(hi-lo, w/2)]
    center_y = [j*0.75*h+lo for j in 0:div(hi-lo, 0.75*h)]
    grid = Array{Float64}[]
    for i in eachindex(center_x)
        for j in eachindex(center_y)
            if (i+j)%2==0; push!(grid, [center_x[i], center_y[j]]); end
        end
    end
    return grid
end

function rotate_grid!(hex_grid::Array{Array{Float64}}, theta::Float64)
    rot = [cos(theta) -sin(theta); sin(theta) cos(theta)]
    center = [0.5,0.5]
    for i in eachindex(hex_grid)
        hex_grid[i] .= rot*(hex_grid[i] .- center) .+ center
    end
    nothing
end

function dilat_grid!(hex_grid::Array{Array{Float64}},ratio::Float64)
    for i in eachindex(hex_grid)
        hex_grid[i] .*= ratio
    end
    nothing
end

function dilat_grid!(hex_grid::HexGrid,ratio::Float64)
    for i in eachindex(hex_grid.grid)
        hex_grid.grid[i] .*= ratio
    end
    hex_grid.node_size *= ratio
    nothing
end

function trunc_grid!(hex_grid::Array{Array{Float64}}, lo_x::Float64 = 0.0,
                     hi_x::Float64 = 1.0, lo_y::Float64 = 0.0,
                     hi_y::Float64 = 1.0)
    filter!(x->(lo_x<=x[1]<=hi_x) && (lo_y<=x[2]<=hi_y),hex_grid)
    nothing
end

function scale_grid!(hex_grid::Array{Array{Float64}}, lo_x::Float64, hi_x::Float64,
                     lo_y::Float64, hi_y::Float64)
    l1 = hi_x - lo_x
    l2 = hi_y - lo_y
    if l1 >= l2
        len = l1
        trunc_grid!(hex_grid,0.0,1.0,0.0,l2/l1)
    end
    if l2 > l1
        len = l2
        trunc_grid!(hex_grid,0.0,l1/l2)
    end
    lo = [lo_x,lo_y]
    for i in eachindex(hex_grid)
        hex_grid[i] .= hex_grid[i] .*len .+lo
    end
    nothing
end

function scale_grid!(hex_grid::HexGrid, lo_x::Float64, hi_x::Float64,
                     lo_y::Float64, hi_y::Float64)
    l1 = hi_x - lo_x
    l2 = hi_y - lo_y
    if l1 >= l2
        len = l1
        trunc_grid!(hex_grid.grid,0.0,1.0,0.0,l2/l1)
    end
    if l2 > l1
        len = l2
        trunc_grid!(hex_grid.grid,0.0,l1/l2)
    end
    lo = [lo_x,lo_y]
    for i in eachindex(hex_grid.grid)
        hex_grid.grid[i] .= hex_grid.grid[i] .*len .+lo
    end
    hex_grid.node_size *= len
    nothing
end

function dist(p1::Array{Float64},p2::Array{Float64})
    d = sqrt((p1[1]-p2[1])^2+(p1[2]-p2[2])^2)
    return d
end

# %% hexmodule struct
mutable struct HexModule

    n_grid::Int64
    grids::Array{HexGrid}

    function HexModule(res::Float64; dilat::Float64 = 1.0, theta::Float64 = 0.0,
                       bor_x::Array{Float64,1} = [0.0,1.0], bor_y::Array{Float64,1} = [0.0,1.0],
                       mod_size::Int64 = 2) # parameter
        ref = HexGrid(res, scale = false)
        grids = HexGrid[]
        for i in 1:mod_size^2
            new_grid = deepcopy(ref)
            trans_grid!(new_grid, i, mod_size)
            dilat_grid!(new_grid, dilat)
            rotate_grid!(new_grid.grid, theta)
            scale_grid!(new_grid,bor_x[1], bor_x[2], bor_y[1], bor_y[2])
            push!(grids, new_grid)
        end
        new(mod_size^2,grids)
    end
end

function trans_grid!(hex_grid::HexGrid, placement::Int64, mod_size::Int64)
    D = dist(hex_grid.grid[1],hex_grid.grid[2])
    d = D/sqrt(3)
    d_trans = d/(mod_size)
    theta_i = pi/3
    theta_j = -pi/3
    i = div(placement-1,mod_size)
    j = (placement-1)%mod_size
    translation_i = d_trans*[cos(theta_i), sin(theta_i)]
    translation_j = d_trans*[cos(theta_j), sin(theta_j)]
    for _ in 1:i
        for index in eachindex(hex_grid.grid)
            hex_grid.grid[index] .+= translation_i
        end
    end
    for _ in 1:j
        for index in eachindex(hex_grid.grid)
            hex_grid.grid[index] .+= translation_j
        end
    end
    hex_grid.node_size = d_trans/2
    nothing
end

# %% Visualization
using Plots
gr()

function scatter_grid(hex_grid::HexGrid)
    x = []
    y = []
    for point in hex_grid.grid
        push!(x,point[1])
        push!(y,point[2])
    end
    display(scatter(x,y,color = :red))
    nothing
end

function scatter_grid!(hex_grid::HexGrid)
    x = []
    y = []
    for point in hex_grid.grid
        push!(x,point[1])
        push!(y,point[2])
    end
    display(scatter!(x,y,legend = false, xlabel = "Variable 1", ylabel = "Variable 2"))
    nothing
end

function scatter_module(hex_mod::HexModule)
    for i in 1:hex_mod.n_grid
        scatter_grid!(hexmod.grids[i])
    end
end

# %% hexmodule test

borx = [0.0,1.0]
bory = [0.0,1.0]
hexmod = HexModule(0.03, dilat = 1.0, theta = pi/10, bor_x = borx,
                   bor_y = bory, mod_size = 2)

plot(size = (600,600))
scatter_grid!(hexmod.grids[1])
#
#
# # %%
#
# @show hexmod.grids[1].grid[1]
# @show hexmod.grids[1].grid[2]
# @show hexmod.grids[1].grid[3]
# @show hexmod.grids[1].grid[4]

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
function street_light(energy::Float64, cost::Float64)

    return 30*cost + min(50000,max(0,(energy - 30)))*cost

end

function VAT(energy::Float64, cost::Float64; rate::Float64 = 0.13, cutoff = 280)
    tax = 0.0

    if energy > cutoff
        tax = cost*rate
    end

    return tax
end

function firefighters(energy::Float64, cost::Float64)

    tax = 0.0

    if energy > 100
        tax = cost*0.0175
    elseif energy > 1750
        tax = 1750*(cost/energy)*0.0175
    end

    return tax
end

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
using BinaryBuilder

# Collection of sources required to build Ipopt
sources = [
    "https://www.coin-or.org/download/source/Ipopt/Ipopt-3.12.8.tgz" =>
    "62c6de314220851b8f4d6898b9ae8cf0a8f1e96b68429be1161f8550bb7ddb03",

    # Ipopt 3.12.8 uses these particular BLAS/LAPACK/MUMPS dependencies
    "http://www.coin-or.org/BuildTools/Blas/blas-20130815.tgz" =>
    "ea87df6dc44829ee0a1733226d130c550b17a0bc51c8dbfcd662fb15520b23b5",
    "http://www.coin-or.org/BuildTools/Lapack/lapack-3.4.2.tgz" =>
    "60a65daaf16ec315034675942618a2230521ea7adf85eea788ee54841072faf0",
    "http://mumps.enseeiht.fr/MUMPS_4.10.0.tar.gz" =>
    "d0f86f91a74c51a17a2ff1be9c9cee2338976f13a6d00896ba5b43a5ca05d933",

    # ASL sources
    "https://github.com/ampl/mp/archive/3.1.0.tar.gz" =>
    "587c1a88f4c8f57bef95b58a8586956145417c8039f59b1758365ccc5a309ae9",
    "https://github.com/staticfloat/mp-extra/archive/v3.1.0-2.tar.gz" =>
    "2f227175437f73d9237d3502aea2b4355b136e29054267ec0678a19b91e9236e",
]

# Bash recipe for building across all platforms
script = raw"""
set -e

# First, install ASL
cd $WORKSPACE/srcdir/mp-3.1.0

# Remove benchmarking library (this is already done on the latest
# ampl/mp master branch, but we don't use that, so backport the removal)
rm -rf thirdparty/benchmark
patch -p1 < $WORKSPACE/srcdir/mp-extra-3.1.0-2/no_benchmark.patch

# Build ASL
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix \
      -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain \
      -DRUN_HAVE_STD_REGEX=0 \
      -DRUN_HAVE_STEADY_CLOCK=0 \
      ../

# Copy over pregenerated files after building arithchk, so as to fake out cmake,
# because cmake will delete our arith.h
make arithchk VERBOSE=1
mkdir -p src/asl
cp -v $WORKSPACE/srcdir/mp-extra-3.1.0-2/expr-info.cc ../src/expr-info.cc
cp -v $WORKSPACE/srcdir/mp-extra-3.1.0-2/arith.h.${target} src/asl/arith.h

# Build and install ASL
make -j${nproc} VERBOSE=1
make install VERBOSE=1

# Install Third-Party sources as Ipopt desires...
cd $WORKSPACE/srcdir/Ipopt-*/ThirdParty
cp -r $WORKSPACE/srcdir/BLAS/*.f ./Blas/
cp -r $WORKSPACE/srcdir/lapack-* ./Lapack/LAPACK
cp -r $WORKSPACE/srcdir/MUMPS_* ./Mumps/MUMPS

# Next, install Ipopt
cd $WORKSPACE/srcdir/Ipopt-*

# The Ipopt buildsystem has a very old config.{sub,guess}.  Update those.
update_configure_scripts

# Build BLAS
(cd ThirdParty/Blas; \
    ./configure --prefix=$prefix --disable-shared --with-pic --host=$target; \
    make -j${nproc}; \
    make install)

# Build LAPACK
(cd ThirdParty/Lapack; \
    ./configure --prefix=$prefix --disable-shared --with-pic --host=$target; \
    make -j${nproc}; \
    make install)

(cd ThirdParty/Mumps; \
    patch -p0 < mumps.patch; \
    patch -p0 < mumps_mpi.patch; \
    mv MUMPS/libseq/mpi.h MUMPS/libseq/mumps_mpi.h; \
    ./configure --prefix=$prefix --disable-shared --with-pic --host=$target; \
    make -j${nproc}; \
    make install)

# Finally, build Ipopt itself.  Ipopt does some unusual things with pkg-config
# paths that don't play well with our definition of PKG_CONFIG_SYSROOT_DIR, so
# we have to define --with-mumps-incdir.  And we do so with extreme prejudice,
# and an extra helping of flag-injection hackiness.  We also need to give it
# asl since it doesn't know how to look for our updated asl automatically.
./configure --prefix=$prefix \
            lt_cv_deplibs_check_method=pass_all \
            --with-mumps-incdir="$(pwd)/ThirdParty/Mumps/MUMPS/include -I$(pwd)/ThirdParty/Mumps/MUMPS/libseq -DCOIN_USE_MUMPS_MPI_H" \
            --with-mumps-lib="-L$(pwd)/ThirdParty/Mumps/MUMPS/.libs -lcoinmumps" \
            --with-asl-lib="$prefix/lib/libasl.a" \
            --with-asl-incdir="$prefix/include/asl" \
            --host=$target

            #--with-blas="$prefix/lib/libcoinblas.a -lgfortran" \
            #--with-lapack="$prefix/lib/libcoinlapack.a" \

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 = [
    Windows(:i686),
    Windows(:x86_64),
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:aarch64, :glibc),
    Linux(:armv7l, :glibc),
    # ppc64le isn't working right now....
    #Linux(:powerpc64le, :glibc),
    MacOS()
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libipopt", :libipopt),
    ExecutableProduct(prefix, "ipopt", :amplexe),
]

# 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, "Ipopt", sources, script, platforms, products, dependencies)


text:
using Test, BenchmarkTools
using StaticKernels

using StaticKernels: Window

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 0.1

@testset "consistency" begin
    a = rand(10, 10)
    a2 = rand(10, 10)

    @testset "window" begin
        b = rand(size(a)...)
        c = CartesianIndex(2, 2)
        k = Kernel{(-1:1, -1:1)}(identity)
        w = Window{axes(k)}(k, a, c)

        # AbstractArray interface
        @test ndims(w) == ndims(k)
        @test eltype(w) == eltype(a)

        # getindex / setindex
        for i in CartesianIndices(w)
            @inferred w[i]
            @test w[i] == a[c + i]
            w[i] = b[c + i]
            @test a[c + i] == w[i]
        end

        aw = a[1:3, 1:3]
        # iteration
        @test_broken sum(v for v in w) == sum(aw)

        # mapreduce
        @test sum(w) == sum(aw)
        @test minimum(w) == minimum(aw)

        # wrong kernel axes
        let a = rand(1)
            k = Kernel{(0:0,)}(w -> w[1])
            @test_throws MethodError map(k, a)
        end
    end

    @testset "extended array" begin
        ae = extend(a, StaticKernels.ExtensionNothing())
        #@test isnothing(ae[0,0])
    end

    @testset "kernel" begin
        k = Kernel{(-1:1,-1:1)}(@inline function(w) sum(Tuple(w)) end)
        k2 = Kernel{(-1:1,-1:1)}(@inline function(w1, w2) sum(Tuple(w1)) + prod(Tuple(w2)) end)
        @test eltype(k, a) == eltype(a)
        @test eltype(k2, a, a2) == eltype(a)
    end

    @testset "map" begin
        k = Kernel{(1:1, -1:-1)}(w -> w[1, -1])
        @test map(k, a) == a

        k = Kernel{(-1:1,0:1)}(w -> w[-1,0] + w[1,1])
        @test map(k, a) == a[1:end-2,1:end-1] .+ a[3:end,2:end]

        k = Kernel{(-1:1,0:1)}((w1, w2) -> w1[-1,0] + w2[1,1])
        @test map(k, a, a2) == a[1:end-2,1:end-1] .+ a2[3:end,2:end]
    end

    @testset "mapreduce" begin
        k = Kernel{(-1:1,0:1)}(w -> w[-1,0] + w[1,1])
        @test mapreduce(k, +, a) ≈ sum(a[1:end-2,1:end-1] .+ a[3:end,2:end])

        k = Kernel{(-1:1,0:1)}((w1, w2) -> w1[-1,0] + w2[1,1])
        @test mapreduce(k, +, a, a2) ≈ sum(a[1:end-2,1:end-1] .+ a2[3:end,2:end])

        k = Kernel{(0:1, 0:0)}(w -> sum(Tuple(w)))
        @test sum(k, a) ≈ sum(a[1:end-1,:]) + sum(a[2:end,:])
    end

    @testset "kernel macro" begin
        b = vec(a)

        k1 = @kernel w -> w[0]
        k2 = Kernel{(0:0,)}(w -> w[0])
        @test map(k1, b) == map(k2, b)

        k1 = @kernel w -> w[1] + w[-2]
        k2 = Kernel{(-2:1,)}(w -> w[1] + w[-2])
        @test map(k1, b) == map(k2, b)

        k1 = @kernel w -> w[1,2] + w[-3,5]
        k2 = Kernel{(-3:1,2:5)}(w -> w[1,2] + w[-3,5])
        @test map(k1, a) == map(k2, a)

        @test_throws Exception @macroexpand @kernel w -> w
        @test_throws Exception @macroexpand @kernel w -> w[1] + w[2,3]
    end

    @testset "gradient" begin
        grad = Kernel{(0:1,0:1)}(w -> (w[1,0] - w[0,0], w[0,1] - w[0,0]))

        grada = map(grad, a)
        gx = axes(grad, a)
        @test diff(a, dims=1)[gx...] ≈ [x[1] for x in grada[gx...]]
        @test diff(a, dims=2)[gx...] ≈ [x[2] for x in grada[gx...]]
    end

    @testset "extensions" begin
        ae = extend(a, StaticKernels.ExtensionNone())
        k = Kernel{(-1:1, -1:1)}(w -> sum(Tuple(w)))
        @test size(map(k, ae)) == size(k, ae)

        ae = extend(a, StaticKernels.ExtensionNothing())
        k = Kernel{(-1:1, -1:1)}(w -> w[1,0] + w[0,0])
        @test_throws MethodError map(k, ae)
        k = Kernel{(-1:1, -1:1)}(w -> something(w[1,0], 0))
        @test size(map(k, ae)) == size(k, ae)
        @test all(map(k, ae)[end, :] .== 0)

        ae = extend(a, StaticKernels.ExtensionReplicate())
        k = Kernel{(-1:1, -1:1)}(w -> w[1,0])
        @test size(map(k, ae)) == size(k, ae)
        @test all(map(k, ae)[end, :] .== ae[end, :])

        ae = extend(a, StaticKernels.ExtensionCircular())
        k = Kernel{(-1:1, -1:1)}(w -> w[1,0])
        @test size(map(k, ae)) == size(k, ae)
        @test all(map(k, ae)[end, :] .== ae[begin, :])

        ae = extend(a, StaticKernels.ExtensionSymmetric())
        k = Kernel{(-1:1, -1:1)}(w -> w[1,0])
        @test size(map(k, ae)) == size(k, ae)
        @test all(map(k, ae)[end, :] .== ae[end-1, :])

        ae = extend(a, StaticKernels.ExtensionConstant(0.))
        k = Kernel{(-1:1, -1:1)}(w -> w[1,0])
        @test eltype(k, ae) == eltype(ae)
        @test size(map(k, ae)) == size(k, ae)
        @test all(map(k, ae)[end, :] .== 0)

        ae = extend(a, StaticKernels.ExtensionConstant(missing))
        k = Kernel{(-1:1, -1:1)}(w -> w[1,0])
        @test eltype(k, ae) == Union{Missing,eltype(ae)}
    end
end

@testset "type stability" begin
    a = rand(100)

    @testset "window" begin
        k = Kernel{(0:0,)}(w -> w)

        @test eltype(k, a) <: Window
        @test isconcretetype(eltype(k, a))
    end

    @testset "kernel evaluation" begin
        k = Kernel{(0:0,)}(w -> w[0])
        @test eltype(k, a) == eltype(a)

        k = Kernel{(0:0,)}((w1, w2) -> w1[0] + w2[0])
        @test eltype(k, a, a) == eltype(a)
        @test eltype(k, a, round.(Int, a)) == Base.promote_type(eltype(a), Int)

        k = Kernel{(0:0,)}(w -> Tuple(w))
        @test eltype(k, a) <: Tuple
        @test isconcretetype(eltype(k, a))
    end
end

@testset "memory allocations" begin
    a = rand(100)
    a2 = rand(100)
    ks = [
        (Kernel{(0:0,)}(w -> w[0]), StaticKernels.ExtensionNone()),
        (Kernel{(0:1,)}(w -> something(w[1], 0.)), StaticKernels.ExtensionNothing()),
        (Kernel{(0:1,)}(w -> w[1]), StaticKernels.ExtensionReplicate()),
        (Kernel{(0:1,)}(w -> w[1]), StaticKernels.ExtensionCircular()),
        (Kernel{(0:1,)}(w -> w[1]), StaticKernels.ExtensionSymmetric()),
        (Kernel{(0:1,)}(w -> w[1]), StaticKernels.ExtensionConstant(0.)),
        (Kernel{(0:1,)}(w -> Tuple(w)), StaticKernels.ExtensionNone())]

    @testset "map! $(x[2])" for x in ks
        k, extension = x
        ae = extend(a, extension)
        b = similar(ae, eltype(k, ae), size(k, ae))
        @test 0 == @ballocated map!($k, $b, $ae)
    end

    @testset "sum $(x[2])" for x in ks[1:6]
        k, extension = x
        ae = extend(a, extension)
        @test eltype(k, ae) == eltype(a)
        @test 0 == @ballocated sum($k, $ae)
    end

    @testset "multiargument" begin
        k = Kernel{(-1:1,)}((w1, w2) -> w1[-1] + w2[1])
        b = similar(a, eltype(k, a, a2), size(k, a))
        @test 0 == @ballocated map!($k, $b, $a, $a2)
        @test 0 == @ballocated mapreduce($k, +, $a, $a2)
    end
end

@testset "performance" begin
    a = rand(1000000)
    a2 = rand(1000, 1000)
    a2ext = extend(rand(1000, 1000), StaticKernels.ExtensionConstant(0))

    @testset "Base.diff" begin
        k = Kernel{(0:1,)}(@inline function(w) @inbounds w[1] - w[0] end)

        @test 1.3 > @belapsed(map($k, $a)) / @belapsed(diff($a))
    end

    @testset "Base.map" begin
        k = Kernel{(0:0,)}(@inline function(w) @inbounds w[0] end)
        b = similar(a, size(k, a))

        @test 1.1 > @belapsed(map!($k, $b, $a)) / @belapsed(map!(identity, $b, $a))
    end

    @testset "window mapreduce" begin
        k1 = Kernel{(-1:1,)}(@inline function(w) @inbounds sum(w) end)
        k2 = Kernel{(-1:1,)}(@inline function(w) @inbounds w[-1] + w[0] + w[1] end)
        b = similar(a, size(k1, a))

        @test 1.1 > @belapsed(map!($k2, $b, $a)) / @belapsed(map!($k1, $b, $a))
    end

    @testset "extension" begin
        k = Kernel{(-1:1,-1:1)}(
            @inline function(w) @inbounds w[0,-1] + w[-1,0] - 4*w[0,0] + w[1,0] + w[0,1] end)
        b2 = similar(a2, size(k, a2))
        b2ext = similar(a2ext, size(k, a2ext))

        @test 1.1 > @belapsed(map!($k, $b2ext, $a2ext)) / @belapsed(map!($k, $b2, $a2))
    end
end


text:
# Copyright 2018 IBM Corp.
#
# Licensed 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.

module ConceptRDF
export presentation_to_rdf, generator_rdf_node

using Serd
using Catlab

using ...Doctrine, ...Ontology
using ..RDFUtils

const R = RDF.Resource

""" Convert concepts into RDF/OWL ontology.

Concepts are repesented as individuals of the OWL class `Concept`.
"""
function presentation_to_rdf(pres::Presentation, prefix::RDF.Prefix;
                             extra_rdf::Union{Function,Nothing}=nothing)
  stmts = RDF.Statement[
    RDF.Prefix("rdf"), RDF.Prefix("rdfs"), RDF.Prefix("owl"),
    RDF.Prefix("list", "http://www.co-ode.org/ontologies/list.owl#"),
    RDF.Prefix("monocl", "https://www.datascienceontology.org/ns/monocl/"),
    prefix
  ]
  for expr in generators(pres)
    append!(stmts, expr_to_rdf(expr, prefix))
    if extra_rdf != nothing && first(expr) != nothing
      append!(stmts, extra_rdf(expr, generator_rdf_node(expr, prefix)))
    end
  end
  stmts
end

""" Generate RDF for object generator.
"""
function expr_to_rdf(ob::Monocl.Ob{:generator}, prefix::RDF.Prefix)
  node = generator_rdf_node(ob, prefix)
  [ RDF.Triple(node, R("rdf","type"), R("monocl","TypeConcept")),
    RDF.Triple(node, R("monocl","id"), RDF.Literal(string(first(ob)))) ]
end

""" Generate RDF for subobject relation.
"""
function expr_to_rdf(sub::Monocl.SubOb, prefix::RDF.Prefix)
  dom_node = generator_rdf_node(dom(sub), prefix)
  codom_node = generator_rdf_node(codom(sub), prefix)
  [ RDF.Triple(dom_node, R("monocl","subtypeOf"), codom_node) ]
end

""" Generate RDF for morphism generator.

The domain and codomain objects are represented as OWL lists.
"""
function expr_to_rdf(hom::Monocl.Hom{:generator}, prefix::RDF.Prefix)
  node = generator_rdf_node(hom, prefix)
  cell_node = name -> R(prefix.name, "$(node.name):$name")
  input_nodes = [ generator_rdf_node(ob, prefix) for ob in collect(dom(hom)) ]
  output_nodes = [ generator_rdf_node(ob, prefix) for ob in collect(codom(hom)) ]
  nin, nout = length(input_nodes), length(output_nodes)
  RDF.Statement[
    [ RDF.Triple(node, R("rdf","type"), R("monocl","FunctionConcept")),
      RDF.Triple(node, R("monocl","id"), RDF.Literal(string(first(hom)))) ];
    owl_inputs_outputs(node, cell_node, nin, nout, index=true) do cell, is_input, i
      gen_node = (is_input ? input_nodes : output_nodes)[i]
      [ RDF.Triple(cell, R("list","hasContents"), gen_node),
        RDF.Triple(cell, R("monocl","isConcept"), gen_node) ]
    end;
  ]
end

""" Generate RDF for submorphism relation.
"""
function expr_to_rdf(sub::Monocl.SubHom, prefix::RDF.Prefix)
  dom_node = generator_rdf_node(dom(sub), prefix)
  codom_node = generator_rdf_node(codom(sub), prefix)
  [ RDF.Triple(dom_node, R("monocl","subfunctionOf"), codom_node) ]
end

""" Create RDF node for generator expression.
"""
function generator_rdf_node(expr::GATExpr{:generator}, prefix::RDF.Prefix)
  @assert first(expr) != nothing
  R(prefix.name, string(first(expr)))
end

end

max_stars_repo_name: cclauss/semanticflowgraph

max_issues_repo_name: cclauss/semanticflowgraph

max_forks_repo_name: cclauss/semanticflowgraph


text:
# You cannot have multiple double quotes like so:
# x = "He said, "Hello, world""
# println(x)
# If we use println(x), it will give us an error

# Neither $ or \ come up literally if they're in a string
# We need to use \ so the double quotes don't give a syntax error
y = "He said, \"Hello, world\""
println(y)
# If we do not use the println() function, the formatting does not work and you get the string "He said, \"Hello, world\""

# We can use \n as "Enter" or a new line in a string
z = "This is a \nnew line"
println(z)

poem = "Errors are red, \nSome things are blue, \nI love coding, \nAnd so should you"
println(poem)

# \t inserts empty space, like a tab (it is 8 characters wide)
tabexample = "I will insert a tab right \t <- here."
println(tabexample)

# \\ can be used if you actually need the backslash in the string
backslash = "This is a backslash: \\"
println(backslash)

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
using MyGraph

G = @buildgraph false 1,6,14 1,3,9 1,2,7 2,3,10 2,4,15 3,6,2 3,4,11 4,5,6 5,6,9

D1 = dijkstra(G;source = 1)
D2 = dijkstra(G;source = 2)
D3 = dijkstra(G;source = 3)
D4 = dijkstra(G;source = 4)
D5 = dijkstra(G;source = 5)
D6 = dijkstra(G;source = 6)