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

## Type aliases for convenience ##

typealias AbstractVector{T} AbstractArray{T,1}
typealias AbstractMatrix{T} AbstractArray{T,2}
typealias AbstractVecOrMat{T} Union{AbstractVector{T}, AbstractMatrix{T}}
typealias RangeIndex Union{Int, Range{Int}, UnitRange{Int}, Colon}

## Basic functions ##

vect() = Array(Any, 0)
vect{T}(X::T...) = T[ X[i] for i=1:length(X) ]

const _oldstyle_array_vcat_ = true

if _oldstyle_array_vcat_
    function oldstyle_vcat_warning(n::Int)
        if n == 1
            before = "[a]"
            after  = "collect(a)"
        elseif n == 2
            before = "[a,b]"
            after  = "[a;b]"
        else
            before = "[a,b,...]"
            after  = "[a;b;...]"
        end
        depwarn("$before concatenation is deprecated; use $after instead", :vect)
    end
    function vect(A::AbstractArray...)
        oldstyle_vcat_warning(length(A))
        vcat(A...)
    end
    function vect(X...)
        for a in X
            if typeof(a) <: AbstractArray
                oldstyle_vcat_warning(length(X))
                break
            end
        end
        vcat(X...)
    end
else
    function vect(X...)
        T = promote_typeof(X...)
        #T[ X[i] for i=1:length(X) ]
        # TODO: this is currently much faster. should figure out why. not clear.
        copy!(Array(T,length(X)), X)
    end
end

size{T,n}(t::AbstractArray{T,n}, d) = d <= n ? size(t)[d] : 1
size(x, d1::Integer, d2::Integer, dx::Integer...) = tuple(size(x, d1), size(x, d2, dx...)...)
eltype{T}(::Type{AbstractArray{T}}) = T
eltype{T,n}(::Type{AbstractArray{T,n}}) = T
elsize{T}(::AbstractArray{T}) = sizeof(T)
ndims{T,n}(::AbstractArray{T,n}) = n
ndims{T,n}(::Type{AbstractArray{T,n}}) = n
ndims{T<:AbstractArray}(::Type{T}) = ndims(super(T))
length(t::AbstractArray) = prod(size(t))::Int
endof(a::AbstractArray) = length(a)
first(a::AbstractArray) = a[first(eachindex(a))]

function first(itr)
    state = start(itr)
    done(itr, state) && throw(ArgumentError("collection must be non-empty"))
    next(itr, state)[1]
end
last(a) = a[end]

function stride(a::AbstractArray, i::Integer)
    if i > ndims(a)
        return length(a)
    end
    s = 1
    for n=1:(i-1)
        s *= size(a, n)
    end
    return s
end

strides(a::AbstractArray) = ntuple(i->stride(a,i), ndims(a))::Dims

function isassigned(a::AbstractArray, i::Int...)
    # TODO
    try
        a[i...]
        true
    catch
        false
    end
end

# used to compute "end" for last index
function trailingsize(A, n)
    s = 1
    for i=n:ndims(A)
        s *= size(A,i)
    end
    return s
end

## Traits for array types ##

abstract LinearIndexing
immutable LinearFast <: LinearIndexing end
immutable LinearSlow <: LinearIndexing end

linearindexing(A::AbstractArray) = linearindexing(typeof(A))
linearindexing{T<:AbstractArray}(::Type{T}) = LinearSlow()
linearindexing{T<:Array}(::Type{T}) = LinearFast()
linearindexing{T<:Range}(::Type{T}) = LinearFast()

linearindexing(A::AbstractArray, B::AbstractArray) = linearindexing(linearindexing(A), linearindexing(B))
linearindexing(A::AbstractArray, B::AbstractArray...) = linearindexing(linearindexing(A), linearindexing(B...))
linearindexing(::LinearFast, ::LinearFast) = LinearFast()
linearindexing(::LinearIndexing, ::LinearIndexing) = LinearSlow()

# The real @inline macro is not available this early in the bootstrap, so this
# internal macro splices the meta Expr directly into the function body.
macro _inline_meta()
    Expr(:meta, :inline)
end
macro _noinline_meta()
    Expr(:meta, :noinline)
end

## Bounds checking ##
@generated function trailingsize{T,N,n}(A::AbstractArray{T,N}, ::Type{Val{n}})
    n > N && return 1
    ex = :(size(A, $n))
    for m = n+1:N
        ex = :($ex * size(A, $m))
    end
    Expr(:block, Expr(:meta, :inline), ex)
end

checkbounds(::Type{Bool}, sz::Integer, i) = throw(ArgumentError("unable to check bounds for indices of type $(typeof(i))"))
checkbounds(::Type{Bool}, sz::Integer, i::Real) = 1 <= i <= sz
checkbounds(::Type{Bool}, sz::Integer, ::Colon) = true
function checkbounds(::Type{Bool}, sz::Integer, r::Range)
    @_inline_meta
    isempty(r) || (checkbounds(Bool, sz, minimum(r)) && checkbounds(Bool, sz, maximum(r)))
end
checkbounds(::Type{Bool}, sz::Integer, I::AbstractArray{Bool}) = length(I) == sz
function checkbounds(::Type{Bool}, sz::Integer, I::AbstractArray)
    @_inline_meta
    b = true
    for i in I
        b &= checkbounds(Bool, sz, i)
    end
    b
end
# Prevent allocation of a GC frame by hiding the BoundsError in a noinline function
throw_boundserror(A, I) = (@_noinline_meta; throw(BoundsError(A, I)))

# Don't define index types on checkbounds to make extending easier
checkbounds(A::AbstractArray, I...) = (@_inline_meta; _internal_checkbounds(A, I...))
# The internal function is named _internal_checkbounds since there had been a
# _checkbounds previously that meant something different.
_internal_checkbounds(A::AbstractArray, I::AbstractArray{Bool}) = size(A) == size(I) || throw_boundserror(A, I)
_internal_checkbounds(A::AbstractArray, I::AbstractVector{Bool}) = length(A) == length(I) || throw_boundserror(A, I)
_internal_checkbounds(A::AbstractArray, I) = (@_inline_meta; checkbounds(Bool, length(A), I) || throw_boundserror(A, I))
function _internal_checkbounds(A::AbstractMatrix, I, J)
    @_inline_meta
    (checkbounds(Bool, size(A,1), I) && checkbounds(Bool, size(A,2), J)) ||
        throw_boundserror(A, (I, J))
end
function _internal_checkbounds(A::AbstractArray, I, J)
    @_inline_meta
    (checkbounds(Bool, size(A,1), I) && checkbounds(Bool, trailingsize(A,Val{2}), J)) ||
        throw_boundserror(A, (I, J))
end
@generated function _internal_checkbounds(A::AbstractArray, I...)
    meta = Expr(:meta, :inline)
    N = length(I)
    Isplat = [:(I[$d]) for d=1:N]
    error = :(throw_boundserror(A, tuple($(Isplat...))))
    args = Expr[:(checkbounds(Bool, size(A,$dim), I[$dim]) || $error) for dim in 1:N-1]
    push!(args, :(checkbounds(Bool, trailingsize(A,Val{$N}), I[$N]) || $error))
    Expr(:block, meta, args...)
end

## Bounds-checking without errors ##
function checkbounds(::Type{Bool}, sz::Dims, I...)
    n = length(I)
    for dim = 1:(n-1)
        checkbounds(Bool, sz[dim], I[dim]) || return false
    end
    s = sz[n]
    for i = n+1:length(sz)
        s *= sz[i]
    end
    checkbounds(Bool, s, I[n])
end

## Constructors ##

# default arguments to similar()
similar{T}(a::AbstractArray{T})               = similar(a, T, size(a))
similar(   a::AbstractArray, T)               = similar(a, T, size(a))
similar{T}(a::AbstractArray{T}, dims::Dims)   = similar(a, T, dims)
similar{T}(a::AbstractArray{T}, dims::Int...) = similar(a, T, dims)
similar(   a::AbstractArray, T, dims::Int...) = similar(a, T, dims)
# similar creates an Array by default
similar(   a::AbstractArray, T, dims::Dims)   = Array(T, dims)

function reshape(a::AbstractArray, dims::Dims)
    if prod(dims) != length(a)
        throw(ArgumentError("dimensions must be consistent with array size (expected $(length(a)), got $(prod(dims)))"))
    end
    copy!(similar(a, dims), a)
end
reshape(a::AbstractArray, dims::Int...) = reshape(a, dims)

## from general iterable to any array

function copy!(dest::AbstractArray, src)
    i = 1
    for x in src
        dest[i] = x
        i += 1
    end
    return dest
end

# if src is not an AbstractArray, moving to the offset might be O(n)
function copy!(dest::AbstractArray, doffs::Integer, src)
    doffs < 1 && throw(BoundsError(dest, doffs))
    st = start(src)
    i, dmax = doffs, length(dest)
    @inbounds while !done(src, st)
        i > dmax && throw(BoundsError(dest, i))
        val, st = next(src, st)
        dest[i] = val
        i += 1
    end
    return dest
end

# copy from an some iterable object into an AbstractArray
function copy!(dest::AbstractArray, doffs::Integer, src, soffs::Integer)
    if (doffs < 1) | (soffs < 1)
        doffs < 1 && throw(BoundsError(dest, doffs))
        throw(ArgumentError(string("source start offset (",soffs,") is < 1")))
    end
    st = start(src)
    for j = 1:(soffs-1)
        if done(src, st)
            throw(ArgumentError(string("source has fewer elements than required, ",
                                       "expected at least ",soffs,", got ",j-1)))
        end
        _, st = next(src, st)
    end
    dn = done(src, st)
    if dn
        throw(ArgumentError(string("source has fewer elements than required, ",
                                      "expected at least ",soffs,", got ",soffs-1)))
    end
    i, dmax = doffs, length(dest)
    @inbounds while !dn
        i > dmax && throw(BoundsError(dest, i))
        val, st = next(src, st)
        dest[i] = val
        i += 1
        dn = done(src, st)
    end
    return dest
end

# this method must be separate from the above since src might not have a length
function copy!(dest::AbstractArray, doffs::Integer, src, soffs::Integer, n::Integer)
    n < 0 && throw(BoundsError(dest, n))
    n == 0 && return dest
    dmax = doffs + n - 1
    if (dmax > length(dest)) | (doffs < 1) | (soffs < 1)
        doffs < 1 && throw(BoundsError(dest, doffs))
        soffs < 1 && throw(ArgumentError(string("source start offset (",soffs,") is < 1")))
        throw(BoundsError(dest, dmax))
    end
    st = start(src)
    for j = 1:(soffs-1)
        if done(src, st)
            throw(ArgumentError(string("source has fewer elements than required, ",
                                       "expected at least ",soffs,", got ",j-1)))
        end
        _, st = next(src, st)
    end
    i = doffs
    @inbounds while i <= dmax && !done(src, st)
        val, st = next(src, st)
        dest[i] = val
        i += 1
    end
    i <= dmax && throw(BoundsError(dest, i))
    return dest
end

## copy between abstract arrays - generally more efficient
## since a single index variable can be used.

copy!(dest::AbstractArray, src::AbstractArray) =
    copy!(linearindexing(dest), dest, linearindexing(src), src)

function copy!(::LinearIndexing, dest::AbstractArray, ::LinearIndexing, src::AbstractArray)
    n = length(src)
    n > length(dest) && throw(BoundsError(dest, n))
    @inbounds for i = 1:n
        dest[i] = src[i]
    end
    return dest
end

function copy!(::LinearIndexing, dest::AbstractArray, ::LinearSlow, src::AbstractArray)
    n = length(src)
    n > length(dest) && throw(BoundsError(dest, n))
    i = 0
    @inbounds for a in src
        dest[i+=1] = a
    end
    return dest
end

function copy!(dest::AbstractArray, doffs::Integer, src::AbstractArray)
    copy!(dest, doffs, src, 1, length(src))
end

function copy!(dest::AbstractArray, doffs::Integer, src::AbstractArray, soffs::Integer)
    soffs > length(src) && throw(BoundsError(src, soffs))
    copy!(dest, doffs, src, soffs, length(src)-soffs+1)
end

function copy!(dest::AbstractArray, doffs::Integer,
               src::AbstractArray, soffs::Integer,
               n::Integer)
    n == 0 && return dest
    n < 0  && throw(BoundsError(src, n))
    soffs+n-1 > length(src)  && throw(BoundsError(src, soffs+n-1))
    doffs+n-1 > length(dest) && throw(BoundsError(dest, doffs+n-1))
    doffs < 1 && throw(BoundsError(dest, doffs))
    soffs < 1 && throw(BoundsError(src, soffs))
    @inbounds for i = 0:(n-1)
        dest[doffs+i] = src[soffs+i]
    end
    return dest
end

copy(a::AbstractArray) = copy!(similar(a), a)

function copy!{R,S}(B::AbstractVecOrMat{R}, ir_dest::Range{Int}, jr_dest::Range{Int},
                    A::AbstractVecOrMat{S}, ir_src::Range{Int}, jr_src::Range{Int})
    if length(ir_dest) != length(ir_src)
        throw(ArgumentError(string("source and destination must have same size (got ",
                                   length(ir_src)," and ",length(ir_dest),")")))
    end
    if length(jr_dest) != length(jr_src)
        throw(ArgumentError(string("source and destination must have same size (got ",
                                   length(jr_src)," and ",length(jr_dest),")")))
    end
    checkbounds(B, ir_dest, jr_dest)
    checkbounds(A, ir_src, jr_src)
    jdest = first(jr_dest)
    for jsrc in jr_src
        idest = first(ir_dest)
        for isrc in ir_src
            B[idest,jdest] = A[isrc,jsrc]
            idest += step(ir_dest)
        end
        jdest += step(jr_dest)
    end
    return B
end

function copy_transpose!{R,S}(B::AbstractVecOrMat{R}, ir_dest::Range{Int}, jr_dest::Range{Int},
                              A::AbstractVecOrMat{S}, ir_src::Range{Int}, jr_src::Range{Int})
    if length(ir_dest) != length(jr_src)
        throw(ArgumentError(string("source and destination must have same size (got ",
                                   length(jr_src)," and ",length(ir_dest),")")))
    end
    if length(jr_dest) != length(ir_src)
        throw(ArgumentError(string("source and destination must have same size (got ",
                                   length(ir_src)," and ",length(jr_dest),")")))
    end
    checkbounds(B, ir_dest, jr_dest)
    checkbounds(A, ir_src, jr_src)
    idest = first(ir_dest)
    for jsrc in jr_src
        jdest = first(jr_dest)
        for isrc in ir_src
            B[idest,jdest] = A[isrc,jsrc]
            jdest += step(jr_dest)
        end
        idest += step(ir_dest)
    end
    return B
end

zero{T}(x::AbstractArray{T}) = fill!(similar(x), zero(T))

## iteration support for arrays by iterating over `eachindex` in the array ##
# Allows fast iteration by default for both LinearFast and LinearSlow arrays

# While the definitions for LinearFast are all simple enough to inline on their
# own, LinearSlow's CartesianRange is more complicated and requires explicit
# inlining.
start(A::AbstractArray) = (@_inline_meta(); itr = eachindex(A); (itr, start(itr)))
next(A::AbstractArray,i) = (@_inline_meta(); (idx, s) = next(i[1], i[2]); (A[idx], (i[1], s)))
done(A::AbstractArray,i) = done(i[1], i[2])

iterstate(i) = i
iterstate(i::Tuple{UnitRange{Int},Int}) = i[2]

# eachindex iterates over all indices. LinearSlow definitions are later.
eachindex(A::AbstractArray) = (@_inline_meta(); eachindex(linearindexing(A), A))

function eachindex(A::AbstractArray, B::AbstractArray)
    @_inline_meta
    eachindex(linearindexing(A,B), A, B)
end
function eachindex(A::AbstractArray, B::AbstractArray...)
    @_inline_meta
    eachindex(linearindexing(A,B...), A, B...)
end
eachindex(::LinearFast, A::AbstractArray) = 1:length(A)
function eachindex(::LinearFast, A::AbstractArray, B::AbstractArray...)
    @_inline_meta
    1:_maxlength(A, B...)
end
_maxlength(A) = length(A)
function _maxlength(A, B, C...)
    @_inline_meta
    max(length(A), _maxlength(B, C...))
end

isempty(a::AbstractArray) = (length(a) == 0)

## Conversions ##

convert{T,N  }(::Type{AbstractArray{T,N}}, A::AbstractArray{T,N}) = A
convert{T,S,N}(::Type{AbstractArray{T,N}}, A::AbstractArray{S,N}) = copy!(similar(A,T), A)
convert{T,S,N}(::Type{AbstractArray{T  }}, A::AbstractArray{S,N}) = convert(AbstractArray{T,N}, A)

convert{T,N}(::Type{Array}, A::AbstractArray{T,N}) = convert(Array{T,N}, A)

full(x::AbstractArray) = x

map(::Type{Integer},  a::Array) = map!(Integer, similar(a,typeof(Integer(one(eltype(a))))), a)
map(::Type{Signed},   a::Array) = map!(Signed, similar(a,typeof(Signed(one(eltype(a))))), a)
map(::Type{Unsigned}, a::Array) = map!(Unsigned, similar(a,typeof(Unsigned(one(eltype(a))))), a)

## range conversions ##

map{T<:Real}(::Type{T}, r::StepRange) = T(r.start):T(r.step):T(last(r))
map{T<:Real}(::Type{T}, r::UnitRange) = T(r.start):T(last(r))
map{T<:AbstractFloat}(::Type{T}, r::FloatRange) = FloatRange(T(r.start), T(r.step), r.len, T(r.divisor))
function map{T<:AbstractFloat}(::Type{T}, r::LinSpace)
    new_len = T(r.len)
    new_len == r.len || error("$r: too long for $T")
    LinSpace(T(r.start), T(r.stop), new_len, T(r.divisor))
end

## unsafe/pointer conversions ##

# note: the following type definitions don't mean any AbstractArray is convertible to
# a data Ref. they just map the array element type to the pointer type for
# convenience in cases that work.
pointer{T}(x::AbstractArray{T}) = unsafe_convert(Ptr{T}, x)
pointer{T}(x::AbstractArray{T}, i::Integer) = (@_inline_meta; unsafe_convert(Ptr{T},x) + (i-1)*elsize(x))

## Approach:
# We only define one fallback method on getindex for all argument types.
# That dispatches to an (inlined) internal _getindex function, where the goal is
# to transform the indices such that we can call the only getindex method that
# we require AbstractArray subtypes must define, either:
#       getindex(::T, ::Int) # if linearindexing(T) == LinearFast()
#       getindex(::T, ::Int, ::Int, #=...ndims(A) indices...=#) if LinearSlow()
# Unfortunately, it is currently impossible to express the latter method for
# arbitrary dimensionalities. We could get around that with ::CartesianIndex{N},
# but that isn't as obvious and would require that the function be inlined to
# avoid allocations.  If the subtype hasn't defined those methods, it goes back
# to the _getindex function where an error is thrown to prevent stack overflows.
#
# We use the same scheme for unsafe_getindex, with the exception that we can
# fallback to the safe version if the subtype hasn't defined the required
# unsafe method.

function getindex(A::AbstractArray, I...)
    @_inline_meta
    _getindex(linearindexing(A), A, I...)
end
function unsafe_getindex(A::AbstractArray, I...)
    @_inline_meta
    _unsafe_getindex(linearindexing(A), A, I...)
end
## Internal defitions
# 0-dimensional indexing is defined to prevent ambiguities. LinearFast is easy:
_getindex(::LinearFast, A::AbstractArray) = (@_inline_meta; getindex(A, 1))
# But LinearSlow must take into account the dimensionality of the array:
_getindex{T}(::LinearSlow, A::AbstractArray{T,0}) = error("indexing not defined for ", typeof(A))
_getindex(::LinearSlow, A::AbstractVector) = (@_inline_meta; getindex(A, 1))
_getindex(l::LinearSlow, A::AbstractArray) = (@_inline_meta; _getindex(l, A, 1))
_unsafe_getindex(::LinearFast, A::AbstractArray) = (@_inline_meta; unsafe_getindex(A, 1))
_unsafe_getindex{T}(::LinearSlow, A::AbstractArray{T,0}) = error("indexing not defined for ", typeof(A))
_unsafe_getindex(::LinearSlow, A::AbstractVector) = (@_inline_meta; unsafe_getindex(A, 1))
_unsafe_getindex(l::LinearSlow, A::AbstractArray) = (@_inline_meta; _unsafe_getindex(l, A, 1))

_getindex(::LinearIndexing, A::AbstractArray, I...) = error("indexing $(typeof(A)) with types $(typeof(I)) is not supported")
_unsafe_getindex(::LinearIndexing, A::AbstractArray, I...) = (@_inline_meta; getindex(A, I...))

## LinearFast Scalar indexing
_getindex(::LinearFast, A::AbstractArray, I::Int) = error("indexing not defined for ", typeof(A))
function _getindex(::LinearFast, A::AbstractArray, I::Real...)
    @_inline_meta
    # We must check bounds for sub2ind; so we can then call unsafe_getindex
    J = to_indexes(I...)
    checkbounds(A, J...)
    unsafe_getindex(A, sub2ind(size(A), J...))
end
_unsafe_getindex(::LinearFast, A::AbstractArray, I::Real) = (@_inline_meta; getindex(A, I))
function _unsafe_getindex(::LinearFast, A::AbstractArray, I::Real...)
    @_inline_meta
    unsafe_getindex(A, sub2ind(size(A), to_indexes(I...)...))
end

# LinearSlow Scalar indexing
@generated function _getindex{T,AN}(::LinearSlow, A::AbstractArray{T,AN}, I::Real...)
    N = length(I)
    if N == AN
        if all(x->x===Int, I)
            :(error("indexing not defined for ", typeof(A)))
        else
            :($(Expr(:meta, :inline)); getindex(A, to_indexes(I...)...))
        end
    elseif N > AN
        # Drop trailing ones
        Isplat = Expr[:(I[$d]) for d = 1:AN]
        Osplat = Expr[:(to_index(I[$d]) == 1) for d = AN+1:N]
        quote
            $(Expr(:meta, :inline))
            (&)($(Osplat...)) || throw_boundserror(A, I)
            getindex(A, $(Isplat...))
        end
    else
        # Expand the last index into the appropriate number of indices
        Isplat = Expr[:(I[$d]) for d = 1:N-1]
        i = 0
        for d=N:AN
            push!(Isplat, :(s[$(i+=1)]))
        end
        sz = Expr(:tuple)
        sz.args = Expr[:(size(A, $d)) for d=N:AN]
        szcheck = Expr[:(size(A, $d) > 0) for d=N:AN]
        quote
            $(Expr(:meta, :inline))
            # ind2sub requires all dimensions to be > 0:
            (&)($(szcheck...)) || throw_boundserror(A, I)
            s = ind2sub($sz, to_index(I[$N]))
            getindex(A, $(Isplat...))
        end
    end
end
@generated function _unsafe_getindex{T,AN}(::LinearSlow, A::AbstractArray{T,AN}, I::Real...)
    N = length(I)
    if N == AN
        :($(Expr(:meta, :inline)); getindex(A, I...))
    elseif N > AN
        # Drop trailing dimensions (unchecked)
        Isplat = Expr[:(I[$d]) for d = 1:AN]
        quote
            $(Expr(:meta, :inline))
            unsafe_getindex(A, $(Isplat...))
        end
    else
        # Expand the last index into the appropriate number of indices
        Isplat = Expr[:(I[$d]) for d = 1:N-1]
        for d=N:AN
            push!(Isplat, :(s[$(d-N+1)]))
        end
        sz = Expr(:tuple)
        sz.args = Expr[:(size(A, $d)) for d=N:AN]
        quote
            $(Expr(:meta, :inline))
            s = ind2sub($sz, to_index(I[$N]))
            unsafe_getindex(A, $(Isplat...))
        end
    end
end

## Setindex! is defined similarly. We first dispatch to an internal _setindex!
# function that allows dispatch on array storage
function setindex!(A::AbstractArray, v, I...)
    @_inline_meta
    _setindex!(linearindexing(A), A, v, I...)
end
function unsafe_setindex!(A::AbstractArray, v, I...)
    @_inline_meta
    _unsafe_setindex!(linearindexing(A), A, v, I...)
end
## Internal defitions
_setindex!(::LinearFast, A::AbstractArray, v) = (@_inline_meta; setindex!(A, v, 1))
_setindex!{T}(::LinearSlow, A::AbstractArray{T,0}, v) = error("indexing not defined for ", typeof(A))
_setindex!(::LinearSlow, A::AbstractVector, v) = (@_inline_meta; setindex!(A, v, 1))
_setindex!(l::LinearSlow, A::AbstractArray, v) = (@_inline_meta; _setindex!(l, A, v, 1))
_unsafe_setindex!(::LinearFast, A::AbstractArray, v) = (@_inline_meta; unsafe_setindex!(A, v, 1))
_unsafe_setindex!{T}(::LinearSlow, A::AbstractArray{T,0}, v) = error("indexing not defined for ", typeof(A))
_unsafe_setindex!(::LinearSlow, A::AbstractVector, v) = (@_inline_meta; unsafe_setindex!(A, v, 1))
_unsafe_setindex!(l::LinearSlow, A::AbstractArray, v) = (@_inline_meta; _unsafe_setindex!(l, A, v, 1))

_setindex!(::LinearIndexing, A::AbstractArray, v, I...) = error("indexing $(typeof(A)) with types $(typeof(I)) is not supported")
_unsafe_setindex!(::LinearIndexing, A::AbstractArray, v, I...) = (@_inline_meta; setindex!(A, v, I...))

## LinearFast Scalar indexing
_setindex!(::LinearFast, A::AbstractArray, v, I::Int) = error("indexed assignment not defined for ", typeof(A))
function _setindex!(::LinearFast, A::AbstractArray, v, I::Real...)
    @_inline_meta
    # We must check bounds for sub2ind; so we can then call unsafe_setindex!
    J = to_indexes(I...)
    checkbounds(A, J...)
    unsafe_setindex!(A, v, sub2ind(size(A), J...))
end
_unsafe_setindex!(::LinearFast, A::AbstractArray, v, I::Real) = (@_inline_meta; setindex!(A, v, I))
function _unsafe_setindex!(::LinearFast, A::AbstractArray, v, I::Real...)
    @_inline_meta
    unsafe_setindex!(A, v, sub2ind(size(A), to_indexes(I...)...))
end

# LinearSlow Scalar indexing
@generated function _setindex!{T,AN}(::LinearSlow, A::AbstractArray{T,AN}, v, I::Real...)
    N = length(I)
    if N == AN
        if all(x->x===Int, I)
            :(error("indexing not defined for ", typeof(A)))
        else
            :($(Expr(:meta, :inline)); setindex!(A, v, to_indexes(I...)...))
        end
    elseif N > AN
        # Drop trailing ones
        Isplat = Expr[:(I[$d]) for d = 1:AN]
        Osplat = Expr[:(to_index(I[$d]) == 1) for d = AN+1:N]
        quote
            $(Expr(:meta, :inline))
            (&)($(Osplat...)) || throw_boundserror(A, I)
            setindex!(A, v, $(Isplat...))
        end
    else
        # Expand the last index into the appropriate number of indices
        Isplat = Expr[:(I[$d]) for d = 1:N-1]
        i = 0
        for d=N:AN
            push!(Isplat, :(s[$(i+=1)]))
        end
        sz = Expr(:tuple)
        sz.args = Expr[:(size(A, $d)) for d=N:AN]
        szcheck = Expr[:(size(A, $d) > 0) for d=N:AN]
        quote
            $(Expr(:meta, :inline))
            # ind2sub requires all dimensions to be > 0:
            (&)($(szcheck...)) || throw_boundserror(A, I)
            s = ind2sub($sz, to_index(I[$N]))
            setindex!(A, v, $(Isplat...))
        end
    end
end
@generated function _unsafe_setindex!{T,AN}(::LinearSlow, A::AbstractArray{T,AN}, v, I::Real...)
    N = length(I)
    if N == AN
        :($(Expr(:meta, :inline)); setindex!(A, v, I...))
    elseif N > AN
        # Drop trailing dimensions (unchecked)
        Isplat = Expr[:(I[$d]) for d = 1:AN]
        quote
            $(Expr(:meta, :inline))
            unsafe_setindex!(A, v, $(Isplat...))
        end
    else
        # Expand the last index into the appropriate number of indices
        Isplat = Expr[:(I[$d]) for d = 1:N-1]
        for d=N:AN
            push!(Isplat, :(s[$(d-N+1)]))
        end
        sz = Expr(:tuple)
        sz.args = Expr[:(size(A, $d)) for d=N:AN]
        quote
            $(Expr(:meta, :inline))
            s = ind2sub($sz, to_index(I[$N]))
            unsafe_setindex!(A, v, $(Isplat...))
        end
    end
end

## get (getindex with a default value) ##

typealias RangeVecIntList{A<:AbstractVector{Int}} Union{Tuple{Vararg{Union{Range, AbstractVector{Int}}}}, AbstractVector{UnitRange{Int}}, AbstractVector{Range{Int}}, AbstractVector{A}}

get(A::AbstractArray, i::Integer, default) = checkbounds(Bool, length(A), i) ? A[i] : default
get(A::AbstractArray, I::Tuple{}, default) = similar(A, typeof(default), 0)
get(A::AbstractArray, I::Dims, default) = checkbounds(Bool, size(A), I...) ? A[I...] : default

function get!{T}(X::AbstractArray{T}, A::AbstractArray, I::Union{Range, AbstractVector{Int}}, default::T)
    ind = findin(I, 1:length(A))
    X[ind] = A[I[ind]]
    X[1:first(ind)-1] = default
    X[last(ind)+1:length(X)] = default
    X
end

get(A::AbstractArray, I::Range, default) = get!(similar(A, typeof(default), length(I)), A, I, default)

function get!{T}(X::AbstractArray{T}, A::AbstractArray, I::RangeVecIntList, default::T)
    fill!(X, default)
    dst, src = indcopy(size(A), I)
    X[dst...] = A[src...]
    X
end

get(A::AbstractArray, I::RangeVecIntList, default) = get!(similar(A, typeof(default), map(length, I)...), A, I, default)

## Concatenation ##

promote_eltype() = Bottom
promote_eltype(v1, vs...) = promote_type(eltype(v1), promote_eltype(vs...))

#TODO: ERROR CHECK
cat(catdim::Integer) = Array(Any, 0)

vcat() = Array(Any, 0)
hcat() = Array(Any, 0)
typed_vcat(T::Type) = Array(T, 0)
typed_hcat(T::Type) = Array(T, 0)

## cat: special cases
vcat{T}(X::T...)         = T[ X[i] for i=1:length(X) ]
vcat{T<:Number}(X::T...) = T[ X[i] for i=1:length(X) ]
hcat{T}(X::T...)         = T[ X[j] for i=1, j=1:length(X) ]
hcat{T<:Number}(X::T...) = T[ X[j] for i=1, j=1:length(X) ]

vcat(X::Number...) = hvcat_fill(Array(promote_typeof(X...),length(X)), X)
hcat(X::Number...) = hvcat_fill(Array(promote_typeof(X...),1,length(X)), X)
typed_vcat(T::Type, X::Number...) = hvcat_fill(Array(T,length(X)), X)
typed_hcat(T::Type, X::Number...) = hvcat_fill(Array(T,1,length(X)), X)

vcat(V::AbstractVector...) = typed_vcat(promote_eltype(V...), V...)
vcat{T}(V::AbstractVector{T}...) = typed_vcat(T, V...)

function typed_vcat(T::Type, V::AbstractVector...)
    n::Int = 0
    for Vk in V
        n += length(Vk)
    end
    a = similar(full(V[1]), T, n)
    pos = 1
    for k=1:length(V)
        Vk = V[k]
        p1 = pos+length(Vk)-1
        a[pos:p1] = Vk
        pos = p1+1
    end
    a
end

hcat(A::AbstractVecOrMat...) = typed_hcat(promote_eltype(A...), A...)
hcat{T}(A::AbstractVecOrMat{T}...) = typed_hcat(T, A...)

function typed_hcat(T::Type, A::AbstractVecOrMat...)
    nargs = length(A)
    nrows = size(A[1], 1)
    ncols = 0
    dense = true
    for j = 1:nargs
        Aj = A[j]
        if size(Aj, 1) != nrows
            throw(ArgumentError("number of rows of each array must match (got $(map(x->size(x,1), A)))"))
        end
        dense &= isa(Aj,Array)
        nd = ndims(Aj)
        ncols += (nd==2 ? size(Aj,2) : 1)
    end
    B = similar(full(A[1]), T, nrows, ncols)
    pos = 1
    if dense
        for k=1:nargs
            Ak = A[k]
            n = length(Ak)
            copy!(B, pos, Ak, 1, n)
            pos += n
        end
    else
        for k=1:nargs
            Ak = A[k]
            p1 = pos+(isa(Ak,AbstractMatrix) ? size(Ak, 2) : 1)-1
            B[:, pos:p1] = Ak
            pos = p1+1
        end
    end
    return B
end

vcat(A::AbstractMatrix...) = typed_vcat(promote_eltype(A...), A...)
vcat{T}(A::AbstractMatrix{T}...) = typed_vcat(T, A...)

function typed_vcat(T::Type, A::AbstractMatrix...)
    nargs = length(A)
    nrows = sum(a->size(a, 1), A)::Int
    ncols = size(A[1], 2)
    for j = 2:nargs
        if size(A[j], 2) != ncols
            throw(ArgumentError("number of columns of each array must match (got $(map(x->size(x,2), A)))"))
        end
    end
    B = similar(full(A[1]), T, nrows, ncols)
    pos = 1
    for k=1:nargs
        Ak = A[k]
        p1 = pos+size(Ak,1)-1
        B[pos:p1, :] = Ak
        pos = p1+1
    end
    return B
end

## cat: general case

function cat(catdims, X...)
    T = promote_type(map(x->isa(x,AbstractArray) ? eltype(x) : typeof(x), X)...)
    cat_t(catdims, T, X...)
end

function cat_t(catdims, typeC::Type, X...)
    catdims = collect(catdims)
    nargs = length(X)
    ndimsX = Int[isa(a,AbstractArray) ? ndims(a) : 0 for a in X]
    ndimsC = max(maximum(ndimsX), maximum(catdims))
    catsizes = zeros(Int,(nargs,length(catdims)))
    dims2cat = zeros(Int,ndimsC)
    for k = 1:length(catdims)
        dims2cat[catdims[k]]=k
    end

    dimsC = Int[d <= ndimsX[1] ? size(X[1],d) : 1 for d=1:ndimsC]
    for k = 1:length(catdims)
        catsizes[1,k] = dimsC[catdims[k]]
    end
    for i = 2:nargs
        for d = 1:ndimsC
            currentdim = (d <= ndimsX[i] ? size(X[i],d) : 1)
            if dims2cat[d] != 0
                dimsC[d] += currentdim
                catsizes[i,dims2cat[d]] = currentdim
            elseif dimsC[d] != currentdim
                throw(DimensionMismatch(string("mismatch in dimension ",d,
                                               " (expected ",dimsC[d],
                                               " got ",currentdim,")")))
            end
        end
    end

    C = similar(isa(X[1],AbstractArray) ? full(X[1]) : [X[1]], typeC, tuple(dimsC...))
    if length(catdims)>1
        fill!(C,0)
    end

    offsets = zeros(Int,length(catdims))
    for i=1:nargs
        cat_one = [ dims2cat[d] == 0 ? (1:dimsC[d]) : (offsets[dims2cat[d]]+(1:catsizes[i,dims2cat[d]]))
                   for d=1:ndimsC ]
        C[cat_one...] = X[i]
        for k = 1:length(catdims)
            offsets[k] += catsizes[i,k]
        end
    end
    return C
end

vcat(X...) = cat(1, X...)
hcat(X...) = cat(2, X...)

typed_vcat(T::Type, X...) = cat_t(1, T, X...)
typed_hcat(T::Type, X...) = cat_t(2, T, X...)

cat{T}(catdims, A::AbstractArray{T}...) = cat_t(catdims, T, A...)

cat(catdims, A::AbstractArray...) = cat_t(catdims, promote_eltype(A...), A...)

# The specializations for 1 and 2 inputs are important
# especially when running with --inline=no, see #11158
vcat(A::AbstractArray) = cat(1, A)
vcat(A::AbstractArray, B::AbstractArray) = cat(1, A, B)
vcat(A::AbstractArray...) = cat(1, A...)
hcat(A::AbstractArray) = cat(2, A)
hcat(A::AbstractArray, B::AbstractArray) = cat(2, A, B)
hcat(A::AbstractArray...) = cat(2, A...)

typed_vcat(T::Type, A::AbstractArray) = cat_t(1, T, A)
typed_vcat(T::Type, A::AbstractArray, B::AbstractArray) = cat_t(1, T, A, B)
typed_vcat(T::Type, A::AbstractArray...) = cat_t(1, T, A...)
typed_hcat(T::Type, A::AbstractArray) = cat_t(2, T, A)
typed_hcat(T::Type, A::AbstractArray, B::AbstractArray) = cat_t(2, T, A, B)
typed_hcat(T::Type, A::AbstractArray...) = cat_t(2, T, A...)

# 2d horizontal and vertical concatenation

function hvcat(nbc::Integer, as...)
    # nbc = # of block columns
    n = length(as)
    mod(n,nbc) != 0 &&
        throw(ArgumentError("number of arrays $n is not a multiple of the requested number of block columns $nbc"))
    nbr = div(n,nbc)
    hvcat(ntuple(i->nbc, nbr), as...)
end

function hvcat{T}(rows::Tuple{Vararg{Int}}, as::AbstractMatrix{T}...)
    nbr = length(rows)  # number of block rows

    nc = 0
    for i=1:rows[1]
        nc += size(as[i],2)
    end

    nr = 0
    a = 1
    for i = 1:nbr
        nr += size(as[a],1)
        a += rows[i]
    end

    out = similar(full(as[1]), T, nr, nc)

    a = 1
    r = 1
    for i = 1:nbr
        c = 1
        szi = size(as[a],1)
        for j = 1:rows[i]
            Aj = as[a+j-1]
            szj = size(Aj,2)
            if size(Aj,1) != szi
                throw(ArgumentError("mismatched height in block row $(i) (expected $szi, got $(size(Aj,1)))"))
            end
            if c-1+szj > nc
                throw(ArgumentError("block row $(i) has mismatched number of columns (expected $nc, got $(c-1+szj))"))
            end
            out[r:r-1+szi, c:c-1+szj] = Aj
            c += szj
        end
        if c != nc+1
            throw(ArgumentError("block row $(i) has mismatched number of columns (expected $nc, got $(c-1))"))
        end
        r += szi
        a += rows[i]
    end
    out
end

hvcat(rows::Tuple{Vararg{Int}}) = []

function hvcat{T<:Number}(rows::Tuple{Vararg{Int}}, xs::T...)
    nr = length(rows)
    nc = rows[1]

    a = Array(T, nr, nc)
    if length(a) != length(xs)
        throw(ArgumentError("argument count does not match specified shape (expected $(length(a)), got $(length(xs)))"))
    end
    k = 1
    @inbounds for i=1:nr
        if nc != rows[i]
            throw(ArgumentError("row $(i) has mismatched number of columns (expected $nc, got $(rows[i]))"))
        end
        for j=1:nc
            a[i,j] = xs[k]
            k += 1
        end
    end
    a
end

function hvcat_fill(a, xs)
    k = 1
    nr, nc = size(a,1), size(a,2)
    for i=1:nr
        @inbounds for j=1:nc
            a[i,j] = xs[k]
            k += 1
        end
    end
    a
end

function typed_hvcat(T::Type, rows::Tuple{Vararg{Int}}, xs::Number...)
    nr = length(rows)
    nc = rows[1]
    for i = 2:nr
        if nc != rows[i]
            throw(ArgumentError("row $(i) has mismatched number of columns (expected $nc, got $(rows[i]))"))
        end
    end
    len = length(xs)
    if nr*nc != len
        throw(ArgumentError("argument count $(len) does not match specified shape $((nr,nc))"))
    end
    hvcat_fill(Array(T, nr, nc), xs)
end

function hvcat(rows::Tuple{Vararg{Int}}, xs::Number...)
    T = promote_typeof(xs...)
    typed_hvcat(T, rows, xs...)
end

# fallback definition of hvcat in terms of hcat and vcat
function hvcat(rows::Tuple{Vararg{Int}}, as...)
    nbr = length(rows)  # number of block rows
    rs = cell(nbr)
    a = 1
    for i = 1:nbr
        rs[i] = hcat(as[a:a-1+rows[i]]...)
        a += rows[i]
    end
    vcat(rs...)
end

function typed_hvcat(T::Type, rows::Tuple{Vararg{Int}}, as...)
    nbr = length(rows)  # number of block rows
    rs = cell(nbr)
    a = 1
    for i = 1:nbr
        rs[i] = hcat(as[a:a-1+rows[i]]...)
        a += rows[i]
    end
    T[rs...;]
end

## Reductions and scans ##

function isequal(A::AbstractArray, B::AbstractArray)
    if A === B return true end
    if size(A) != size(B)
        return false
    end
    if isa(A,Range) != isa(B,Range)
        return false
    end
    for i in eachindex(A,B)
        if !isequal(A[i], B[i])
            return false
        end
    end
    return true
end

function lexcmp(A::AbstractArray, B::AbstractArray)
    nA, nB = length(A), length(B)
    for i = 1:min(nA, nB)
        res = lexcmp(A[i], B[i])
        res == 0 || return res
    end
    return cmp(nA, nB)
end

function (==)(A::AbstractArray, B::AbstractArray)
    if size(A) != size(B)
        return false
    end
    if isa(A,Range) != isa(B,Range)
        return false
    end
    for i in eachindex(A,B)
        if !(A[i]==B[i])
            return false
        end
    end
    return true
end

sub2ind(dims::Tuple{Vararg{Integer}}) = 1
sub2ind(dims::Tuple{Vararg{Integer}}, I::Integer...) = _sub2ind(dims,I)
@generated function _sub2ind{N,M}(dims::NTuple{N,Integer}, I::NTuple{M,Integer})
    meta = Expr(:meta,:inline)
    ex = :(I[$M] - 1)
    for i = M-1:-1:1
        if i > N
            ex = :(I[$i] - 1 + $ex)
        else
            ex = :(I[$i] - 1 + dims[$i]*$ex)
        end
    end
    Expr(:block, meta,:($ex + 1))
end

@generated function ind2sub{N}(dims::NTuple{N,Integer}, ind::Integer)
    meta = Expr(:meta,:inline)
    N==0 && return :($meta; ind==1 ? () : throw(BoundsError()))
    exprs = Expr[:(ind = ind-1)]
    for i = 1:N-1
        push!(exprs,:(ind2 = div(ind,dims[$i])))
        push!(exprs,Expr(:(=),symbol(:s,i),:(ind-dims[$i]*ind2+1)))
        push!(exprs,:(ind=ind2))
    end
    push!(exprs,Expr(:(=),symbol(:s,N),:(ind+1)))
    Expr(:block,meta,exprs...,Expr(:tuple,[symbol(:s,i) for i=1:N]...))
end

# TODO in v0.5: either deprecate line 1 or add line 2
ind2sub(a::AbstractArray, ind::Integer) = ind2sub(size(a), ind)
# sub2ind(a::AbstractArray, I::Integer...) = sub2ind(size(a), I...)

function sub2ind{T<:Integer}(dims::Tuple{Vararg{Integer}}, I::AbstractVector{T}...)
    N = length(dims)
    M = length(I[1])
    indices = Array{T}(length(I[1]))
    copy!(indices,I[1])

    s = dims[1]
    for j=2:length(I)
        Ij = I[j]
        for i=1:M
            indices[i] += s*(Ij[i]-1)
        end
        s *= (j <= N ? dims[j] : 1)
    end
    return indices
end

function ind2sub{N,T<:Integer}(dims::NTuple{N,Integer}, ind::AbstractVector{T})
    M = length(ind)
    t = NTuple{N,Vector{T}}(ntuple(n->Array{T}(M),N))
    copy!(t[1],ind)
    for j = 1:N-1
        d = dims[j]
        tj = t[j]
        tj2 = t[j+1]
        for i = 1:M
            ind2 = div(tj[i]-1, d)
            tj[i] -= d*ind2
            tj2[i] = ind2+1
        end
    end
    return t
end

function ind2sub!{T<:Integer}(sub::Array{T}, dims::Tuple{Vararg{T}}, ind::T)
    ndims = length(dims)
    for i=1:ndims-1
        ind2 = div(ind-1,dims[i])+1
        sub[i] = ind - dims[i]*(ind2-1)
        ind = ind2
    end
    sub[ndims] = ind
    return sub
end

## iteration utilities ##

# generic map on any iterator
function map(f, iters...)
    result = []
    len = length(iters)
    states = [start(iters[idx]) for idx in 1:len]
    nxtvals = cell(len)
    cont = true
    for idx in 1:len
        if done(iters[idx], states[idx])
            cont = false
            break
        end
    end
    while cont
        for idx in 1:len
            nxtvals[idx],states[idx] = next(iters[idx], states[idx])
        end
        push!(result, f(nxtvals...))
        for idx in 1:len
            if done(iters[idx], states[idx])
                cont = false
                break
            end
        end
    end
    result
end

## map over arrays ##

## transform any set of dimensions
## dims specifies which dimensions will be transformed. for example
## dims==1:2 will call f on all slices A[:,:,...]
mapslices(f, A::AbstractArray, dims) = mapslices(f, A, [dims...])
function mapslices(f, A::AbstractArray, dims::AbstractVector)
    if isempty(dims)
        return map(f,A)
    end

    dimsA = [size(A)...]
    ndimsA = ndims(A)
    alldims = [1:ndimsA;]

    otherdims = setdiff(alldims, dims)

    idx = cell(ndimsA)
    fill!(idx, 1)
    Asliceshape = tuple(dimsA[dims]...)
    itershape   = tuple(dimsA[otherdims]...)
    for d in dims
        idx[d] = 1:size(A,d)
    end

    r1 = f(reshape(A[idx...], Asliceshape))

    # determine result size and allocate
    Rsize = copy(dimsA)
    # TODO: maybe support removing dimensions
    if !isa(r1, AbstractArray) || ndims(r1) == 0
        r1 = [r1]
    end
    Rsize[dims] = [size(r1)...; ones(Int,max(0,length(dims)-ndims(r1)))]
    R = similar(r1, tuple(Rsize...))

    ridx = cell(ndims(R))
    fill!(ridx, 1)
    for d in dims
        ridx[d] = 1:size(R,d)
    end

    R[ridx...] = r1

    first = true
    nidx = length(otherdims)
    for I in CartesianRange(itershape)
        if first
            first = false
        else
            for i in 1:nidx
                idx[otherdims[i]] = ridx[otherdims[i]] = I.I[i]
            end
            R[ridx...] = f(reshape(A[idx...], Asliceshape))
        end
    end

    return R
end

# using promote_type
function promote_to!{T,F}(f::F, offs, dest::AbstractArray{T}, A::AbstractArray)
    # map to dest array, checking the type of each result. if a result does not
    # match, do a type promotion and re-dispatch.
    @inbounds for i = offs:length(A)
        el = f(A[i])
        S = typeof(el)
        if S === T || S <: T
            dest[i] = el::T
        else
            R = promote_type(T, S)
            if R !== T
                new = similar(dest, R)
                copy!(new,1, dest,1, i-1)
                new[i] = el
                return promote_to!(f, i+1, new, A)
            end
            dest[i] = el
        end
    end
    return dest
end

function map_promote(f, A::AbstractArray)
    if isempty(A); return similar(A, Bottom); end
    first = f(A[1])
    dest = similar(A, typeof(first))
    dest[1] = first
    return promote_to!(f, 2, dest, A)
end

## 1 argument
map!{F}(f::F, A::AbstractArray) = map!(f, A, A)
function map!{F}(f::F, dest::AbstractArray, A::AbstractArray)
    for i = 1:length(A)
        dest[i] = f(A[i])
    end
    return dest
end

function map_to!{T,F}(f::F, offs, dest::AbstractArray{T}, A::AbstractArray)
    # map to dest array, checking the type of each result. if a result does not
    # match, widen the result type and re-dispatch.
    @inbounds for i = offs:length(A)
        el = f(A[i])
        S = typeof(el)
        if S === T || S <: T
            dest[i] = el::T
        else
            R = typejoin(T, S)
            new = similar(dest, R)
            copy!(new,1, dest,1, i-1)
            new[i] = el
            return map_to!(f, i+1, new, A)
        end
    end
    return dest
end

function map(f, A::AbstractArray)
    if isempty(A)
        return isa(f,Type) ? similar(A,f) : similar(A)
    end
    first = f(A[1])
    dest = similar(A, typeof(first))
    dest[1] = first
    return map_to!(f, 2, dest, A)
end

## 2 argument
function map!{F}(f::F, dest::AbstractArray, A::AbstractArray, B::AbstractArray)
    for i = 1:length(A)
        dest[i] = f(A[i], B[i])
    end
    return dest
end

function map_to!{T,F}(f::F, offs, dest::AbstractArray{T}, A::AbstractArray, B::AbstractArray)
    @inbounds for i = offs:length(A)
        el = f(A[i], B[i])
        S = typeof(el)
        if (S !== T) && !(S <: T)
            R = typejoin(T, S)
            new = similar(dest, R)
            copy!(new,1, dest,1, i-1)
            new[i] = el
            return map_to!(f, i+1, new, A, B)
        end
        dest[i] = el::T
    end
    return dest
end

function map(f, A::AbstractArray, B::AbstractArray)
    shp = promote_shape(size(A),size(B))
    if prod(shp) == 0
        return similar(A, promote_type(eltype(A),eltype(B)), shp)
    end
    first = f(A[1], B[1])
    dest = similar(A, typeof(first), shp)
    dest[1] = first
    return map_to!(f, 2, dest, A, B)
end

## N argument

ith_all(i, ::Tuple{}) = ()
ith_all(i, as) = (as[1][i], ith_all(i, tail(as))...)

function map_n!{F}(f::F, dest::AbstractArray, As)
    n = length(As[1])
    for i = 1:n
        dest[i] = f(ith_all(i, As)...)
    end
    return dest
end

map!{F}(f::F, dest::AbstractArray, As::AbstractArray...) = map_n!(f, dest, As)

function map_to_n!{T,F}(f::F, offs, dest::AbstractArray{T}, As)
    @inbounds for i = offs:length(As[1])
        el = f(ith_all(i, As)...)
        S = typeof(el)
        if (S !== T) && !(S <: T)
            R = typejoin(T, S)
            new = similar(dest, R)
            copy!(new,1, dest,1, i-1)
            new[i] = el
            return map_to_n!(f, i+1, new, As)
        end
        dest[i] = el::T
    end
    return dest
end

function map(f, As::AbstractArray...)
    shape = mapreduce(size, promote_shape, As)
    if prod(shape) == 0
        return similar(As[1], promote_eltype(As...), shape)
    end
    first = f(map(a->a[1], As)...)
    dest = similar(As[1], typeof(first), shape)
    dest[1] = first
    return map_to_n!(f, 2, dest, As)
end

# multi-item push!, unshift! (built on top of type-specific 1-item version)
# (note: must not cause a dispatch loop when 1-item case is not defined)
push!(A, a, b) = push!(push!(A, a), b)
push!(A, a, b, c...) = push!(push!(A, a, b), c...)
unshift!(A, a, b) = unshift!(unshift!(A, b), a)
unshift!(A, a, b, c...) = unshift!(unshift!(A, c...), a, b)

## hashing collections ##

const hashaa_seed = UInt === UInt64 ? 0x7f53e68ceb575e76 : 0xeb575e76
const hashrle_seed = UInt == UInt64 ? 0x2aab8909bfea414c : 0xbfea414c
function hash(a::AbstractArray, h::UInt)
    h += hashaa_seed
    h += hash(size(a))

    state = start(a)
    done(a, state) && return h
    x2, state = next(a, state)
    done(a, state) && return hash(x2, h)

    x1 = x2
    while !done(a, state)
        x1 = x2
        x2, state = next(a, state)
        if isequal(x2, x1)
            # For repeated elements, use run length encoding
            # This allows efficient hashing of sparse arrays
            runlength = 2
            while !done(a, state)
                x2, state = next(a, state)
                isequal(x1, x2) || break
                runlength += 1
            end
            h += hashrle_seed
            h = hash(runlength, h)
        end
        h = hash(x1, h)
    end
    !isequal(x2, x1) && (h = hash(x2, h))
    return h
end

max_stars_repo_name: lidiapopescu/JuliaDT

max_issues_repo_name: lidiapopescu/JuliaDT

max_forks_repo_name: JuliaComputing/JuliaDT


text:
# # Linear shell
#
# ![](linear_shell.png)
#-
# ## Introduction
#
# In this example we show how shell elements can be analyzed in Ferrite.jl. The shell implemented here comes from the book
# "The finite elment method - Linear static and dynamic finite element analysis" by Hughes (1987), and a brief description of it is
# given at the end of this tutorial.  The first part of the tutorial explains how to set up the problem.

# ## Setting up the problem
using Ferrite
using ForwardDiff

function main() #wrap everything in a function...

# First we generate a flat rectangular mesh. There is currently no built-in function for generating
# shell meshes in Ferrite, so we have to create our own simple mesh generator (see the
# function `generate_shell_grid` further down in this file).
#+
nels = (10,10)
size = (10.0, 10.0)
grid = generate_shell_grid(nels, size)

# Here we define the bi-linear interpolation used for the geometrical description of the shell.
# We also create two quadrature rules for the in-plane and out-of-plane directions. Note that we use
# under integration for the inplane integration, to avoid shear locking.
#+
ip = Lagrange{2,RefCube,1}()
qr_inplane = QuadratureRule{2,RefCube}(1)
qr_ooplane = QuadratureRule{1,RefCube}(2)
cv = CellScalarValues(qr_inplane, ip)

# Next we distribute displacement dofs,`:u = (x,y,z)` and rotational dofs, `:θ = (θ₁,  θ₂)`.
#+
dh = DofHandler(grid)
push!(dh, :u, 3, ip)
push!(dh, :θ, 2, ip)
close!(dh)

# In order to apply our boundary conditions, we first need to create some edge- and vertex-sets. This
# is done with `addedgeset!` and `addvertexset!` (similar to `addfaceset!`)
#+
addedgeset!(grid, "left",  (x) -> x[1] ≈ 0.0)
addedgeset!(grid, "right", (x) -> x[1] ≈ size[1])
addvertexset!(grid, "corner", (x) -> x[1] ≈ 0.0 && x[2] ≈ 0.0 && x[3] ≈ 0.0)

# Here we define the boundary conditions. On the left edge, we lock the displacements in the x- and z- directions, and all the rotations.
#+
ch = ConstraintHandler(dh)
add!(ch,  Dirichlet(:u, getedgeset(grid, "left"), (x, t) -> (0.0, 0.0), [1,3])  )
add!(ch,  Dirichlet(:θ, getedgeset(grid, "left"), (x, t) -> (0.0, 0.0), [1,2])  )

# On the right edge, we also lock the displacements in the x- and z- directions, but apply a precribed roation.
#+
add!(ch,  Dirichlet(:u, getedgeset(grid, "right"), (x, t) -> (0.0, 0.0), [1,3])  )
add!(ch,  Dirichlet(:θ, getedgeset(grid, "right"), (x, t) -> (0.0, pi/10), [1,2])  )

# In order to not get rigid body motion, we lock the y-displacement in one fo the corners.
#+
add!(ch,  Dirichlet(:θ, getvertexset(grid, "corner"), (x, t) -> (0.0), [2])  )

close!(ch)
update!(ch, 0.0)

# Next we define relevant data for the shell, such as shear correction factor and stiffness matrix for the material.
# In this linear shell, plane stress is assumed, ie $\\sigma_{zz} = 0 $. Therefor, the stiffness matrix is 5x5 (opposed to the normal 6x6).
#+
κ = 5/6 # Shear correction factor
E = 210.0
ν = 0.3
a = (1-ν)/2
C = E/(1-ν^2) * [1 ν 0   0   0;
                ν 1 0   0   0;
                0 0 a*κ 0   0;
                0 0 0   a*κ 0;
                0 0 0   0   a*κ]

data = (thickness = 1.0, C = C); #Named tuple

# We now assemble the problem in standard finite element fashion
#+
nnodes = getnbasefunctions(ip)
ndofs_shell = ndofs_per_cell(dh)

K = create_sparsity_pattern(dh)
f = zeros(Float64, ndofs(dh))

ke = zeros(ndofs_shell, ndofs_shell)
fe = zeros(ndofs_shell)

celldofs = zeros(Int, ndofs_shell)
cellcoords = zeros(Vec{3,Float64}, nnodes)

assembler = start_assemble(K, f)
for cellid in 1:getncells(grid)
    fill!(ke, 0.0)

    celldofs!(celldofs, dh, cellid)
    getcoordinates!(cellcoords, grid, cellid)

    #Call the element routine
    integrate_shell!(ke, cv, qr_ooplane, cellcoords, data)

    assemble!(assembler, celldofs, fe, ke)
end

# Apply BC and solve.
#+
apply!(K, f, ch)
a = K\f

# Output results.
#+
vtk_grid("linear_shell", dh) do vtk
    vtk_point_data(vtk, dh, a)
end

end; #end main functions

# Below is the function that creates the shell mesh. It simply generates a 2d-quadrature mesh, and appends
# a third coordinate (z-direction) to the node-positions.
function generate_shell_grid(nels, size)
    _grid = generate_grid(Quadrilateral, nels, Vec((0.0,0.0)), Vec(size))
    nodes = [(n.x[1], n.x[2], 0.0) |> Vec{3} |> Node  for n in _grid.nodes]
    cells = [Quadrilateral3D(cell.nodes) for cell in _grid.cells]

    grid = Grid(cells, nodes)

    return grid
end;

# ## The shell element
#
# The shell presented here comes from the book "The finite elment method - Linear static and dynamic finite element analysis" by Hughes (1987).
# The shell is a so called degenerate shell element, meaning it is based on a continuum element.
# A brief describtion of the shell is given here.

#md # !!! note
#md #     This element might experience various locking phenomenas, and should only be seen as a proof of concept.

# ##### Fiber coordinate system
# The element uses two coordinate systems. The first coordianate system, called the fiber system, is created for each
# element node, and is used as a reference frame for the rotations. The function below implements an algorthim that return the
# fiber directions, $\boldsymbol{e}^{f}_{a1}$, $\boldsymbol{e}^{f}_{a2}$ and $\boldsymbol{e}^{f}_{a3}$, at each node $a$.
function fiber_coordsys(Ps::Vector{Vec{3,Float64}})

    ef1 = Vec{3,Float64}[]
    ef2 = Vec{3,Float64}[]
    ef3 = Vec{3,Float64}[]
    for P in Ps
        a = abs.(P)
        j = 1
        if a[1] > a[3]; a[3] = a[1]; j = 2; end
        if a[2] > a[3]; j = 3; end

        e3 = P
        e2 = Tensors.cross(P, basevec(Vec{3}, j))
        e2 /= norm(e2)
        e1 = Tensors.cross(e2, P)

        push!(ef1, e1)
        push!(ef2, e2)
        push!(ef3, e3)
    end
    return ef1, ef2, ef3

end;

# ##### Lamina coordinate system
# The second coordinate system is the so called Lamina Coordinate system. It is
# created for each integration point, and is defined to be tangent to the
# mid-surface. It is in this system that we enforce that plane stress assumption,
# i.e. $\sigma_{zz} = 0$. The function below returns the rotation matrix, $\boldsymbol{q}$, for this coordinate system.
function lamina_coordsys(dNdξ, ζ, x, p, h)

    e1 = zero(Vec{3})
    e2 = zero(Vec{3})

    for i in 1:length(dNdξ)
        e1 += dNdξ[i][1] * x[i] + 0.5*h*ζ * dNdξ[i][1] * p[i]
        e2 += dNdξ[i][2] * x[i] + 0.5*h*ζ * dNdξ[i][1] * p[i]
    end

    e1 /= norm(e1)
    e2 /= norm(e2)

    ez = Tensors.cross(e1,e2)
    ez /= norm(ez)

    a = 0.5*(e1 + e2)
    a /= norm(a)

    b = Tensors.cross(ez,a)
    b /= norm(b)

    ex = sqrt(2)/2 * (a - b)
    ey = sqrt(2)/2 * (a + b)

    return Tensor{2,3}(hcat(ex,ey,ez))
end;

# ##### Geometrical description
# A material point in the shell is defined as
# ```math
# \boldsymbol x(\xi, \eta, \zeta) = \sum_{a=1}^{N_{\text{nodes}}} N_a(\xi, \eta) \boldsymbol{\bar{x}}_{a} + ζ \frac{h}{2} \boldsymbol{\bar{p}_a}
# ```
# where $\boldsymbol{\bar{x}}_{a}$ are nodal positions on the mid-surface, and $\boldsymbol{\bar{p}_a}$ is an vector that defines the fiber direction
# on the reference surface. $N_a$ arethe shape functions.
#
# Based on the defintion of the position vector, we create an function for obtaining the Jacobian-matrix,
# ```math
# J_{ij} = \frac{\partial x_i}{\partial \xi_j},
# ```
function getjacobian(q, N, dNdξ, ζ, X, p, h)

    J = zeros(3,3)
    for a in 1:length(N)
        for i in 1:3, j in 1:3
            _dNdξ = (j==3) ? 0.0 : dNdξ[a][j]
            _dζdξ = (j==3) ? 1.0 : 0.0
            _N = N[a]

            J[i,j] += _dNdξ * X[a][i]  +  (_dNdξ*ζ + _N*_dζdξ) * h/2 * p[a][i]
        end
    end

    return (q' * J) |> Tensor{2,3,Float64}
end;

# ##### Strains
# Small deformation is assumed,
# ```math
# \varepsilon_{ij}= \frac{1}{2}(\frac{\partial u_{i}}{\partial x_j} + \frac{\partial u_{j}}{\partial x_i})
# ```
# The displacement field is calculated as:
# ```math
# \boldsymbol u = \sum_{a=1}^{N_{\text{nodes}}} N_a \bar{\boldsymbol u}_{a} +
#  N_a ζ\frac{h}{2}(\theta_{a2} \boldsymbol e^{f}_{a1} - \theta_{a1} \boldsymbol e^{f}_{a2})
#
# ```
# The gradient of the displacement (in the lamina coordinate system), then becomes:
# ```math
# \frac{\partial u_{i}}{\partial x_j} = \sum_{m=1}^3 q_{im} \sum_{a=1}^{N_{\text{nodes}}} \frac{\partial N_a}{\partial x_j} \bar{u}_{am} +
#  \frac{\partial(N_a ζ)}{\partial x_j} \frac{h}{2} (\theta_{a2} e^{f}_{am1} - \theta_{a1} e^{f}_{am2})
# ```
function strain(dofvec::Vector{T}, N, dNdx, ζ, dζdx, q, ef1, ef2, h) where T

    u = reinterpret(Vec{3,T}, dofvec[1:12])
    θ = reinterpret(Vec{2,T}, dofvec[13:20])

    dudx = zeros(T, 3, 3)
    for m in 1:3, j in 1:3
        for a in 1:length(N)
            dudx[m,j] += dNdx[a][j] * u[a][m] + h/2 * (dNdx[a][j]*ζ + N[a]*dζdx[j]) * (θ[a][2]*ef1[a][m] - θ[a][1]*ef2[a][m])
        end
    end

    dudx = q*dudx
    ε = [dudx[1,1], dudx[2,2], dudx[1,2]+dudx[2,1], dudx[2,3]+dudx[3,2], dudx[1,3]+dudx[3,1]]
    return ε

end;

# ##### Main element routine
# Below is the main routine that calculates the stiffness matrix of the shell element.
# Since it is a so called degenerate shell element, the code is similar to that for an standard continuum element.
function integrate_shell!(ke, cv, qr_ooplane, X, data)
    nnodes = getnbasefunctions(cv)
    ndofs = nnodes*5
    h = data.thickness

    #Create the directors in each node.
    #Note: For a more general case, the directors should
    #be input parameters for the element routine.
    p = zeros(Vec{3}, nnodes)
    for i in 1:nnodes
        a = Vec{3}((0.0, 0.0, 1.0))
        p[i] = a/norm(a)
    end

    ef1, ef2, ef3 = fiber_coordsys(p)

    for iqp in 1:getnquadpoints(cv)

        dNdξ = cv.dNdξ[:,iqp]
        N = cv.N[:,iqp]

        for oqp in 1:length(qr_ooplane.weights)

            ζ = qr_ooplane.points[oqp][1]

            q = lamina_coordsys(dNdξ, ζ, X, p, h)
            J = getjacobian(q, N, dNdξ, ζ, X, p, h)
            Jinv = inv(J)

            dζdx = Vec{3}((0.0, 0.0, 1.0)) ⋅ Jinv
            dNdx = [Vec{3}((dNdξ[i][1], dNdξ[i][2], 0.0)) ⋅ Jinv for i in 1:nnodes]

            #For simplicity, use automatic differentiation to construct the B-matrix from the strain.
            B = ForwardDiff.jacobian(
                (a) -> strain(a, N, dNdx, ζ, dζdx, q, ef1, ef2, h), zeros(Float64, ndofs) )

            dV = det(J) * cv.qr_weights[iqp] * qr_ooplane.weights[oqp]
            ke .+= B'*data.C*B * dV
        end
    end
end;

# Run everything:
main()

ext: jl


text:
connector_comp_name(i::Int) = Symbol("ConnectorComp$i")

# Return the datatype to use for instance variables/parameters
function _instance_datatype(md::ModelDef, def::AbstractDatumDef)
    dtype = def.datatype == Number ? number_type(md) : def.datatype
    dims = dim_names(def)
    num_dims = dim_count(def)

    ti = get_time_index_position(def)

    if num_dims == 0
        T = ScalarModelParameter{dtype}

    elseif ti === nothing     # there's no time dimension
        T = Array{dtype, num_dims}

    else
        if isuniform(md)
            first, stepsize = first_and_step(md)
            first === nothing && @warn "_instance_datatype: first === nothing"
            T = TimestepArray{FixedTimestep{first, stepsize}, Union{dtype, Missing}, num_dims, ti}
        else
            times = time_labels(md)
            T = TimestepArray{VariableTimestep{(times...,)}, Union{dtype, Missing}, num_dims, ti}
        end
    end

    # @info "_instance_datatype returning $T"
    return T
end

# Create the Ref or Array that will hold the value(s) for a Parameter or Variable
function _instantiate_datum(md::ModelDef, def::AbstractDatumDef)
    dtype = _instance_datatype(md, def)
    dims = dim_names(def)
    num_dims = length(dims)

    # Scalar datum
    if num_dims == 0
        value = dtype(0)

    # Array datum, with :time dimension
    elseif dims[1] == :time

        if num_dims == 1
            value = dtype(dim_count(md, :time))
        else
            counts = dim_counts(md, Vector{Symbol}(dims))
            value = dtype <: AbstractArray ? dtype(undef, counts...) : dtype(counts...)
        end

    # Array datum, without :time dimension
    else
        # TBD: Handle unnamed indices properly
        counts = dim_counts(md, Vector{Symbol}(dims))
        value = dtype <: AbstractArray ? dtype(undef, counts...) : dtype(counts...)
    end

    return value
end

"""
    _instantiate_component_vars(md::ModelDef, comp_def::ComponentDef)

Instantiate a component `comp_def` in the model `md` and its variables (but not its
parameters). Return the resulting ComponentInstanceVariables.
"""
function _instantiate_component_vars(md::ModelDef, comp_def::ComponentDef)
    var_defs = variables(comp_def)

    names  = Symbol[nameof(def) for def in var_defs]
    values = Any[_instantiate_datum(md, def) for def in var_defs]
    types  = DataType[_instance_datatype(md, def) for def in var_defs]
    paths  = repeat(Any[comp_def.comp_path], length(names))

    return ComponentInstanceVariables(names, types, values, paths)
end

# Create ComponentInstanceVariables for a composite component from the list of exported vars
function _combine_exported_vars(comp_def::AbstractCompositeComponentDef, var_dict::Dict{ComponentPath, Any})
    names  = Symbol[]
    values = Any[]

    types = DataType[typeof(val) for val in values]
    paths = repeat(Any[comp_def.comp_path], length(names))
    ci_vars = ComponentInstanceVariables(names, types, values, paths)
    # @info "ci_vars: $ci_vars"]
    return ci_vars
end

function _combine_exported_pars(comp_def::AbstractCompositeComponentDef, par_dict::Dict{Tuple{ComponentPath, Symbol}, Any})
    names  = Symbol[]
    values = Any[]
    paths = repeat(Any[comp_def.comp_path], length(names))
    types = DataType[typeof(val) for val in values]
    return ComponentInstanceParameters(names, types, values, paths)
end

function _instantiate_vars(comp_def::ComponentDef, md::ModelDef, var_dict::Dict{ComponentPath, Any})
    var_dict[comp_def.comp_path] = _instantiate_component_vars(md, comp_def)
end

# Creates the top-level vars for the model
function _instantiate_vars(md::ModelDef, var_dict::Dict{ComponentPath, Any})
    _instantiate_vars(md, md, var_dict)
end

# Recursively instantiate all variables and store refs in the given dict.
function _instantiate_vars(comp_def::AbstractCompositeComponentDef, md::ModelDef, var_dict::Dict{ComponentPath, Any})
    comp_path = comp_def.comp_path
    # @info "_instantiate_vars composite $comp_path"

    for cd in compdefs(comp_def)
        _instantiate_vars(cd, md, var_dict)
    end
    var_dict[comp_path] = _combine_exported_vars(comp_def, var_dict)
end

# Do nothing if called on a leaf component
_collect_params(comp_def::ComponentDef, var_dict, par_dict) = nothing

# Recursively collect all parameters with connections to allocated storage for variables
function _collect_params(comp_def::AbstractCompositeComponentDef,
                         var_dict::Dict{ComponentPath, Any},
                         par_dict::Dict{Tuple{ComponentPath, Symbol}, Any})
    # depth-first search of composites
    for cd in compdefs(comp_def)
        _collect_params(cd, var_dict, par_dict)
    end

    # @info "Collecting params for $(comp_def.comp_id)"

    # Iterate over connections to create parameters, referencing storage in vars
    for ipc in internal_param_conns(comp_def)
        src_vars = var_dict[ipc.src_comp_path]
        var_value_obj = get_property_obj(src_vars, ipc.src_var_name)
        par_dict[(ipc.dst_comp_path, ipc.dst_par_name)] = var_value_obj
        # @info "internal conn: $(ipc.src_comp_path):$(ipc.src_var_name) => $(ipc.dst_comp_path):$(ipc.dst_par_name)"
    end

    for ext in external_param_conns(comp_def)
        param = external_param(comp_def, ext.external_param)
        par_dict[(ext.comp_path, ext.param_name)] = (param isa ScalarModelParameter ? param : value(param))
        # @info "external conn: $(ext.comp_name).$(ext.param_name) => $(param)"
    end

    # Make the external parameter connections for the hidden ConnectorComps.
    # Connect each :input2 to its associated backup value.
    for (i, backup) in enumerate(comp_def.backups)
        conn_comp = compdef(comp_def, connector_comp_name(i))
        conn_path = conn_comp.comp_path

        param = external_param(comp_def, backup)
        par_dict[(conn_path, :input2)] = (param isa ScalarModelParameter ? param : value(param))
    end
end

function _instantiate_params(comp_def::ComponentDef, par_dict::Dict{Tuple{ComponentPath, Symbol}, Any})
    # @info "Instantiating params for $(comp_def.comp_path)"
    comp_path = comp_def.comp_path
    names = parameter_names(comp_def)
    vals  = Any[par_dict[(comp_path, name)] for name in names]
    types = DataType[typeof(val) for val in vals]
    paths = repeat([comp_def.comp_path], length(names))

    return ComponentInstanceParameters(names, types, vals, paths)
end

function _instantiate_params(comp_def::AbstractCompositeComponentDef, par_dict::Dict{Tuple{ComponentPath, Symbol}, Any})
    _combine_exported_pars(comp_def, par_dict)
end

# Return a built leaf or composite LeafComponentInstance
function _build(comp_def::ComponentDef,
                var_dict::Dict{ComponentPath, Any},
                par_dict::Dict{Tuple{ComponentPath, Symbol}, Any},
                time_bounds::Tuple{Int, Int})
    # @info "_build leaf $(comp_def.comp_id)"
    # @info "  var_dict $(var_dict)"
    # @info "  par_dict $(par_dict)"

    pars = _instantiate_params(comp_def, par_dict)
    vars = var_dict[comp_def.comp_path]

    return LeafComponentInstance(comp_def, vars, pars, time_bounds)
end

function _build(comp_def::AbstractCompositeComponentDef,
                var_dict::Dict{ComponentPath, Any},
                par_dict::Dict{Tuple{ComponentPath, Symbol}, Any},
                time_bounds::Tuple{Int, Int})
    # @info "_build composite $(comp_def.comp_id)"
    # @info "  var_dict $(var_dict)"
    # @info "  par_dict $(par_dict)"

    comps = [_build(cd, var_dict, par_dict, time_bounds) for cd in compdefs(comp_def)]
    return CompositeComponentInstance(comps, comp_def, time_bounds)
end

function _build(md::ModelDef)
    # @info "_build(md)"
    add_connector_comps!(md)

    # check if all parameters are set
    not_set = unconnected_params(md)

    if ! isempty(not_set)
        params = join(not_set, "\n  ")
        error("Cannot build model; the following parameters are not set:\n  $params")
    end

    var_dict = Dict{ComponentPath, Any}()                 # collect all var defs and
    par_dict = Dict{Tuple{ComponentPath, Symbol}, Any}()  # store par values as we go

    _instantiate_vars(md, var_dict)
    _collect_params(md, var_dict, par_dict)

    # @info "var_dict: $var_dict"
    # @info "par_dict: $par_dict"

    t = dimension(md, :time)
    time_bounds = (firstindex(t), lastindex(t))

    propagate_time!(md, t)

    ci = _build(md, var_dict, par_dict, time_bounds)
    mi = ModelInstance(ci, md)
    return mi
end

function build(m::Model)
    # fix paths and propagate imports
    fix_comp_paths!(m.md)

    # Reference a copy in the ModelInstance to avoid changes underfoot
    md = deepcopy(m.md)
    m.mi = _build(md)
    m.md.dirty = false
    return nothing
end

"""
    create_marginal_model(base::Model, delta::Float64=1.0)

Create a `MarginalModel` where `base` is the baseline model and `delta` is the
difference used to create the `marginal` model.  Return the resulting `MarginaModel`
which shares the internal `ModelDef` between the `base` and `marginal`.
"""
function create_marginal_model(base::Model, delta::Float64=1.0)
    # Make sure the base has a ModelInstance before we copy since this
    # copies the ModelDef to avoid being affected by later changes.
    if ! is_built(base)
        build(base)
    end

    # Create a marginal model, which shares the internal ModelDef between base and marginal
    mm = MarginalModel(base, delta)
end

function Base.run(mm::MarginalModel; ntimesteps::Int=typemax(Int))
    run(mm.base, ntimesteps=ntimesteps)
    run(mm.marginal, ntimesteps=ntimesteps)
end

function build(mm::MarginalModel)
    build(mm.base)
    build(mm.marginal)
end


text:
OK = '✓'
FAIL= '✖'
RESET = "\033[0m"
PAD = " " ^ 2
FAINT_COLOR = "\033[90m"
PASS_COLOR = "\033[32m"
PASS_LIGHT_COLOR = "\033[92m"
FAILED_COLOR = "\033[31m"
FAILED_LIGHT_COLOR = "\033[91m"
HOUR = 3600_000
MINUTES = 60_000
SECOND = 1000
DESCRIPTION_ERROR_MESSAGE = "PANIC! Got an error while called describe"

toMilis(n) = int(n * 10e2)

function testColor()
  [println("\033[$(i)m Test_$i \033[0m") for i in 1:109]
end

function report(desc::Description)
  desc.depth == 1 && println()
  println(PAD ^ desc.depth, desc.name)
end

function summaryReport(passes::Array{Test,1},errors::Array{Union(Error,DescriptionError),1},elapsed)
  println()
  println(PASS_LIGHT_COLOR, PAD, length(passes), " passing ",FAINT_COLOR, "(", showTime(elapsed), ")")
  length(errors) == 0 || println(FAILED_LIGHT_COLOR, PAD, length(errors), " failing")
  println(RESET)
  length(errors) == 0  || fullFailedReport(errors)
end

function showTime(n)
  if n > HOUR
    string(n/HOUR)[1:3] * "h"
  elseif n > MINUTES
    string(n/MINUTES)[1:3] * "m"
  elseif n > SECOND
    string(n/SECOND)[1:3] * "s"
  else
    string(n) * "ms"
  end
end

function colorTime(elapsed::Int)
  if elapsed != 0
    color = elapsed > 500 ? 91 : (elapsed > 100 ? 93 : 90)
    " \e[$(color)m($(showTime(elapsed)))$RESET"
  else
    ""
  end
end

function passReport(test::Test, elapsed::Int)
  println(PAD ^ test.desc.depth,PAD, PASS_LIGHT_COLOR, OK, "\033[90m ", test.name, colorTime(elapsed), RESET)
end

function failedReport(err::DescriptionError)
  err.desc.depth == 1 && println()
  println(
    PAD ^ err.desc.depth,FAILED_COLOR, FAIL, " ",
    DESCRIPTION_ERROR_MESSAGE, " - ", err.desc.name,
    RESET
  )
end

function failedReport(err::Error)
  print(PAD ^ err.test.desc.depth,PAD,FAILED_COLOR, FAIL," ",err.test.name)
  print(RESET,'\n')
end

function fullFailedReport(errors::Array{Union(Error,DescriptionError),1})
  for i in 1:length(errors)
    err = errors[i]
    print(PAD, i, ") ")
    if isa(err, DescriptionError)
      print(DESCRIPTION_ERROR_MESSAGE, " - ", err.desc.name)
    else
      print(err.test.desc.name, " - ", err.test.name)
    end
    println(":", FAILED_COLOR)
    dump(err.err)
    println(RESET)
  end
end


text:
module ZulipReminderBot

using Base64
using HTTP
using JSON3
using StructTypes
using Dates
using DBInterface
using Strapping
using TimeZones
using Setfield

include("dotenv.jl")
include("miniorm.jl")
include("db_utils.jl")
include("migrations.jl")
include("zulipclient.jl")
include("parser.jl")
include("processors.jl")
include("server.jl")

export setupbot!
export DotEnv

precompile(cron_worker, (Channel{TimedMessage}, Channel{TimedMessage}))
precompile(HTTP.Handlers.serve, (Function, String, Int64))
precompile(JSON3.read, (Vector{UInt8}, ))
precompile(JSON3.read, (JSON3.VectorString{Vector{UInt8}}, ))
precompile(HTTP.request, (String, String, Vector{Pair{String, String}}, String))
precompile(HTTP.request, (String, String, Vector{Pair{String, String}}, SubString{String}))
precompile(query, (ZulipClient, String, Dict{Symbol, String}))

end # module


text:
module LibDeflate

using libdeflate_jll

# Error message strings.
const BAD_HEADER = "Bad header"
const UNTERMINATED_NULL_STRING = "Unterminated null string"
const HEADER_CRC16_CHECKSUM_DOES_NOT_MATCH = "Header CRC16 checksum does not match"
const PAYLOAD_CRC132_CHECKSUM_DOES_NOT_MATCH = "Payload CRC132 checksum does not match"
const OUTPUT_DATA_TOO_LONG = "Output data too long"
const INPUT_DATA_TOO_SHORT = "Input data too short"
const EXTRA_DATA_TOO_LONG = "Extra data too long"
const EXTRA_DATA_INVALID = "Extra data invalid"
const OUTPUT_BUFFER_TOO_SMALL = "Output buffer too small"

# Must be mutable for the GC to be able to interact with it
"""
    Decompressor()

Create an object which can decompress using the DEFLATE algorithm.
The same decompressor cannot be used by multiple threads at the same time.
To parallelize decompression, create multiple instances of `Decompressor`
and use one for each thread.

See also: [`decompress!`](@ref), [`unsafe_decompress!`](@ref)
"""
mutable struct Decompressor
    actual_nbytes_ret::UInt
    ptr::Ptr{Nothing}
end

Base.unsafe_convert(::Type{Ptr{Nothing}}, x::Decompressor) = x.ptr

function Decompressor()
    decompressor = Decompressor(0, ccall((:libdeflate_alloc_decompressor,
                   libdeflate), Ptr{Nothing}, ()))
    finalizer(free_decompressor, decompressor)
    return decompressor
end

function free_decompressor(decompressor::Decompressor)
    ccall((:libdeflate_free_decompressor, libdeflate),
           Nothing, (Ptr{Nothing},), decompressor)
    return nothing
end

"""
    Compressor(compresslevel::Int=6)

Create an object which can compress using the DEFLATE algorithm. `compresslevel`
can be from 1 (fast) to 12 (slow), and defaults to 6. The same compressor cannot
be used by multiple threads at the same time. To parallelize compression, create
multiple instances of `Compressor` and use one for each thread.

See also: [`compress!`](@ref), [`unsafe_compress!`](@ref)
"""
mutable struct Compressor
    level::Int
    ptr::Ptr{Nothing}
end

Base.unsafe_convert(::Type{Ptr{Nothing}}, x::Compressor) = x.ptr

function Compressor(compresslevel::Integer=6)
    compresslevel in 1:12 || throw(ArgumentError("Compresslevel must be in 1:12"))
    ptr = ccall((:libdeflate_alloc_compressor, libdeflate), Ptr{Nothing},
                (UInt,), compresslevel)
    compressor = Compressor(compresslevel, ptr)
    finalizer(free_compressor, compressor)
    return compressor
end

function free_compressor(compressor::Compressor)
    ccall((:libdeflate_free_compressor, libdeflate), Nothing, (Ptr{Nothing},), compressor)
    return nothing
end

# Compression and decompression functions
# Types and constants
const LIBDEFLATE_SUCCESS            = Cint(0)
const LIBDEFLATE_BAD_DATA           = Cint(1)
const LIBDEFLATE_SHORT_INPUT        = Cint(2)
const LIBDEFLATE_INSUFFICIENT_SPACE = Cint(3)

const BAD_DATA = "libdeflate 1: bad data"
const SHORT_INPUT = "libdeflate 2: short input"
const INSUFFICIENT_SPACE = "libdeflate 3: insufficient space"

"""
    LibDeflateError(message::String)

`LibDeflate` failed with `message`.
"""
struct LibDeflateError <: Exception
    msg::String
end

function translate_return_code(code)

    if code == LIBDEFLATE_BAD_DATA
        return BAD_DATA
    end

    if code == LIBDEFLATE_SHORT_INPUT
        return SHORT_INPUT
    end

    if code == LIBDEFLATE_INSUFFICIENT_SPACE
        return INSUFFICIENT_SPACE
    end

    throw(LibDeflateError("Unknown code"))
end

@noinline function check_return_code(code)
    iszero(code) && return nothing
    message = translate_return_code(code)
    throw(LibDeflateError(message))
end

# Raw C call - do not export this
function _unsafe_decompress!(decompressor::Decompressor,
                             outptr::Ptr{UInt8}, outlen::Integer,
                             inptr::Ptr{UInt8}, inlen::Integer, nptr::Ptr)
    status = ccall((:libdeflate_deflate_decompress, libdeflate), UInt,
                  (Ptr{Nothing}, Ptr{UInt8}, UInt, Ptr{UInt8}, UInt, Ptr{UInt}),
                   decompressor, inptr, inlen, outptr, outlen, nptr)
    check_return_code(status)
    return nothing
end

"""
    unsafe_decompress!(s::IteratorSize, ::Decompressor, outptr, n_out, inptr, n_in)

Decompress `n_in` bytes from `inptr` to `outptr` using the DEFLATE algorithm,
returning the number of decompressed bytes.
`s` gives whether you know the decompressed size or not.

If `s` isa `Base.HasLength`, the number of decompressed bytes is given as `n_out`.
This is more efficient, but will fail if the number is not correct.

If `s` isa `Base.SizeUnknown`, pass the number of available space in the output
to `n_out`.

See also: [`decompress!`](@ref)
"""
function unsafe_decompress! end

function unsafe_decompress!(::Base.HasLength, decompressor::Decompressor,
                            outptr::Ptr{UInt8}, n_out::Integer,
                            inptr::Ptr{UInt8}, n_in::Integer)
    _unsafe_decompress!(decompressor, outptr, n_out, inptr, n_in, C_NULL)
    return n_out
end

function unsafe_decompress!(::Base.SizeUnknown, decompressor::Decompressor,
                            outptr::Ptr{UInt8}, n_out::Integer,
                            inptr::Ptr{UInt8}, n_in::Integer)
    GC.@preserve decompressor begin
        retptr = pointer_from_objref(decompressor)
        _unsafe_decompress!(decompressor, outptr, n_out, inptr, n_in, retptr)
    end
    return decompressor.actual_nbytes_ret % Int
end

"""
    decompress!(::Decompressor, outdata, indata, [n_out::Integer]) -> Int

Use the passed `Decompressor` to decompress the byte vector `indata` into the
first bytes of `outdata` using the DEFLATE algorithm.
If the decompressed size is known beforehand, pass it as `n_out`. This will increase
performance, but will fail if it is wrong.

Return the number of bytes written to `outdata`.
"""
function decompress! end

# Decompress method with length known (preferred)
function decompress!(decompressor::Decompressor,
                     outdata::Vector{UInt8}, indata::Vector{UInt8}, n_out::Integer)
    if length(outdata) < n_out
        throw(ValueError("n_out must be less than or equal to length of outdata"))
    end
    GC.@preserve outdata indata unsafe_decompress!(Base.HasLength(),
        decompressor, pointer(outdata), n_out, pointer(indata), length(indata))
end

# Decompress method with length unknown (not preferred)
function decompress!(decompressor::Decompressor,
                     outdata::Vector{UInt8}, indata::Vector{UInt8})
    GC.@preserve outdata indata unsafe_decompress!(Base.SizeUnknown(),
        decompressor, pointer(outdata), length(outdata), pointer(indata), length(indata))
end

"""
    unsafe_compress(::Compressor, outptr, n_out, inptr, n_in)

Use the passed `Compressor` to compress `n_in` bytes from the pointer `inptr`
to the pointer `n_out`. If the compressed size is larger than the available
space `n_out`, throw an error.

See also: [`compress!`](@ref)
"""
function unsafe_compress!(compressor::Compressor, outptr::Ptr{UInt8}, n_out::Integer,
                          inptr::Ptr{UInt8}, n_in::Integer)
    bytes = ccall((:libdeflate_deflate_compress, libdeflate), UInt,
            (Ptr{Nothing}, Ptr{UInt8}, UInt, Ptr{UInt8}, UInt),
            compressor, inptr, n_in, outptr, n_out)

    if iszero(bytes)
        throw(LibDeflateError(OUTPUT_BUFFER_TOO_SMALL))
    end
    return bytes % Int
end

"""
    compress!(::Compressor, outdata, indata) -> Int

Use the passed `Compressor` to compress the byte vector `indata` into the first
bytes of `outdata` using the DEFLATE algorithm.

The output must fit in `outdata`. Return the number of bytes written to `outdata`.
"""
function compress!(compressor::Compressor,
                   outdata::Vector{UInt8}, indata::Vector{UInt8})
    GC.@preserve outdata indata unsafe_compress!(compressor, pointer(outdata),
        length(outdata), pointer(indata), length(indata))
end

"""
    unsafe_crc32(inptr, n_in, start) -> UInt32

Calculate the crc32 checksum of the first `n_in` of the pointer `inptr`,
with seed `start` (default is 0).
Note that crc32 is a different and slower algorithm than the `crc32c` provided
in the Julia standard library.

See also: [`crc32`](@ref)
"""
function unsafe_crc32(inptr::Ptr{UInt8}, n_in::Integer, start::UInt32=UInt32(0))
    return ccall((:libdeflate_crc32, libdeflate),
               UInt32, (UInt, Ptr{UInt8}, UInt), start, inptr, n_in)
end

"""
    crc32(data, start=UInt32(0)) -> UInt32

Calculate the crc32 checksum of the byte vector `data` and seed `start`.
Note that crc32 is a different and slower algorithm than the `crc32c` provided
in the Julia standard library.
"""
function crc32(data::Vector{UInt8}, start::UInt32=UInt32(0))
    GC.@preserve data unsafe_crc32(pointer(data), length(data), start)
end

include("gzip.jl")

export Decompressor,
       Compressor,

       LibDeflateError,
       unsafe_decompress!,
       decompress!,
       unsafe_gzip_decompress!,
       gzip_decompress!,
       unsafe_compress!,
       compress!,
       unsafe_gzip_compress!,
       gzip_compress!,
       unsafe_crc32,
       crc32,

       unsafe_parse_gzip_header,
       is_valid_extra_data

end # module


text:
#==============================================================================#
#
#               Routines for sampling trajectories from the
#               diffusion laws defined by this package.
#
#==============================================================================#
"""
    abstract type AbstractSDESolver end

Supertype of flags indicating ODE solvers
"""
abstract type AbstractSDESolver end

"""
    struct EulerMaruyama <: AbstractSDESolver end

Flag for indicating use of the Euler-Maruyama scheme for sampling diffusions.
IMPORTANT: this is the only diffusion path sampler implemented in this package.
There are no plans for implementing any other SDE solver in the forseeable
future.
"""
struct EulerMaruyama <: AbstractSDESolver end

"""
    struct Wiener{D,Tdevice}
    end

A struct defining the Wiener process. `D` indicates the dimension of the process
if it cannot be inferred from the DataType of the Trajectory. If the dimension
can be inferred then it takes precedence over the value of `D`.
"""
struct Wiener{D,T}
    Wiener{D}() where D = new{D,Float64}()
    Wiener{D,T}() where {D,T} = new{D,T}()
    Wiener(D::Integer=false, T::Type=Float64) = new{D,T}()
    Wiener(::DiffusionProcess{T,TP,TW}) where {T,TP,TW} = new{TW,T}()
end
wiener(args...) = Wiener(args...)

ismutable(el) = ismutable(typeof(el))
ismutable(::Type) = Val(false)
ismutable(::Type{K}) where K<:Array = Val(true)
ismutable(::Trajectory{T,Vector{K}}) where {T,K} = ismutable(K)

Base.zero(w::Wiener{D,T}) where {D,T} = zero(SVector{D,T})

#===============================================================================
                        Sampling Wiener processes
===============================================================================#
"""
    Base.rand(w::Wiener, tt, y1)

Samples Wiener process on `tt`, started from `y1` and returns a new object with
a sampled trajectory.
"""
Base.rand(w::Wiener, tt, y1=zero(w)) = rand(Random.GLOBAL_RNG, w, tt, y1)

function Base.rand(
        rng::Random.AbstractRNG,
        w::Wiener{D},
        tt,
        y1::K=zero(w),
    ) where {D,K}
    path = trajectory(tt, K, D, ismutable(K))
    rand!(rng, w, path, y1)
end

"""
    Random.rand!(
        rng::Random.AbstractRNG,
        ::Wiener,
        path::Trajectory{T,Vector{K}},
        y1=zero(K)
    ) where {T,K}

Samples Wiener process with immutable states, started from `y1` and saves the
data in `path`. `rng` is used as a random number generator.
"""
function Random.rand!(
        rng::Random.AbstractRNG,
        ::Wiener,
        path::Trajectory{T,<:Vector{K}},
        y1=zero(K)
    ) where {T,K}
    tt, yy = path.t, path.x

    N = length(path)
    yy[1] = y1
    for i in 2:N
        rootdt = sqrt(tt[i] - tt[i-1])
        yy[i] = yy[i-1] + rootdt*randn(rng, K)
    end
    path
end

"""
    Random.rand!(
        rng::Random.AbstractRNG,
        ::Wiener{D},
        path::Trajectory{T,Vector{Vector{K}}},
        y1=zeros(K,D)
    ) where {K,T,D}

Samples Wiener process with mutable states, started from `y1` and saves the data
in `path`. `rng` is used as a random number generator.
"""
function Random.rand!(
        rng::Random.AbstractRNG,
        ::Wiener{D},
        path::Trajectory{T,<:Vector{<:Array{K}}},
        y1=zeros(K,D)
    ) where {K,T,D}
    yy, tt = path.x, path.t
    N = length(path)
    yy[1] = y1
    for i in 2:N
        rootdt = sqrt(tt[i] - tt[i-1])
        for j in eachindex(yy[i])
            yy[i][j] = yy[i-1][j] + rootdt*randn(rng, K)
        end
    end
    path
end

"""
    Random.rand!(
        w::Wiener{D},
        path::Trajectory{T,Vector{K}},
        y1=zero(K,D,ismutable(K))
    ) where {T,K,D}

Samples Wiener process started from `y1` and saves the data in `path`. Uses
a default random number generator.
"""
function Random.rand!(
        w::Wiener{D},
        path::Trajectory{T,Vector{K}},
        y1=zero(K, D, ismutable(K))
    ) where {T,K,D}
    rand!(Random.GLOBAL_RNG, w, path, y1)
end

#===============================================================================
                    Euler-Maruyama scheme for solving SDEs
===============================================================================#

value(x) = x
value!(y, x) = (y.= x)
value(x::SVector{N,ForwardDiff.Dual{K,T,M}}) where {N,K,T,M} = SVector{N,T}(
    map(x->x.value, x)
)
value!(y, x::AbstractArray{<:ForwardDiff.Dual}) = map!(x->x.value, y, x)

struct _DEFAULT_F end
const __DEFAULT_F = _DEFAULT_F()
Base.getindex(f::_DEFAULT_F, i::Int) = f
Base.setindex!(f::_DEFAULT_F, v, i::Int) = nothing
@inline (::_DEFAULT_F)(args...) = nothing
const _FINAL = Val(:final)

"""
    solve!(XX, WW, P, y1)

Compute a trajectory, started from `y1` and following the diffusion law `P`,
from the sampled Wiener process `WW`. Save the sampled path in `XX`. Return
prematurely with a `false` massage if the numerical scheme has led to the solver
violating the state-space restrictions.
"""
function solve!(XX, WW, P, y1; f=__DEFAULT_F)
    solve!(EulerMaruyama(), ismutable(XX), XX, WW, P, y1; f=f)
end

"""
    solve!(XX, WW, P, y1, buffer)

Same as `solve!(XX, WW, P, y1)`, but additionally provides a pre-allocated
buffer for performing in-place computations.
"""
function solve!(XX, WW, P, y1, buffer; f=__DEFAULT_F)
    solve!(EulerMaruyama(), ismutable(XX), XX, WW, P, y1, buffer; f=f)
end

function solve!(
        ::EulerMaruyama, ::Val{false}, XX, WW, P, y1::Ky1;
        f=__DEFAULT_F
    ) where Ky1
    yy, ww, tt = XX.x, WW.x, XX.t
    N = length(XX)

    yy[1] = value(y1)
    y = y1 # more appropriate name
    f_accum = f(P, y)

    for i in 2:N
        dt = tt[i] - tt[i-1]
        dW = ww[i] - ww[i-1]
        f_accum = f(f_accum, P, y, tt[i-1], dt, dW, i-1)

        y = y + _b((tt[i-1], i-1), y, P)*dt + _σ((tt[i-1], i-1), y, P)*dW

        yy[i] = value(y) # strip duals
        bound_satisfied(P, yy[i]) || return false, nothing
    end
    f_accum = f(f_accum, P, y, _FINAL)
    true, f_accum
end

function solve!(
        ::EulerMaruyama, ::Val{true}, XX, WW, P, y1::K,
        buffer=StandardEulerBuffer{K}(P);
        f=__DEFAULT_F
    ) where K
    yy, ww, tt = XX.x, WW.x, XX.t
    N = length(XX)

    value!(yy[1], y1)
    y::K = y1 # name alias
    f_accum = f(buffer, P, y)

    for i in 2:N
        dt = tt[i] - tt[i-1]
        buffer.dW .= ww[i] .- ww[i-1]
        f(buffer, P, y, tt[i-1], dt, i-1)

        _b!(buffer, (tt[i-1], i-1), y, P)
        _σ!(buffer, (tt[i-1], i-1), y, P)

        mul!(buffer.y, buffer.σ, buffer.dW)
        y .= y .+ buffer.y .+ buffer.b .* dt

        value!(yy[i], y) # strip duals
        bound_satisfied(P, yy[i]) || return false, nothing
    end
    f(buffer, P, y, _FINAL)
    true, f_accum
end

#=
function solve!(
        ::EulerMaruyama, ::Val{true}, XX, WW, P, y1::K,
        buffer=StandardEulerBuffer{K}(P)
    ) where K
    yy, ww, tt = XX.x, WW.x, XX.t
    N = length(XX)

    value!(yy[1], y1)
    y = y1 # name alias
    y_temp = get_y(buffer)
    b_temp = get_b(buffer)
    σ_temp = get_σ(buffer)
    dW = get_dW(buffer)

    for i in 2:N
        dt = tt[i] - tt[i-1]

        _b!(buffer, (tt[i-1], i-1), y, P)
        _σ!(buffer, (tt[i-1], i-1), y, P)

        dW .= ww[i] .- ww[i-1]
        mul!(y_temp, σ_temp, dW)
        y .+= y_temp .+ b_temp .* dt

        value!(yy[i], y) # strip duals
        bound_satisfied(P, yy[i]) || return false, nothing
    end
    true, y
end
=#

#===============================================================================
            Computing Wiener path reproducing given trajectory
===============================================================================#

function invsolve!(XX, WW, P)
    invsolve!(EulerMaruyama(), ismutable(XX), XX, WW, P)
end

function invsolve!(XX, WW, P, buffer)
    invsolve!(EulerMaruyama(), ismutable(XX), XX, WW, P, buffer)
end

function invsolve!(::EulerMaruyama, ::Val{false}, XX, WW, P)
    yy, ww, tt = XX.x, WW.x, XX.t
    N = length(XX)

    ww[N] = zero(ww[N])

    for i in N-1:-1:1
        yᵢ₊₁ = yy[i+1]
        yᵢ = yy[i]
        dt = tt[i+1] - tt[i]

        ww[i] = ww[i+1] - (
            nonhypo_σ((tt[i], i), yᵢ, P)\
            nonhypo(yᵢ₊₁ - yᵢ - _b((tt[i], i), yᵢ, P)*dt, P)
        )
    end
    for i in N:-1:1
        ww[i] = ww[i] - ww[1]
    end
    WW
end

#===============================================================================
        Sampling Diffusion processes using the Euler Maruyama scheme
===============================================================================#

function Base.rand(P::DiffusionProcess, tt, y1=zero(P); f=__DEFAULT_F)
    rand(Random.GLOBAL_RNG, P, tt, y1, ismutable(y1); f=f)
end

function Base.rand(
        rng::Random.AbstractRNG,
        P::DiffusionProcess{T,DP,DW},
        tt,
        y1::K,
        v::Val{true};
        f=__DEFAULT_F
    ) where {T,DP,DW,K}
    Wnr = Wiener{DW}()
    w0 = zeros(eltype(y1), DW)
    WW = trajectory(tt, typeof(w0), DW, v) #rand(rng, Wnr, tt, w0) # it samples wiener twice
    XX = trajectory(tt, K, DP, v)
    success, f_accum = false, nothing
    buffer = StandardEulerBuffer{K}(P)
    while !success
        rand!(rng, Wnr, WW, w0)
        success, f_accum = solve!(XX, WW, P, y1, buffer; f=f)
    end
    typeof(f) != _DEFAULT_F && return XX, f_accum
    XX
end

StaticArrays.similar_type(::Type{K}, s::Size) where K <: Number = SVector{s[1],K}

function Base.rand(
        rng::Random.AbstractRNG,
        P::DiffusionProcess{T,DP,DW},
        tt,
        y1::K,
        v::Val{false};
        f=__DEFAULT_F
    ) where {T,DP,DW,K}
    w0 = (
        default_wiener_type(P) <: Number ?
        zero(eltype(K)) :
        zero(similar_type(K, Size(DW)))
    )
    Wnr = Wiener()
    WW = trajectory(tt, typeof(w0), DW, v) #rand(rng, Wnr, tt, w0)
    XX = trajectory(tt, K, DP, v)
    success, f_accum = false, nothing
    while !success
        rand!(rng, Wnr, WW, w0)
        success, f_accum = solve!(XX, WW, P, y1; f=f)
    end
    typeof(f) != _DEFAULT_F && return XX, f_accum
    XX
end

#===============================================================================
                    Computing gradients of functionals
===============================================================================#

"""
    grad_θ(θ, y1, W, X, Law, f)

Compute ∇`f` with respect to parameters `θ` for a fixed Wiener path `W`. `X` is
a container where the the trajectory computed for the Wiener path `W` under the
law `Law`(`θ`) will be stored. `y1` is the starting position.
"""
function grad_θ(θ, y1, W, X, Law, f)
    function foo(θ°)
        P = Law(θ°...)
        _, foo_result = solve!(X, W, P, y1_θ; f=f)
        foo_result
    end
    K = eltype(y1)
    y1_θ_type = similar_type(θ, Dual{Tag{typeof(foo),K},K,length(θ)}, Size(y1))
    y1_θ = y1_θ_type(y1)

    ForwardDiff.gradient(foo, θ)
end

"""
    grad_y1(y1, W, X, P, f)

Compute ∇`f` with respect to the starting position `y1` for a fixed Wiener path
`W`. `X` is a container where the the trajectory computed for the Wiener path
`W` under the law `P` will be stored.
"""
function grad_y1(y1, W, X, P, f)
    function foo(y1°)
        _, foo_result = solve!(X, W, P, y1°; f=f)
        foo_result
    end
    ForwardDiff.gradient(foo, y1)
end


text:
__precompile__()

module OrbitalInsertion

using DifferentialEquations, ForwardDiff, DiffResults, Printf, LinearAlgebra, NLopt, Roots, Interpolations, Plots

abstract type System end

struct EMsystem <:System end
struct SEsystem <:System end

struct ProblemInfo{T<:Float64}
    mu  ::T
    rp  ::T
    rs  ::T
    r12 ::T
    w   ::T
    mup ::T
    x1  ::T
    x2  ::T
end
struct ProblemParameters{T<:Float64,J<:Int64,K<:Bool}
    tspan   ::T
    hgoal   ::T
    hmmin   ::T
    hmmax   ::T
    sType   ::K
    points  ::J
    stop    ::K
    maxt0   ::T
    rest    ::Tuple{J,J}
    lim     ::T
    cTol    ::T
    tRest   ::K
    odeAlg  ::K
    atol    ::T
    rtol    ::T
    retorno ::K
    detMin  ::T
    x0      ::Array{T,1}
    y       ::Array{T,1}
    optTol  ::T
    neval   ::J
    col     ::K
end
mutable struct SystemInfo{T<:Float64}
    mu  ::T
    Rp  ::T
    Rs  ::T
    r12 ::T
    m1  ::T
    G   ::T
end

export createConstants, createSystem, createPar, search, filterGlobalSol, point, plotSols
export EMsystem, SEsystem
export ProblemInfo, ProblemParameters, SystemInfo, ProblemParameters

include("funcs_insert.jl")
include("systems.jl")
include("restrictions.jl")

end # module


text:
__precompile__()

module Currencies

using CurrenciesBase

# Exports
export AbstractMonetary, Monetary
export currency, decimals, majorunit, @usingcurrencies
export currencyinfo, iso4217num, iso4217alpha, shortsymbol, longsymbol
export valuate, ExchangeRateTable, ecbrates
export Basket, StaticBasket, DynamicBasket
export simplefv, compoundfv
export newcurrency!, @usingcustomcurrency
export format

# DeclarativeFormatting (not specific to Currencies; under development)
include("DeclarativeFormatting/DeclarativeFormatting.jl")

include("Baskets/Basket.jl")
using .Baskets

include("Valuation/Valuation.jl")
using .Valuation

# Interface (display/formatting, convenience macro)
include("CurrencyFormatting/CurrencyFormatting.jl")
using .CurrencyFormatting

# Deprecations
include("deprecated.jl")

end # module


text:
using Napari
using PyCall
using TestImages
using Test

try
    global skimagedata = pyimport("skimage.data")
    global astronaut = skimagedata.astronaut()
    global moon = skimagedata.moon()
    global camera = skimagedata.camera()
    global MagickMock = pyimport("unittest.mock").MagicMock
catch err
    @info "Napari.jl tests use data from skimage.data, but they could not be loaded"
    rethrow(err)
end

@testset "Napari.jl" begin
    @testset "PyCall and Python environment" begin
       @test py"""
       from skimage import data
       import napari

       astronaut = data.astronaut()
       viewer = napari.view_image( astronaut )
       viewer.close()
       """ === nothing
    end

    @testset "PyObject Interface" begin
        viewer = napari.Viewer()
        @test viewer isa PyObject
        @test viewer.add_image( astronaut ) isa PyObject
        @test viewer.add_image( moon ) isa PyObject
        @test viewer.add_image( camera ) isa PyObject
        viewer.close()
        viewer = napari.view_image( astronaut , name = "Eileen Collins" )
        @test viewer isa PyObject
        viewer.close()
    end

    @testset "@namekw" begin
        viewer = @namekw napari.view_image( astronaut )
        @test viewer isa PyObject
        @test ( @namekw viewer.add_image( moon ) ) isa PyObject
        @test ( @namekw viewer.add_image( camera ) ) isa PyObject
        viewer.close()
    end

    @testset "Add and View Macros" begin
        viewer = @view_image astronaut
        @test viewer isa PyObject
        viewer.close()
        viewer = @view_image( astronaut, name ="Eileen Collins" )
        @test viewer isa PyObject
        @test @add_image(viewer, moon) isa PyObject
        @test @add_image(viewer, camera) isa PyObject
    end

    @testset "TestImages" begin
        viewer = MagickMock(napari.Viewer)
        for f in TestImages.remotefiles[1:10]
            @info f
            @test viewer.add_image( testimage(f), name=f) isa PyObject
        end
        for f in TestImages.remotefiles[11:20]
            @info f
            @test viewer.add_image( testimage(f), name=f) isa PyObject
        end
        for f in TestImages.remotefiles[21:30]
            @info f
            @test viewer.add_image( testimage(f), name=f) isa PyObject
        end
        for f in TestImages.remotefiles[31:40]
            @info f
            @test viewer.add_image( testimage(f), name=f) isa PyObject
        end
        for f in TestImages.remotefiles[41:end]
            @info f
            @test viewer.add_image( testimage(f), name=f) isa PyObject
        end
        @info "Shepp Logan"
        @test viewer.add_image( TestImages.shepp_logan(512, 512) ) isa PyObject
    end
end

function view_all_test_images()
    viewer = napari.Viewer()
    for (i,f) in enumerate(TestImages.remotefiles)
        @add_image(viewer, testimage(f), name="$i: $f")
    end
end


text:
module Chemostat_Nanchen2006

    import BSON
    import DrWatson
    import Chemostat
    const Ch = Chemostat
    const ChU = Ch.Utils
    const ChSS = Ch.SteadyState
    const ChLP = Ch.LP

    using ProjAssistant
    @gen_top_proj

    include("NanchenData/NanchenData.jl")
    include("Utils/Utils.jl")
    include("BegData/BegData.jl")
    include("iJR904/iJR904.jl")

end


text:
module LifConfig

const Float = Float64
const Time = Float

export Float, Time

end


text:
module TestQrValue
   using Base.Test
   using Brim
   #=testdir = dirname(@__FILE__)=#
   A = [1 0 0; 0 1 0; 0 0 1]
   S = [1 0 0; 0 1 0; 0 0 1; 1 0 0; 0 1 0; 0 0 1]
   M = Modular(A, S)
   # Approximate values, calculated by hand
   @test M |> Qr == 1.0
end


text:
using EffectSizes
using Test
using EffectSizes: bootstrapsample, twotailedquantile, AbstractConfidenceInterval,
    ConfidenceInterval, BootstrapConfidenceInterval

@testset "AbstractConfidenceInterval" begin
    l = .1
    u = .9
    q = .95
    ci = ConfidenceInterval(l, u, q)
    @test confint(ci) == (l, u)
    @test lower(ci) == l
    @test upper(ci) == u
    @test quantile(ci) == q
    @test_throws ArgumentError ConfidenceInterval(1, -1, .95)
    @test_throws ArgumentError ConfidenceInterval(1., 0.999, .95)

    @testset "ConfidenceInterval" begin
        l = .1
        u = .9
        q = .95
        ci = ConfidenceInterval(l, u, q)
        @test confint(ci) == (l, u)
        @test lower(ci) == l
        @test upper(ci) == u
        @test quantile(ci) == q
        @test_throws ArgumentError ConfidenceInterval(1, -1, .95)
        @test_throws ArgumentError ConfidenceInterval(1., 0.999, .95)
        @test_throws DomainError ConfidenceInterval(l, u, -0.1)
        @test_throws DomainError ConfidenceInterval(l, u, 1.1)

        xs = randn(90)
        ys = randn(110)
        es = 0.2
        q = 0.9
        ci = ConfidenceInterval(xs, ys, es, quantile=q)
        @test ci isa AbstractConfidenceInterval
        @test quantile(ci) == q
        @test lower(ci) < upper(ci)
        @test ci == ConfidenceInterval(xs, ys, es, quantile=q)
        ci2 = ConfidenceInterval(xs, ys, es, quantile=0.8)
        @test ci !== ci2
        @test lower(ci2) > lower(ci)
        @test upper(ci2) < upper(ci)
        @test_throws DomainError ConfidenceInterval(xs, ys, es, quantile=1.1)
    end

    @testset "BootstrapConfidenceInterval" begin
        l = .1
        u = .9
        q = .95
        b = 10^4
        ci = BootstrapConfidenceInterval(l, u, q, b)
        @test confint(ci) == (l, u)
        @test lower(ci) == l
        @test upper(ci) == u
        @test quantile(ci) == q
        @test ci.bootstrap == b
        @test_throws ArgumentError BootstrapConfidenceInterval(1, -1, q, b)
        @test_throws ArgumentError BootstrapConfidenceInterval(1., 0.999, q, b)
        @test_throws DomainError BootstrapConfidenceInterval(l, u, -0.1, b)
        @test_throws DomainError BootstrapConfidenceInterval(l, u, 1.1, b)
        @test_throws DomainError BootstrapConfidenceInterval(l, u, q, -1)

        xs = randn(90)
        ys = randn(110)
        q = 0.8
        f(xs, ys) = sum([sum(xs), sum(ys)] ./ 1000)
        ci = BootstrapConfidenceInterval(f, xs, ys, 100, quantile=q)
        @test ci isa AbstractConfidenceInterval
        @test quantile(ci) == q
        @test lower(ci) < upper(ci)
        @test_throws DomainError BootstrapConfidenceInterval(f, xs, ys, 100, quantile=1.1)
        @test_throws DomainError BootstrapConfidenceInterval(f, xs, ys, -1, quantile=q)
    end
end

@testset "bootstrapsample" begin
    for _ = 1:100
        xs = rand(100)
        ys = bootstrapsample(xs)
        @test length(ys) == length(xs)
        @test Set(ys) ≤ Set(xs)
    end
end

@testset "twotailedquantile" begin
    l, u = twotailedquantile(.95)
    @test l ≈ 0.025
    @test u ≈ 0.975
    @test l + u == 1

    l, u = twotailedquantile(.9)
    @test l ≈ 0.05
    @test u ≈ 0.95
    @test l + u == 1

    for q = rand(100)
        l, u = twotailedquantile(q)
        @test l + u == 1
    end

    @test_throws DomainError twotailedquantile(-.1)
    @test_throws DomainError twotailedquantile(1.1)
end


text:
using ModelingToolkit
using Base.Test

d = DiscreteDomain(5) * Interval(0.0,5.0)
typeof(d) <: ProductDomain
typeof(d.domains[1]) <: DiscreteDomain
typeof(d.domains[2]) <: Interval


text:
function cleanurl(url::Union{SubString{String},String})::String

  indexinit::Int64 = 1
  indexend::Int64 = length(url)

  if url[1] == '/'
  	indexinit=2
  end

  if url[indexend] == '/'
    indexend = indexend-1
  end

  return url[indexinit:indexend]

end