text:
"""
    DArray(init, dims, [procs, dist])

Construct a distributed array.

The parameter `init` is a function that accepts a tuple of index ranges.
This function should allocate a local chunk of the distributed array and initialize it for the specified indices.

`dims` is the overall size of the distributed array.

`procs` optionally specifies a vector of process IDs to use.
If unspecified, the array is distributed over all worker processes only. Typically, when running in distributed mode,
i.e., nprocs() > 1, this would mean that no chunk of the distributed array exists on the process hosting the
interactive julia prompt.

`dist` is an integer vector specifying how many chunks the distributed array should be divided into in each dimension.

For example, the `dfill` function that creates a distributed array and fills it with a value `v` is implemented as:

### Example
```jl
dfill(v, args...) = DArray(I->fill(v, map(length,I)), args...)
```
"""
mutable struct DArray{T,N,A} <: AbstractArray{T,N}
    id::Tuple
    dims::NTuple{N,Int}
    pids::Array{Int,N}                          # pids[i]==p ⇒ processor p has piece i
    indices::Array{NTuple{N,UnitRange{Int}},N}  # indices held by piece i
    cuts::Vector{Vector{Int}}                   # cuts[d][i] = first index of chunk i in dimension d
    localpart::Union{A,Nothing}
    release::Bool

    function DArray{T,N,A}(id, dims, pids, indices, cuts, lp) where {T,N,A}
        # check invariants
        if dims != map(last, last(indices))
            throw(ArgumentError("dimension of DArray (dim) and indices do not match"))
        end
        release = (myid() == id[1])

        d = d_from_weakref_or_d(id)
        if d === nothing
            d = new(id, dims, pids, indices, cuts, lp, release)
        end

        if release
            push!(refs, id)
            registry[id] = WeakRef(d)

#            println("Installing finalizer for : ", d.id, ", : ", object_id(d), ", isbits: ", isbits(d))
            finalizer(close, d)
        end
        d
    end

    DArray{T,N,A}() where {T,N,A} = new()
end

function d_from_weakref_or_d(id)
    d = get(registry, id, nothing)
    isa(d, WeakRef) && return d.value
    return d
end

Base.eltype(::Type{DArray{T}}) where {T} = T
empty_localpart(T,N,A) = A(Array{T}(undef, ntuple(zero, N)))

const SubDArray{T,N,D<:DArray} = SubArray{T,N,D}
const SubOrDArray{T,N} = Union{DArray{T,N}, SubDArray{T,N}}

localtype(::Type{DArray{T,N,S}}) where {T,N,S} = S
localtype(::Type{SubDArray{T,N,D}}) where {T,N,D} = localtype(D)
localtype(A::SubOrDArray) = localtype(typeof(A))
localtype(A::AbstractArray) = typeof(A)

Base.hash(d::DArray, h::UInt) = Base.hash(d.id, h)

## core constructors ##

function DArray(id, init, dims, pids, idxs, cuts)
    localtypes = Vector{DataType}(undef,length(pids))

    @sync begin
        for i = 1:length(pids)
            @async begin
                local typA
                if isa(init, Function)
                    typA = remotecall_fetch(construct_localparts, pids[i], init, id, dims, pids, idxs, cuts)
                else
                    # constructing from an array of remote refs.
                    typA = remotecall_fetch(construct_localparts, pids[i], init[i], id, dims, pids, idxs, cuts)
                end
                localtypes[i] = typA
            end
        end
    end

    if length(unique(localtypes)) != 1
        @sync for p in pids
            @async remotecall_fetch(release_localpart, p, id)
        end
        throw(ErrorException("Constructed localparts have different `eltype`: $(localtypes)"))
    end
    A = first(localtypes)

    if myid() in pids
        d = registry[id]
        d = isa(d, WeakRef) ? d.value : d
    else
        T = eltype(A)
        N = length(dims)
        d = DArray{T,N,A}(id, dims, pids, idxs, cuts, empty_localpart(T,N,A))
    end
    d
end

function construct_localparts(init, id, dims, pids, idxs, cuts; T=nothing, A=nothing)
    localpart = isa(init, Function) ? init(idxs[localpartindex(pids)]) : fetch(init)
    if A == nothing
        A = typeof(localpart)
    end
    if T == nothing
        T = eltype(A)
    end
    N = length(dims)
    d = DArray{T,N,A}(id, dims, pids, idxs, cuts, localpart)
    registry[id] = d
    A
end

function ddata(;T::Type=Any, init::Function=I->nothing, pids=workers(), data::Vector=[])
    pids=sort(vec(pids))
    id = next_did()
    npids = length(pids)
    ldata = length(data)
    idxs, cuts = chunk_idxs([npids], [npids])

    if ldata > 0
        @assert rem(ldata,npids) == 0
        if ldata == npids
            T = eltype(data)
            s = DestinationSerializer(pididx->data[pididx], pids)
            init = I->localpart(s)
        else
            # call the standard distribute function
            return distribute(data)
        end
    end

    @sync for i = 1:length(pids)
        @async remotecall_fetch(construct_localparts, pids[i], init, id, (npids,), pids, idxs, cuts; T=T, A=T)
    end

    if myid() in pids
        d = registry[id]
        d = isa(d, WeakRef) ? d.value : d
    else
        d = DArray{T,1,T}(id, (npids,), pids, idxs, cuts, nothing)
    end
    d
end

function gather(d::DArray{T,1,T}) where T
    a=Array{T}(undef, length(procs(d)))
    @sync for (i,p) in enumerate(procs(d))
        @async a[i] = remotecall_fetch(localpart, p, d)
    end
    a
end

function DArray(init, dims, procs, dist)
    np = prod(dist)
    procs = reshape(procs[1:np], ntuple(i->dist[i], length(dist)))
    idxs, cuts = chunk_idxs([dims...], dist)
    id = next_did()

    return DArray(id, init, dims, procs, idxs, cuts)
end

function DArray(init, dims, procs)
    if isempty(procs)
        throw(ArgumentError("no processors given"))
    end
    return DArray(init, dims, procs, defaultdist(dims, procs))
end
DArray(init, dims) = DArray(init, dims, workers()[1:min(nworkers(), maximum(dims))])

# Create a DArray from a collection of references
# The refs must have the same layout as the parts distributed.
# i.e.
#    size(refs) must specify the distribution of dimensions across processors
#    prod(size(refs)) must equal number of parts
# FIXME : Empty parts are currently not supported.
function DArray(refs)
    dimdist = size(refs)
    id = next_did()

    npids = [r.where for r in refs]
    nsizes = Array{Tuple}(undef, dimdist)
    @sync for i in 1:length(refs)
        let i=i
            @async nsizes[i] = remotecall_fetch(sz_localpart_ref, npids[i], refs[i], id)
        end
    end

    nindices = Array{NTuple{length(dimdist),UnitRange{Int}}}(undef, dimdist...)

    for i in 1:length(nindices)
        subidx = CartesianIndices(dimdist)[i]
        nindices[i] = ntuple(length(subidx)) do x
            idx_in_dim = subidx[x]
            startidx = 1
            for j in 1:(idx_in_dim-1)
                prevsubidx = ntuple(y -> y == x ? j : subidx[y], length(subidx))
                prevsize = nsizes[prevsubidx...]
                startidx += prevsize[x]
            end
            startidx:startidx+(nsizes[i][x])-1
        end
    end

    lastidxs = hcat([Int[last(idx_in_d)+1 for idx_in_d in idx] for idx in nindices]...)
    ncuts = Array{Int,1}[pushfirst!(sort(unique(lastidxs[x,:])), 1) for x in 1:length(dimdist)]
    ndims = tuple([sort(unique(lastidxs[x,:]))[end]-1 for x in 1:length(dimdist)]...)

    DArray(id, refs, ndims, reshape(npids, dimdist), nindices, ncuts)
end

macro DArray(ex0::Expr)
    if ex0.head !== :comprehension
        throw(ArgumentError("invalid @DArray syntax"))
    end
    ex = ex0.args[1]
    if ex.head !== :generator
        throw(ArgumentError("invalid @DArray syntax"))
    end
    ex.args[1] = esc(ex.args[1])
    ndim = length(ex.args) - 1
    ranges = map(r->esc(r.args[2]), ex.args[2:end])
    for d = 1:ndim
        var = ex.args[d+1].args[1]
        ex.args[d+1] = :( $(esc(var)) = ($(ranges[d]))[I[$d]] )
    end
    return :( DArray((I::Tuple{Vararg{UnitRange{Int}}})->($ex0),
                tuple($(map(r->:(length($r)), ranges)...))) )
end

# new DArray similar to an existing one
DArray(init, d::DArray) = DArray(next_did(), init, size(d), procs(d), d.indices, d.cuts)

sz_localpart_ref(ref, id) = size(fetch(ref))

Base.similar(d::DArray, T::Type, dims::Dims) = DArray(I->Array{T}(undef, map(length,I)), dims, procs(d))
Base.similar(d::DArray, T::Type) = similar(d, T, size(d))
Base.similar(d::DArray{T}, dims::Dims) where {T} = similar(d, T, dims)
Base.similar(d::DArray{T}) where {T} = similar(d, T, size(d))

Base.size(d::DArray) = d.dims

chunktype(d::DArray{T,N,A}) where {T,N,A} = A

## chunk index utilities ##

# decide how to divide each dimension
# returns size of chunks array
function defaultdist(dims, pids)
    dims = [dims...]
    chunks = ones(Int, length(dims))
    np = length(pids)
    f = sort!(collect(keys(factor(np))), rev=true)
    k = 1
    while np > 1
        # repeatedly allocate largest factor to largest dim
        if np % f[k] != 0
            k += 1
            if k > length(f)
                break
            end
        end
        fac = f[k]
        (d, dno) = findmax(dims)
        # resolve ties to highest dim
        dno = findlast(isequal(d), dims)
        if dims[dno] >= fac
            dims[dno] = div(dims[dno], fac)
            chunks[dno] *= fac
        end
        np = div(np, fac)
    end
    return chunks
end

# get array of start indices for dividing sz into nc chunks
function defaultdist(sz::Int, nc::Int)
    if sz >= nc
        chunk_size = div(sz,nc)
        remainder = rem(sz,nc)
        grid = zeros(Int64, nc+1)
        for i = 1:(nc+1)
            grid[i] += (i-1)*chunk_size + 1
            if i<= remainder
                grid[i] += i-1
            else
                grid[i] += remainder
            end
        end
        return grid
    else
        return [[1:(sz+1);]; zeros(Int, nc-sz)]
    end
end

# compute indices array for dividing dims into chunks
function chunk_idxs(dims, chunks)
    cuts = map(defaultdist, dims, chunks)
    n = length(dims)
    idxs = Array{NTuple{n,UnitRange{Int}}}(undef, chunks...)
    for cidx in CartesianIndices(tuple(chunks...))
        idxs[cidx.I...] = ntuple(i -> (cuts[i][cidx[i]]:cuts[i][cidx[i] + 1] - 1), n)
    end
    return (idxs, cuts)
end

function localpartindex(pids::Array{Int})
    mi = myid()
    for i = 1:length(pids)
        if pids[i] == mi
            return i
        end
    end
    return 0
end
localpartindex(d::DArray) = localpartindex(procs(d))

"""
    localpart(d::DArray)

Get the local piece of a distributed array.
Returns an empty array if no local part exists on the calling process.

d[:L], d[:l], d[:LP], d[:lp] are an alternative means to get localparts.
This syntaxt can also be used for assignment. For example,
`d[:L]=v` will assign `v` to the localpart of `d`.
"""
function localpart(d::DArray{T,N,A}) where {T,N,A}
    lpidx = localpartindex(d)
    if lpidx == 0
        return empty_localpart(T,N,A)::A
    end

    return d.localpart::A
end

localpart(d::DArray, localidx...) = localpart(d)[localidx...]

_localindex(i::Integer, offset) = i - offset
_localindex(i::AbstractRange, offset) = (first(i)-offset):step(i):(last(i)-offset)
_localindex(i::AbstractUnitRange, offset) = (first(i)-offset):(last(i)-offset)

"""
    makelocal(A::DArray, I...)

Equivalent to `Array(view(A, I...))` but optimised for the case that the data is local.
Can return a view into `localpart(A)`
"""
function makelocal(A::DArray{<:Any, <:Any, AT}, I::Vararg{Any, N}) where {N, AT}
    Base.@_inline_meta
    J = map(i->Base.unalias(A, i), to_indices(A, I))
    J = map(j-> isa(j, Base.Slice) ? j.indices : j, J)
    @boundscheck checkbounds(A, J...)

    lidcs = localindices(A)
    if Base.checkbounds_indices(Bool, lidcs, J)
        # data we want is local
        viewidcs = ntuple(i -> _localindex(J[i], first(lidcs[i]) - 1), ndims(A))
        view(localpart(A), viewidcs...)
    else
        # Make more efficient (?maybe) by allocating new memory
        # only for the remote part
        viewidcs = ntuple(i -> _localindex(J[i], 0), ndims(A))
        arr = similar(AT, map(length, viewidcs)...)
        copyto!(arr, view(A, viewidcs...))
    end
end

# shortcut to set/get localparts of a distributed object
function Base.getindex(d::DArray, s::Symbol)
    @assert s in [:L, :l, :LP, :lp]
    return localpart(d)
end

function Base.setindex!(d::DArray{T,N,A}, new_lp::A, s::Symbol) where {T,N,A}
    @assert s in [:L, :l, :LP, :lp]
    d.localpart = new_lp
    new_lp
end

# fetch localpart of d at pids[i]
Base.fetch(d::DArray{T,N,A}, i) where {T,N,A} = remotecall_fetch(localpart, d.pids[i], d)

"""
    localindices(d)

A tuple describing the indices owned by the local process.
Returns a tuple with empty ranges if no local part exists on the calling process.
"""
function localindices(d::DArray)
    lpidx = localpartindex(d)
    if lpidx == 0
        return ntuple(i -> 1:0, ndims(d))
    end
    return d.indices[lpidx]
end

# Equality
function Base.:(==)(d::DArray{<:Any,<:Any,A}, a::AbstractArray) where A
    if size(d) != size(a)
        return false
    else
        b = asyncmap(procs(d)) do p
            remotecall_fetch(p) do
                localpart(d) == A(a[localindices(d)...])
            end
        end
        return all(b)
    end
end
function Base.:(==)(d::SubDArray, a::AbstractArray)
    cd = copy(d)
    t = cd == a
    close(cd)
    return t
end
Base.:(==)(a::AbstractArray, d::DArray) = d == a
Base.:(==)(a::AbstractArray, d::SubDArray) = d == a
Base.:(==)(d1::DArray, d2::DArray) = invoke(==, Tuple{DArray, AbstractArray}, d1, d2)
function Base.:(==)(d1::SubDArray, d2::DArray)
    cd1 = copy(d1)
    t = cd1 == d2
    close(cd1)
    return t
end
function Base.:(==)(d1::DArray, d2::SubDArray)
    cd2 = copy(d2)
    t = d1 == cd2
    close(cd2)
    return t
end
function Base.:(==)(d1::SubDArray, d2::SubDArray)
    cd1 = copy(d1)
    t = cd1 == d2
    close(cd1)
    return t
end

"""
    locate(d::DArray, I::Int...)

Determine the index of `procs(d)` that hold element `I`.
"""
function locate(d::DArray, I::Int...)
    ntuple(ndims(d)) do i
        fi = searchsortedlast(d.cuts[i], I[i])
        if fi >= length(d.cuts[i])
            throw(ArgumentError("element not contained in array"))
        end
        return fi
    end
end

chunk(d::DArray{T,N,A}, i...) where {T,N,A} = remotecall_fetch(localpart, d.pids[i...], d)::A

## convenience constructors ##

"""
     dzeros(dims, ...)

Construct a distributed array of zeros.
Trailing arguments are the same as those accepted by `DArray`.
"""
dzeros(dims::Dims, args...) = DArray(I->zeros(map(length,I)), dims, args...)
dzeros(::Type{T}, dims::Dims, args...) where {T} = DArray(I->zeros(T,map(length,I)), dims, args...)
dzeros(::Type{T}, d1::Integer, drest::Integer...) where {T} = dzeros(T, convert(Dims, tuple(d1, drest...)))
dzeros(d1::Integer, drest::Integer...) = dzeros(Float64, convert(Dims, tuple(d1, drest...)))
dzeros(d::Dims) = dzeros(Float64, d)

"""
    dones(dims, ...)

Construct a distributed array of ones.
Trailing arguments are the same as those accepted by `DArray`.
"""
dones(dims::Dims, args...) = DArray(I->ones(map(length,I)), dims, args...)
dones(::Type{T}, dims::Dims, args...) where {T} = DArray(I->ones(T,map(length,I)), dims, args...)
dones(::Type{T}, d1::Integer, drest::Integer...) where {T} = dones(T, convert(Dims, tuple(d1, drest...)))
dones(d1::Integer, drest::Integer...) = dones(Float64, convert(Dims, tuple(d1, drest...)))
dones(d::Dims) = dones(Float64, d)

"""
     dfill(x, dims, ...)

Construct a distributed array filled with value `x`.
Trailing arguments are the same as those accepted by `DArray`.
"""
dfill(v, dims::Dims, args...) = DArray(I->fill(v, map(length,I)), dims, args...)
dfill(v, d1::Integer, drest::Integer...) = dfill(v, convert(Dims, tuple(d1, drest...)))

"""
     drand(dims, ...)

Construct a distributed uniform random array.
Trailing arguments are the same as those accepted by `DArray`.
"""
drand(r, dims::Dims, args...) = DArray(I -> rand(r, map(length,I)), dims, args...)
drand(r, d1::Integer, drest::Integer...) = drand(r, convert(Dims, tuple(d1, drest...)))
drand(d1::Integer, drest::Integer...) = drand(Float64, convert(Dims, tuple(d1, drest...)))
drand(d::Dims, args...)  = drand(Float64, d, args...)

"""
     drandn(dims, ...)

Construct a distributed normal random array.
Trailing arguments are the same as those accepted by `DArray`.
"""
drandn(dims::Dims, args...) = DArray(I->randn(map(length,I)), dims, args...)
drandn(d1::Integer, drest::Integer...) = drandn(convert(Dims, tuple(d1, drest...)))

## conversions ##

"""
     distribute(A[; procs, dist])

Convert a local array to distributed.

`procs` optionally specifies an array of process IDs to use. (defaults to all workers)
`dist` optionally specifies a vector or tuple of the number of partitions in each dimension
"""
function distribute(A::AbstractArray;
    procs = workers()[1:min(nworkers(), maximum(size(A)))],
    dist = defaultdist(size(A), procs))
    np = prod(dist)
    procs_used = procs[1:np]
    idxs, _ = chunk_idxs([size(A)...], dist)

    s = verified_destination_serializer(reshape(procs_used, size(idxs)), size(idxs)) do pididx
        A[idxs[pididx]...]
    end
    return DArray(I->localpart(s), size(A), procs_used, dist)
end

"""
    distribute(A, DA)

Distribute a local array `A` like the distributed array `DA`.

"""
function distribute(A::AbstractArray, DA::DArray)
    size(DA) == size(A) || throw(DimensionMismatch("Distributed array has size $(size(DA)) but array has $(size(A))"))

    s = verified_destination_serializer(procs(DA), size(DA.indices)) do pididx
        A[DA.indices[pididx]...]
    end
    return DArray(I->localpart(s), DA)
end

DArray{T,N,S}(A::S) where {T,N,S<:AbstractArray} = distribute(convert(AbstractArray{T,N}, A))

function Array{S,N}(d::DArray{T,N}) where {S,T,N}
    a = Array{S}(undef, size(d))
    @sync begin
        for i = 1:length(d.pids)
            @async a[d.indices[i]...] = chunk(d, i)
        end
    end
    return a
end

function Array{S,N}(s::SubDArray{T,N}) where {S,T,N}
    I = s.indices
    d = s.parent
    if isa(I,Tuple{Vararg{UnitRange{Int}}}) && S<:T && T<:S && !isempty(s)
        l = locate(d, map(first, I)...)
        if isequal(d.indices[l...], I)
            # SubDArray corresponds to a chunk
            return chunk(d, l...)
        end
    end
    a = Array{S}(undef, size(s))
    copyto!(a, s)
end

function Base.copyto!(a::Array, s::SubDArray)
    N = ndims(a)
    a[[1:size(a,i) for i=1:N]...] = s
    return a
end

if VERSION < v"1.2"
    # This is an internal API that has changed
    reindex(A, I, J) = Base.reindex(A, I, J)
else
    reindex(A, I, J) = Base.reindex(I, J)
end

function DArray(SD::SubArray{T,N}) where {T,N}
    D = SD.parent
    DArray(size(SD), procs(D)) do I
        lindices = reindex(SD, SD.indices, I)
        convert(Array, D[lindices...])
    end
end

function Base.reshape(A::DArray{T,1,S}, d::Dims) where {T,S<:Array}
    if prod(d) != length(A)
        throw(DimensionMismatch("dimensions must be consistent with array size"))
    end
    return DArray(d) do I
        sz = map(length,I)
        d1offs = first(I[1])
        nd = length(I)

        B = Array{T}(undef, sz)
        nr = size(B,1)
        sztail = size(B)[2:end]

        for i=1:div(length(B),nr)
            i2 = CartesianIndices(sztail)[i]
            globalidx = [ I[j][i2[j-1]] for j=2:nd ]

            a = LinearIndices(d)[d1offs, globalidx...]

            B[:,i] = Array(A[a:(a+nr-1)])
        end
        B
    end
end

## indexing ##
const _allowscalar = Ref(true)
allowscalar(flag = true) = (_allowscalar[] = flag)
_scalarindexingallowed() = _allowscalar[] || throw(ErrorException("scalar indexing disabled"))

getlocalindex(d::DArray, idx...) = localpart(d)[idx...]
function getindex_tuple(d::DArray{T}, I::Tuple{Vararg{Int}}) where T
    chidx = locate(d, I...)
    idxs = d.indices[chidx...]
    localidx = ntuple(i -> (I[i] - first(idxs[i]) + 1), ndims(d))
    pid = d.pids[chidx...]
    return remotecall_fetch(getlocalindex, pid, d, localidx...)::T
end

function Base.getindex(d::DArray, i::Int)
    _scalarindexingallowed()
    return getindex_tuple(d, Tuple(CartesianIndices(d)[i]))
end
function Base.getindex(d::DArray, i::Int...)
    _scalarindexingallowed()
    return getindex_tuple(d, i)
end

Base.getindex(d::DArray) = d[1]
if VERSION > v"1.1-"
Base.getindex(d::SubDArray, I::Int...) = invoke(getindex, Tuple{SubArray{<:Any,N},Vararg{Int,N}} where N, d, I...)
end
Base.getindex(d::SubOrDArray, I::Union{Int,UnitRange{Int},Colon,Vector{Int},StepRange{Int,Int}}...) = view(d, I...)

function Base.isassigned(D::DArray, i::Integer...)
    try
        getindex_tuple(D, i)
        true
    catch e
        if isa(e, BoundsError) || isa(e, UndefRefError)
            return false
        else
            rethrow(e)
        end
    end
end

Base.copy(d::SubDArray) = copyto!(similar(d), d)
Base.copy(d::SubDArray{<:Any,2}) = copyto!(similar(d), d)

function Base.copyto!(dest::SubOrDArray, src::AbstractArray)
    asyncmap(procs(dest)) do p
        remotecall_fetch(p) do
            ldest = localpart(dest)
            copyto!(ldest, view(src, localindices(dest)...))
        end
    end
    return dest
end

function Base.deepcopy(src::DArray)
    dest = similar(src)
    asyncmap(procs(src)) do p
        remotecall_fetch(p) do
            dest[:L] = deepcopy(src[:L])
        end
    end
    return dest
end

# local copies are obtained by convert(Array, ) or assigning from
# a SubDArray to a local Array.

function Base.setindex!(a::Array, d::DArray,
        I::Union{UnitRange{Int},Colon,Vector{Int},StepRange{Int,Int}}...)
    n = length(I)
    @sync for i = 1:length(d.pids)
        K = d.indices[i]
        @async a[[I[j][K[j]] for j=1:n]...] = chunk(d, i)
    end
    return a
end

# We also want to optimize setindex! with a SubDArray source, but this is hard
# and only works on 0.5.

# Similar to Base.indexin, but just create a logical mask. Note that this
# must return a logical mask in order to support merging multiple masks
# together into one linear index since we need to know how many elements to
# skip at the end. In many cases range intersection would be much faster
# than generating a logical mask, but that loses the endpoint information.
indexin_mask(a, b::Number) = a .== b
indexin_mask(a, r::AbstractRange{Int}) = [i in r for i in a]
indexin_mask(a, b::AbstractArray{Int}) = indexin_mask(a, IntSet(b))
indexin_mask(a, b::AbstractArray) = indexin_mask(a, Set(b))
indexin_mask(a, b) = [i in b for i in a]

import Base: tail
# Given a tuple of indices and a tuple of masks, restrict the indices to the
# valid regions. This is, effectively, reversing Base.setindex_shape_check.
# We can't just use indexing into MergedIndices here because getindex is much
# pickier about singleton dimensions than setindex! is.
restrict_indices(::Tuple{}, ::Tuple{}) = ()
function restrict_indices(a::Tuple{Any, Vararg{Any}}, b::Tuple{Any, Vararg{Any}})
    if (length(a[1]) == length(b[1]) == 1) || (length(a[1]) > 1 && length(b[1]) > 1)
        (vec(a[1])[vec(b[1])], restrict_indices(tail(a), tail(b))...)
    elseif length(a[1]) == 1
        (a[1], restrict_indices(tail(a), b))
    elseif length(b[1]) == 1 && b[1][1]
        restrict_indices(a, tail(b))
    else
        throw(DimensionMismatch("this should be caught by setindex_shape_check; please submit an issue"))
    end
end
# The final indices are funky - they're allowed to accumulate together.
# An easy (albeit very inefficient) fix for too many masks is to use the
# outer product to merge them. But we can do that lazily with a custom type:
function restrict_indices(a::Tuple{Any}, b::Tuple{Any, Any, Vararg{Any}})
    (vec(a[1])[vec(ProductIndices(b, map(length, b)))],)
end
# But too many indices is much harder; this requires merging the indices
# in `a` before applying the final mask in `b`.
function restrict_indices(a::Tuple{Any, Any, Vararg{Any}}, b::Tuple{Any})
    if length(a[1]) == 1
        (a[1], restrict_indices(tail(a), b))
    else
        # When one mask spans multiple indices, we need to merge the indices
        # together. At this point, we can just use indexing to merge them since
        # there's no longer special handling of singleton dimensions
        (view(MergedIndices(a, map(length, a)), b[1]),)
    end
end

struct ProductIndices{I,N} <: AbstractArray{Bool, N}
    indices::I
    sz::NTuple{N,Int}
end
Base.size(P::ProductIndices) = P.sz
# This gets passed to map to avoid breaking propagation of inbounds
Base.@propagate_inbounds propagate_getindex(A, I...) = A[I...]
Base.@propagate_inbounds Base.getindex(P::ProductIndices{J,N}, I::Vararg{Int, N}) where {J,N} =
    Bool((&)(map(propagate_getindex, P.indices, I)...))

struct MergedIndices{I,N} <: AbstractArray{CartesianIndex{N}, N}
    indices::I
    sz::NTuple{N,Int}
end
Base.size(M::MergedIndices) = M.sz
Base.@propagate_inbounds Base.getindex(M::MergedIndices{J,N}, I::Vararg{Int, N}) where {J,N} =
    CartesianIndex(map(propagate_getindex, M.indices, I))
# Additionally, we optimize bounds checking when using MergedIndices as an
# array index since checking, e.g., A[1:500, 1:500] is *way* faster than
# checking an array of 500^2 elements of CartesianIndex{2}. This optimization
# also applies to reshapes of MergedIndices since the outer shape of the
# container doesn't affect the index elements themselves. We can go even
# farther and say that even restricted views of MergedIndices must be valid
# over the entire array. This is overly strict in general, but in this
# use-case all the merged indices must be valid at some point, so it's ok.
const ReshapedMergedIndices{T,N,M<:MergedIndices} = Base.ReshapedArray{T,N,M}
const SubMergedIndices{T,N,M<:Union{MergedIndices, ReshapedMergedIndices}} = SubArray{T,N,M}
const MergedIndicesOrSub = Union{MergedIndices, ReshapedMergedIndices, SubMergedIndices}
import Base: checkbounds_indices
@inline checkbounds_indices(::Type{Bool}, inds::Tuple{}, I::Tuple{MergedIndicesOrSub,Vararg{Any}}) =
    checkbounds_indices(Bool, inds, (parent(parent(I[1])).indices..., tail(I)...))
@inline checkbounds_indices(::Type{Bool}, inds::Tuple{Any}, I::Tuple{MergedIndicesOrSub,Vararg{Any}}) =
    checkbounds_indices(Bool, inds, (parent(parent(I[1])).indices..., tail(I)...))
@inline checkbounds_indices(::Type{Bool}, inds::Tuple, I::Tuple{MergedIndicesOrSub,Vararg{Any}}) =
    checkbounds_indices(Bool, inds, (parent(parent(I[1])).indices..., tail(I)...))

# The tricky thing here is that we want to optimize the accesses into the
# distributed array, but in doing so, we lose track of which indices in I we
# should be using.
#
# I’ve come to the conclusion that the function is utterly insane.
# There are *6* flavors of indices with four different reference points:
# 1. Find the indices of each portion of the DArray.
# 2. Find the valid subset of indices for the SubArray into that portion.
# 3. Find the portion of the `I` indices that should be used when you access the
#    `K` indices in the subarray.  This guy is nasty.  It’s totally backwards
#    from all other arrays, wherein we simply iterate over the source array’s
#    elements.  You need to *both* know which elements in `J` were skipped
#    (`indexin_mask`) and which dimensions should match up (`restrict_indices`)
# 4. If `K` doesn’t correspond to an entire chunk, reinterpret `K` in terms of
#    the local portion of the source array
function Base.setindex!(a::Array, s::SubDArray,
        I::Union{UnitRange{Int},Colon,Vector{Int},StepRange{Int,Int}}...)
    Inew = Base.to_indices(a, I)
    Base.setindex_shape_check(s, Base.index_lengths(Inew...)...)
    n = length(Inew)
    d = s.parent
    J = Base.to_indices(d, s.indices)
    @sync for i = 1:length(d.pids)
        K_c = d.indices[i]
        K = map(intersect, J, K_c)
        if !any(isempty, K)
            K_mask = map(indexin_mask, J, K_c)
            idxs = restrict_indices(Inew, K_mask)
            if isequal(K, K_c)
                # whole chunk
                @async a[idxs...] = chunk(d, i)
            else
                # partial chunk
                @async a[idxs...] =
                    remotecall_fetch(d.pids[i]) do
                        view(localpart(d), [K[j].-first(K_c[j]).+1 for j=1:length(J)]...)
                    end
            end
        end
    end
    return a
end

function Base.fill!(A::DArray, x)
    @sync for p in procs(A)
        @async remotecall_fetch((A,x)->(fill!(localpart(A), x); nothing), p, A, x)
    end
    return A
end

using Random

function Random.rand!(A::DArray, ::Type{T}) where T
    asyncmap(procs(A)) do p
        remotecall_wait((A, T)->rand!(localpart(A), T), p, A, T)
    end
end

max_issues_repo_name: JuliaParallel/DistributedArrays

max_forks_repo_name: JuliaParallel/DistributedArrays


text:
"""
    OpenStack{X,K}

    OpenStack{X,K}(filename, types, sizes, eachchunk, haschunks, write)

A wrapper for any stack-like opened dataSet that can be indexed
with `Symbol` keys to retrieve `AbstractArray` layers.

It is usually hidden from users, wrapped in a regular `RasterStack`
passed as the function argument in `open(stack)` when the stack is
contained in a single file..

`X` is a backend singleton like `NCDfile`, and `K` is a tuple
of `Symbol` keys.
"""
struct OpenStack{X,K,DS}
    dataset::DS
end
OpenStack{X,K}(dataset::DS) where {X,K,DS} = OpenStack{X,K,DS}(dataset)

# OpenStack has `X` and `K` parameter that is not recoverable from fields.
ConstructionBase.constructorof(::Type{<:OpenStack{X,K}}) where {X,K} = OpenStack{X,K}

dataset(os::OpenStack) = os.dataset
Base.keys(os::OpenStack{<:Any,K}) where K = K
Base.values(os::OpenStack{<:Any}) = (fs[k] for k in keys(fs))
# Indexing OpenStack returns memory-backed Raster.
Base.getindex(os::OpenStack, key::Symbol) = dataset(os)[key]

ext: jl


text:
col_1_to_30 = [ -2.5   0.5          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   3.0  -2.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
1.0  -3.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   3.0  -2.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
0.0   0.0         -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.5
0.0   0.0         -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
0.0   0.0          0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.5
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0
0.0   0.0          1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -5.0   5.0  -2.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0
0.0   0.0          0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.5  -2.0   0.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.5  -2.0  -3.5   4.5  -3.0   1.5   0.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -5.0   5.0  -2.0   0.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.5  -2.0   0.5   0.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.5  -2.0  -3.5   4.5  -3.0   1.5   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -5.0   5.0  -2.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.5  -2.0   0.5   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.5  -2.0  -3.5   0.0
0.0   4.5          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -5.25
0.0   1.8735e-16   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.75
0.0  -1.249e-15    0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0   0.0   0.0  -3.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0   0.0   2.25
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0   0.0   9.36751e-17
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0  -6.245e-16
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
0.0   0.0          0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0]

col_31_to_40 = [
    0.0          0.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
    0.0          0.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
    0.0          0.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
    0.5          0.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
    0.0          0.5           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
    0.0          0.0           0.5           0.0          0.0           0.0    0.0   0.0    0.0   0.0
    0.0          0.0           0.0           0.5          0.0           0.0    0.0   0.0    0.0   0.0
    0.0          0.0           0.0           0.0          0.5           0.0    0.0   0.0    0.0   0.0
    0.0          0.0           0.0           0.0          0.0           0.5    0.0   0.0    0.0   0.0
    0.0          0.0           0.0           0.0          0.0           0.0    0.5   0.0    0.0   0.0
    0.0          0.0           0.0           0.0          0.0           0.0    0.0   0.5    0.0   0.0
    0.0          0.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
    0.5          0.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
    0.0          0.5           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
    0.0          0.0           0.5           0.0          0.0           0.0    0.0   0.0    0.0   0.0
    0.0          0.0           0.0           0.5          0.0           0.0    0.0   0.0    0.0   0.0
    0.0          0.0           0.0           0.0          0.5           0.0    0.0   0.0    0.0   0.0
    0.0          0.0           0.0           0.0          0.0           0.5    0.0   0.0    0.0   0.0
    0.0          0.0           0.0           0.0          0.0           0.0    0.5   0.0    0.0   0.0
    0.0          0.0           0.0           0.0          0.0           0.0    0.0   0.5    0.0   0.0
    0.0          0.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
    1.0          0.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
    0.0          1.0           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
    0.0          0.0           1.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
    0.0          0.0           0.0           1.0          0.0           0.0    0.0   0.0    0.0   0.0
    0.0          0.0           0.0           0.0          1.0           0.0    0.0   0.0    0.0   0.0
    0.0          0.0           0.0           0.0          0.0           1.0    0.0   0.0    0.0   0.0
    0.0          0.0           0.0           0.0          0.0           0.0    1.0   0.0    0.0   0.0
    0.0          0.0           0.0           0.0          0.0           0.0    0.0   1.0    0.0   3.0
   -3.0          0.75          0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
   -3.0         -0.75          0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
    7.5         -7.5           0.0           0.0          0.0           0.0    0.0   0.0    0.0   0.0
   -4.5          6.75         -5.25         -3.0          0.75          0.0    0.0   0.0    0.0   0.0
   -1.8735e-16   2.81025e-16   0.75         -3.0         -0.75          0.0    0.0   0.0    0.0   0.0
    1.249e-15   -1.8735e-15   -3.0           7.5         -7.5           0.0    0.0   0.0    0.0   0.0
    0.0          0.0           2.25         -4.5          6.75         -5.25  -3.0   0.75   0.0   0.0
    0.0          0.0           9.36751e-17  -1.8735e-16   2.81025e-16   0.75  -3.0  -0.75   0.0   0.0
    0.0          0.0          -6.245e-16     1.249e-15   -1.8735e-15   -3.0    7.5  -7.5    0.0   0.0
    0.0          0.0           0.0           0.0          0.0           0.25  -0.5   0.75  -2.5   2.0
    0.0          0.0           0.0           0.0          0.0           0.5   -1.0   1.5    1.0  -2.0
]

B_data = [col_1_to_30 col_31_to_40]


text:
############################################################################
### Types
############################################################################

struct SemDiffEmpty <: SemDiff end

############################################################################
### Constructor
############################################################################

# SemDiffEmpty(;kwargs...) = SemDiffEmpty()

############################################################################
### Recommended methods
############################################################################

update_observed(diff::SemDiffEmpty, observed::SemObs; kwargs...) = diff

############################################################################
### Pretty Printing
############################################################################

function Base.show(io::IO, struct_inst::SemDiffEmpty)
    StructuralEquationModels.print_type_name(io, struct_inst)
end


text:
# NOTE: unlike `Base.show` we always ask for explicit MIME type

"""
    pprint([io::IO=stdout], [mimetype], x; kw...)

Pretty print an object `x` to `io`, default `io` is `stdout`.

!!! note

    `pprint` will detect if an object type has overloaded `Base.show`,
    and use that if possible, overloading `Base.show` to `GarishPrint`
    for custom type should use [`pprint_struct`](@ref) to avoid
    recursive call into `Base.show`.

# Keyword Arguments

- `indent::Int`: indent size, default is `2`.
- `compact::Bool`: whether print withint one line, default is `get(io, :compact, false)`.
- `limit::Bool`: whether the print is limited, default is `get(io, :compact, false)`.
- `displaysize::Tuple{Int, Int}`: the displaysize hint of printed string, note this is not stricted obeyed,
default is `displaysize(io)`.
- `show_indent::Bool`: whether print indentation hint, default is `true`.
- `color::Bool`: whether print with color, default is `true`.

## Keyword Arguments for option struct defined by Configurations

- `include_defaults::Bool`: whether print the default values,
    default is `false` to provide more compact printing in REPL.

## Color Preference

color preference is available as keyword arguments to override the
default color scheme. These arguments may take any of the values
`:normal`, `:default`, `:bold`, `:black`, `:blink`, `:blue`,
`:cyan`, `:green`, `:hidden`, `:light_black`, `:light_blue`, `:light_cyan`, `:light_green`,
`:light_magenta`, `:light_red`, `:light_yellow`, `:magenta`, `:nothing`, `:red`, `:reverse`,
`:underline`, `:white`, or `:yellow` or an integer between 0 and 255 inclusive. Note that
not all terminals support 256 colors.

The default color scheme can be checked via `GarishPrint.default_colors_256()` for 256 color,
and `GarishPrint.default_colors_ansi()` for ANSI color. The 256 color will be used when
the terminal is detected to support 256 color.

- `fieldname`: field name of a struct.
- `type`: the color of a type.
- `operator`: the color of an operator, e.g `+`, `=>`.
- `literal`: the color of literals.
- `constant`: the color of constants, e.g `π`.
- `number`: the color of numbers, e.g `1.2`, `1`.
- `string`: the color of string.
- `comment`: comments, e.g `# some comments`
- `undef`: the const binding to `UndefInitializer`
- `linenumber`: line numbers.

# Notes

The color print and compact print can also be turned on/off by
setting `IOContext`, e.g `IOContext(io, :color=>false)` will print
without color, and `IOContext(io, :compact=>true)` will print within
one line. This is also what the standard Julia `IO` objects follows
in printing by default.
"""
function pprint end

pprint(x; kw...) = pprint(stdout, x; kw...)
pprint(io::IO, x;kw...) = pprint(io, MIME"text/plain"(), x; kw...)
# wrap other IO with GarishIO
pprint(io::IO, m::MIME, x; kw...) = pprint(GarishIO(io; kw...), m, x)

# Type
pprint(io::GarishIO, mime::MIME, @specialize(x::Type)) = show(io, mime, x)
pprint(io::GarishIO, ::MIME"text/plain", @specialize(x::Type)) = print_token(io, :type, x)

# Struct
function pprint(io::GarishIO, mime::MIME, @nospecialize(x))
    if fallback_to_default_show(io, x) && isstructtype(typeof(x))
        return pprint_struct(io, mime, x)
    elseif io.state.level > 0 # print show inside
        show_text_within(io, mime, x)
    else # fallback to show unless it is a struct type
        show(IOContext(io), mime, x)
    end
end

function show_text_within(io::GarishIO, mime::MIME, x)
    buf = IOBuffer()
    indentation = io.indent * io.state.level + io.state.offset
    new_displaysize = (io.displaysize[1], io.displaysize[2] - indentation)

    buf_io = GarishIO(buf, io;  limit=true, displaysize=new_displaysize, state=PrintState())
    show(buf_io, mime, x)
    raw = String(take!(buf))
    buf_lines = split(raw, '\n')
    for k in 1:length(buf_lines)-1
        line = buf_lines[k]
        if !(io.state.noindent_in_first_line && k == 1)
            print_indent(io)
        end
        println(io.bland_io, line)
    end

    # pretty print last line
    last_line = buf_lines[end]
    length(buf_lines) == 1 || print_indent(io)
    if endswith(last_line, ')')
        print(io.bland_io, last_line)
    else
        print(io.bland_io, last_line)
        if length(buf_lines) > 1
            println(io)
            print_indent(io) # force put a new line at the end
        end
    end
    return
end

function fallback_to_default_show(io::IO, x)
    # NOTE: Base.show(::IO, ::MIME"text/plain", ::Any) forwards to
    # Base.show(::IO, ::Any)

    # check if we are gonna call Base.show(::IO, ::MIME"text/plain", ::Any)
    mt = methods(Base.show, (typeof(io), MIME"text/plain", typeof(x)))
    length(mt.ms) == 1 && any(mt.ms) do method
        method.sig == Tuple{typeof(Base.show), IO, MIME"text/plain", Any}
    end || return false

    # check if we are gonna call Base.show(::IO, ::Any)
    mt = methods(Base.show, (typeof(io), typeof(x)))
    length(mt.ms) == 1 && return any(mt.ms) do method
        method.sig == Tuple{typeof(Base.show), IO, Any}
    end
end


text:
#READ DATA FIRST
read_vertex("./smallest/khazanah.csv", "./smallest/permintaan.csv")
read_edge("./smallest/trayek.csv")
read_mode("./smallest/kendaraan.csv")

#test f
f(1) == f(E(1))

#INVENTORY FEASIBILITY
@test inventory_feasibility(V())

#exhaustiveness of edge list
@test EM_check(E(),M())
@test EV_check(E(),V())

#optimality of test pack
tes_model_IP = raw_model_IP(V(),E(),M(),T())
optimize!(tes_model_IP)
@test termination_status(tes_model_IP) == MOI.OPTIMAL

ip_sol = convert_raw(tes_model_IP.obj_dict[:o], tes_model_IP.obj_dict[:p])
@test typeof(collect(keys(period_view(ip_sol)))) == Vector{Int64}
@test typeof(collect(keys(lin_view(ip_sol)))) == Vector{Int64}

#optimality of relaxation of test pack
tes_model_LP = raw_model_LP(V(),E(),M(),T())
optimize!(tes_model_LP)
@test termination_status(tes_model_LP) == MOI.OPTIMAL

lp_sol = convert_raw(tes_model_LP.obj_dict[:o], tes_model_LP.obj_dict[:p])
@test typeof(collect(keys(period_view(lp_sol)))) == Vector{Int64}
@test typeof(collect(keys(lin_view(lp_sol)))) == Vector{Int64}


text:
function softmax4d!(y::CuArray{T}, x::CuArray{T};
                 handle=cudnnhandle(),
                 algorithm=CUDNN_SOFTMAX_ACCURATE, # or CUDNN_SOFTMAX_FAST
                 mode=CUDNN_SOFTMAX_MODE_INSTANCE, # or CUDNN_SOFTMAX_MODE_CHANNEL
                 alpha=1.0, beta=0.0) where T
    @cuda(cudnn, cudnnSoftmaxForward,
          (Cptr, Cuint, Cuint, Ptr{T}, Cptr, Ptr{T}, Ptr{T}, Cptr, Ptr{T}),
          handle, algorithm, mode, Ref(T(alpha)), TD(x), x.ptr, Ref(T(beta)), TD(y), y.ptr)
    return y
end

softmax4d(x::CuArray{T}) where T = softmax4d!(similar(x), x)

function softmax4d_grad!(dx::CuArray{T}, y::CuArray{T}, dy::CuArray{T};
                         handle=cudnnhandle(),
                         algorithm=CUDNN_SOFTMAX_ACCURATE, # or CUDNN_SOFTMAX_FAST
                         mode=CUDNN_SOFTMAX_MODE_INSTANCE, # or CUDNN_SOFTMAX_MODE_CHANNEL
                         alpha=1.0, beta=0.0) where T
    @cuda(cudnn, cudnnSoftmaxBackward,
          (Cptr, Cuint, Cuint,
           Ptr{T}, Cptr, Ptr{T},
           Cptr, Ptr{T},
           Ptr{T}, Cptr, Ptr{T}),
          handle, algorithm, mode,
          Ref(T(alpha)), TD(y), y.ptr,
          TD(dy), dy.ptr,
          Ref(T(beta)), TD(dx), dx.ptr)
    return dx
end

softmax4d_grad(y::CuArray{T}, dy::CuArray{T}) where T = softmax4d_grad!(similar(dy), y, dy)


text:
import LinearAlgebra: ldiv!, \

export Identity

#### Type-handling
"""
    Adivtype(A, b)
Determine type of the division of an element of `b` against an element of `A`:
`typeof(one(eltype(b))/one(eltype(A)))`
"""
Adivtype(A, b) = typeof(one(eltype(b))/one(eltype(A)))

"""
    zerox(A, b)
Build a zeros vector `Vector{T}`, where `T` is `Adivtype(A,b)`.
"""
zerox(A, b) = zeros(Adivtype(A, b), size(A, 2))

"""
No-op preconditioner
"""
struct Identity end

\(::Identity, x) = copy(x)
ldiv!(::Identity, x) = x
ldiv!(y, ::Identity, x) = copyto!(y, x)


text:
# Error type and decoding functionality

export VEError

"""
    VEError(code)
    VEError(code, meta)

Create a VEDA error object with error code `code`. The optional `meta` parameter indicates
whether extra information, such as error logs, is known.
"""
struct VEError <: Exception
    code::VEDAresult
    meta::Any

    VEError(code, meta=nothing) = new(code, meta)
end

Base.convert(::Type{VEDAresult}, err::VEError) = err.code

Base.:(==)(x::VEError,y::VEError) = x.code == y.code

"""
    name(err::VEError)

Gets the string representation of an error code.

```jldoctest
julia> err = VEError(VEDA.VEDAresult_enum(1))
VEError(VEDA_ERROR_INVALID_VALUE)

julia> name(err)
"ERROR_INVALID_VALUE"
```
"""
function name(err::VEError)
    str_ref = Ref{Cstring}()
    vedaGetErrorName(err, str_ref)
    unsafe_string(str_ref[])[6:end]
end

"""
    description(err::VEError)

Gets the string description of an error code.
"""
function description(err::VEError)
    if err.code == -1%UInt32
        "Cannot use the VEDA stub libraries."
    else
        str_ref = Ref{Cstring}()
        vedaGetErrorString(err, str_ref)
        unsafe_string(str_ref[])
    end
end

function Base.showerror(io::IO, err::VEError)
    try
        print(io, "VEDA error: $(description(err)) (code $(reinterpret(Int32, err.code)), $(name(err)))")
    catch
        # we might throw before the library is initialized
        print(io, "VEDA error (code $(reinterpret(Int32, err.code)), $(err.code))")
    end

    if err.meta != nothing
        print(io, "\n")
        print(io, err.meta)
    end
end

Base.show(io::IO, ::MIME"text/plain", err::VEError) = print(io, "VEError($(err.code))")

@enum_without_prefix VEDAresult_enum VEDA_

## API call wrapper

#@inline function initialize_context()
#    prepare_cuda_state()
#    return
#end

# outlined functionality to avoid GC frame allocation
@noinline function throw_api_error(res)
    #if res == ERROR_OUT_OF_MEMORY
    #    throw(OutOfGPUMemoryError())
    #else
        throw(VEError(res))
    #end
end

macro check(ex)
    quote
        res = $(esc(ex))
        if res != SUCCESS
            throw_api_error(res)
        end
        nothing
    end
end


text:
using Documenter, ClusterTrees
# using DocumenterMarkdown
makedocs(sitename="ClusterTrees Documentation", clean=false)

deploydocs(
    repo = "github.com/krcools/ClusterTrees.jl.git",
)


text:
function B_P(cur_reactor::AbstractReactor)
  cur_B = (2/5) * cur_reactor.pi

  cur_B *= cur_reactor.I_P

  cur_B /= perimeter(cur_reactor)

  cur_B
end

function b_p(cur_gamma::SymEngine.Basic, cur_rho::SymEngine.Basic, use_default::Bool=false)
  use_default && return b_p_default(cur_gamma, cur_rho)
  symbols(:b_theta)
end

function b_p(cur_gamma::AbstractFloat, cur_rho::AbstractFloat)
  b_p_default(cur_gamma, cur_rho)
end

function b_p_default(cur_gamma::AbstractSymbol, cur_rho::AbstractSymbol)
  cur_repeat = 1 + cur_gamma

  cur_b = cur_gamma * cur_rho ^ 2
  cur_b -= cur_repeat

  cur_b *= -exp( cur_gamma * cur_rho ^ 2 )
  cur_b -= cur_repeat

  cur_b /= cur_rho
  cur_b /= exp(cur_gamma) - cur_repeat

  cur_b
end

function int_b_p(cur_gamma::AbstractFloat)
  isfinite(cur_gamma) || return NaN

  cur_func = (cur_rho) -> cur_rho * b_p(cur_gamma, cur_rho) ^ 2

  cur_integral = norm_int(cur_func)

  cur_integral
end


text:
using Test
using NeuralThermalODE

# Test data file.
# The CSV format is:
# | Time | Temperature |
# ----------------------
# | INT  | FLOAT       |
# | ...  | ...         |
test_data = read_data("test_data/test_data.csv")

# Domain.
const N = 10
const tmat = 0.00635 # Steel thickness.
const dx = tmat / N  # Element size.

# Boundary conditions from the test data.
const sinkT = test_data[end, 2] # Sink temperature on exposed side.
const ambT = test_data[1, 2]    # Sink temperature on unexposed side.
u20 = ambT * ones(N)            # Initial temperature.

const tspan = (Float64(test_data[1, 1]), Float64(test_data[2, 1]))

# For unexposed HTC formula for ambient see:
# /home/nik/Dropbox/PhD/Academic/General_Data/HTC_effect/Unexposed
# This might need to be investigated to see if 4.0 is good enough.
const pr = BuildThermalODE(dx, 450.0, 7850.0, 42.0, sinkT, ambT, 4.0)

# Solve ODE.
p = read_checkpoint("test_data/initial_guess.csv")

# Pass additional arguments through closure.
prob = ODEProblem((du, u20, p, t) -> heat_transfer(pr, du, u20, p, t), u20, tspan, p, saveat = 1.0)

data = Iterators.repeated((), 5)
opt = ADAM(20.0)

target_loss = 10.0

# Track loss function.
tot_loss = []

training_cnt = -1
early_exit = false

algo = Tsit5()

cb = function ()
    global training_cnt += 1

    # Save loss function.
    loss = loss_rd(Vector(test_data[2]), p, prob, algo)
    append!(tot_loss, loss)

    if loss <= target_loss
        early_exit = true
        try
            Flux.stop()
        catch
        end
    end

    sol = solve(remake(prob, p = Flux.data(p)), algo, saveat = 1.0)
end

cb()
Flux.train!(() -> loss_rd(Vector(test_data[2]), p, prob, algo), [p], data, opt, cb = cb)

res = Flux.data(p)[1:end-1]
@test isapprox(res[1], 706.0, atol = 1.0)

max_stars_repo_name: NikStoyanov/NeuralThermalODE

max_issues_repo_name: NikStoyanov/NeuralThermalODE

max_forks_repo_name: NikStoyanov/NeuralThermalODE


text:
include("evaluate_lambda.jl")

using ApproxFun, PyPlot

x = Fun(Segment{BigFloat}())

αBF = big(-0.5)
βBF = big(0.0)

δ = 1.0
α = Float64(αBF)
ℓmax = 1000

λ = zeros(Float64, ℓmax)
λfast = zeros(Float64, ℓmax)
λBF = zeros(BigFloat, ℓmax)

for ℓ in 1:ℓmax
    λ[ℓ] = evaluate_lambda_rec(ℓ, δ, α)
    λfast[ℓ] = evaluate_lambda_asy(ℓ, δ, α)
    Pℓ = legendre(BigFloat, ℓ)
    f = Pℓ(1-big(δ)^2/2*(1-x)/2)
    g = (f-1)/(1-x)
    #w = Fun(JacobiWeight(big(0.0), big(α), space(x)),[big(1.0)])
    #λBF[ℓ] = sum(g*w)
    n = ncoefficients(g)
    c = zeros(BigFloat, n)
    αG = αBF
    if space(g) == space(x)
        αG = αBF
    else
        αG = αBF-1
    end
    c[1] = big(2.0)^(βBF+αG+1)*gamma(βBF+1)*gamma(αG+1)/gamma(βBF+αG+2)
    if n > 1
        c[2] = c[1]*(βBF-αG)/(βBF+αG+2)
        for i=1:n-2
            c[i+2] = (2(βBF-αG)*c[i+1]-(βBF+αG-i+2)*c[i])/(βBF+αG+i+2)
        end
    end
    λBF[ℓ] = ApproxFun.dotu(g.coefficients, c)
    println("Finished ℓ: ",ℓ)
end

λBF64 = Float64.(λBF)

clf()
semilogy(1:ℓmax, abs.((λ-λBF64)./λBF64), ".k", 1:ℓmax, abs.((λfast-λBF64)./λBF64), "xk")
xlabel("\$\\ell\$"); ylabel("Relative Error"); grid()
legend(["REC","ASY"])
savefig("spectrum_error.pdf")

L = [collect(1:10); round.([Int], logspace(0,4,80))[22:end]]
tREC = zeros(Float64, length(L))
tASY = zeros(Float64, length(L))

j = 1
for ℓ in L
    tREC[j] = @elapsed evaluate_lambda_rec(ℓ, δ, α)
    tASY[j] = @elapsed evaluate_lambda_asy(ℓ, δ, α)
    println("Finished ℓ: ",ℓ)
    gc()
    sleep(0.1)
    j += 1
end

clf()
loglog(L, tREC, ".k", L, tASY, "xk")
loglog(L, 6.25e-8*L.*log.(L), "-k")
loglog(L, 1e-8*L.^2, "--k")
ylim(1e-5,1e0)
xlabel("\$\\ell\$"); ylabel("Execution Time (s)"); grid()
legend(["REC","ASY","\$\\mathcal{O}(\\ell\\log\\ell)\$","\$\\mathcal{O}(\\ell^2)\$"])
savefig("spectrum_time.pdf")


text:
# Author: Gabrio Rizzuti, rizzuti.gabrio@gatech.edu
# Date: August, 2020
# Copyright: Georgia Institute of Technology, 2020
#
# Run VI via invertible networks to sample from Gaussian mixture in 2D

## Module loading

using DrWatson
@quickactivate "UncertaintyQuantificationAVP"
using UncertaintyQuantificationAVP
using Flux.Optimise, InvertibleNetworks
using LinearAlgebra, PyPlot, Statistics
using Random; seed = 1; Random.seed!(seed)

## Setting negative log density for Gaussian mixture

# Gaussian modes parameters
nmodes = 20
R = 1f0
muComplex = Complex{Float32}.(R*exp.(1im*2f0*pi*collect(range(0f0, 1f0; length=nmodes+5))))
mu = [real.(muComplex[1:end-5]) imag.(muComplex[1:end-5])]
sigma2 = 0.025f0*ones(Float32, nmodes)

# Sampling function
sample(n) = sample_gaussmixture(n, mu, sigma2)
Xtrue = sample(1024)

# Negative log density
ε = 1f-10
neglogden(X) = negLogDensity_gaussmixture(X, mu, sigma2; ε = ε)

## Setting loss function

function loss(Z::Array{Float32, 4})
    X, logdet = T.forward(Z)
    nlogp, ΔX = neglogden(X)
    f = sum(nlogp)/batchsize-logdet
    T.backward(ΔX/batchsize, X)
    return f
end

## Network initialization

# Generator
depth = 20
n_hidden = 32
batchsize = 2^6
net_type = NetworkHINT1D
T = net_type(2, batchsize, n_hidden, depth; permute="full")
# net_type = NetworkGlow1D
# T = net_type(2, batchsize, n_hidden, depth)
initialize_id!(T)

## Training

# Training data
trainsize = 2^11
Ztrain = randn(Float32, 1, 1, 2, trainsize)

# N of iterations
nepochs = 2^7

# Stepsize
lr = 1f-3
decay = 0.9f0
nbatches = Int64(round(trainsize/batchsize))
decay_step = nbatches*2^4
# decay_step = nbatches*Int64(round(nepochs/log(1f-2)*log(decay)))
clip = 0f0
opt = Optimiser(ExpDecay(lr, decay, decay_step, clip), ADAM(lr))

# Train model
fval = Array{Float32, 1}(undef, 0)
training!(T, fval, loss, Ztrain, nepochs, batchsize, opt; gradclip=true)

## Testing

testsize = trainsize
Ztest = randn(Float32, 1, 1, 2, testsize)
Xtest, _ = T.forward(Ztest)

## Saving results

params_exp = Dict(
    :seed => seed, #seed
    :nmodes => nmodes, :mu => mu, :sigma2 => sigma2, # Gaussian mixture pars
    :batchsize => batchsize, :trainsize => trainsize, :Ztrain => Ztrain, :nepochs => nepochs, # training pars
    :lr => lr, :decay => decay, :decay_step => decay_step, :clip => clip, # optimizer pars
    :net_type => string(net_type), :depth => depth, :n_hidden => n_hidden) # network arch
results_exp = copy(params_exp)
results_exp[:fval] = fval
results_exp[:net_params] = get_params(T)
results_exp[:Ztest] = Ztest
results_exp[:Xtest] = Xtest
results_exp[:Xtrue] = Xtrue
wsave(datadir("GaussMixture", savename("VI", params_exp, "bson")), results_exp)

# ## Plotting
#
# figure()
# plot(Xtrue[1, 1, 1, :], Xtrue[1, 1, 2, :], ".")
# plot(Xtest[1, 1, 1, :], Xtest[1, 1, 2, :], "r*")
# title("VI samples")


text:
using Statistics: mean
using ForwardDiff
using Dierckx

function toDict(res)
    return Dict(
        fieldnames(typeof(res)) .=> getfield.(Ref(res), fieldnames(typeof(res)))
    )
end

function naturalsort(x::Vector{String})
    f = text -> all(isnumeric, text) ? Char(parse(Int, text)) : text
    sorter = key -> join(f(c) for c in eachmatch(r"[0-9]+|[^0-9]+", key))
    sort(x, by=sorter)
end

# ------------------------------ interpolations ------------------------------ #

"""
Simple linear interpolation
"""
# ξ(x) = interpolate(x, BSpline(Quadratic(Reflect(OnCell()))))
ξ(x) = interpolate(x, BSpline(Linear()))

"""
  interpolate_wrt_to(x, y)

Interpolate `y` with respect to `x`.

Given arrays `a=1:101` and `time=0:100`, `ξ(a)` gives an interpolation
object that can only be queried by values ∈ [1, 101] (the indices range of `a`).
`interpolate_wrt_to(time, a)` on the other hand gives an interpolation object that
can be queried by values ∈ [0, 100], the values ranges of `time`.
"""
function interpolate_wrt_to(x, y)
    @assert length(x) == length(y) "The x and data array should have the same length"
    return Spline1D(x, y; w=ones(length(x)), k=1, bc="nearest", s=0.0)
end

"""
Upsamble a set of variables through interpolation
"""
function upsample(data...; δp=0.001)
    n = length(data)
    P = range(0, 1, length(data[1]))

    # upsample data
    itp = Interpolations.scale(
        interpolate(hcat(data...), (BSpline(Cubic(Natural(OnGrid()))), NoInterp())), P, 1:n
    )

    tfine = 0:δp:1
    upsampled = []
    for i in 1:n
        push!(upsampled, [itp(t, i) for t in tfine])
    end
    return upsampled
end

# -------------------------------- derivatives ------------------------------- #
"""
Interpolate an array and get the derivative at
time values.

From: https://discourse.julialang.org/t/differentiation-without-explicit-function-np-gradient/57784
"""
function ∂(time, x)
    itp = interpolate((time,), x, Gridded(Linear()))
    return map((t) -> ForwardDiff.derivative(itp, t), time)
end

"""
Prepend a 0 to a vector, useful when takin diff
but want the resulting vector to have the right length
"""
ρ(x) = vcat([0], x)

"""
Shorthand for diff givin a vector of right length
"""
Δ(x) = ρ(diff(x))

# ----------------------------------- misc ----------------------------------- #
"""
Round a Float and turn it into an integer
"""
int(x::Number) = (Int ∘ round)(x)

"""
Unwrap circular variables (in radians)
"""
function unwrap(v, inplace=false)
    # currently assuming an array
    unwrapped = inplace ? v : copy(v)
    for i in 2:length(v)
        while unwrapped[i] - unwrapped[i - 1] >= π
            unwrapped[i] -= 2π
        end
        while unwrapped[i] - unwrapped[i - 1] <= -π
            unwrapped[i] += 2π
        end
    end
    return unwrapped
end

unwrap!(v) = unwrap(v, true)

# --------------------------------- smoothing -------------------------------- #
"""
Moving average with window size `n`
"""
function movingaverage(g, n)
    return [i < n ? mean(g[begin:i]) : mean(g[(i - n + 1):i]) for i in 1:length(g)]
end

# --------------------------------- geometry --------------------------------- #

"""
  closest_point_idx(X::Vector{Number}, x::Number, Y::Vector{Number}, y::Number)

Get the index of the point in a vector (X, Y) closest to a point (x,y)
"""
function closest_point_idx(X::Vector{T}, x::T, Y::Vector{T}, y::T) where {T<:Number}
    return argmin(sqrt.((X .- x) .^ 2 .+ (Y .- y) .^ 2))
end

euclidean(x0, x1, y0, y1) = √((x0 - x1)^2 + (y0 - y1)^2)

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
#function image_func(M::Map)
#  return M.header.image
#end
#
#function preim_func(M::Map)
#  return M.header.preim
#end

function show(io::IO, M::Map)
  print(io, "Map with following data\n")
  print(io, "Domain:\n")
  print(io, "=======\n")
  print(io, "$(domain(M))\n\n")
  print(io, "Codomain:\n")
  print(io, "=========\n")
  print(io, "$(codomain(M))\n")
end

function preimage(M::Map, a)
  if isdefined(M.header, :preimage)
    p = M.header.preimage(a)#::elem_type(D)
    @assert parent(p) == domain(M)
    return p
  end
  error("No preimage function known")
end

function image{D, C}(M::Map{D, C}, a)
  if isdefined(M.header, :image)
    return M.header.image(a)::elem_type(C)
  else
    error("No image function known")
  end
end

function show(io::IO, M::InverseMap)
  println(io, "inverse of")
  println(io, " ", M.origin)
end

function inv(a::Map)
  return InverseMap(a)
end

function show(io::IO, M::CoerceMap)
  println(io, "Coerce: $(domain(M)) -> $(codomain(M))")
end

\(f::Map, x) = preimage(f, x)


text: Dict{String,Any}("answer8" => Dict{String,Any}("value" => "δ","type" => "string"),"answer4" => Dict{String,Any}("value" => "δ","type" => "string"))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
#!/usr/bin/env julia

using Weber

version = v"0.0.1"
sid,skip = @read_args("Duration Discrimination ($version).")

#======================================================================
Experiment Setup
======================================================================#
# 2AFC duration discrimination. Each stimulus consists of 2 short tones
# separated by given duration. The task is to indicate which of 2 intervals
# is longer.

const atten_dB = 30 # adjust to calibrate levels
const n_trials_per_block = 60
const n_blocks = 6

const tone_freq = 1kHz
const tone_length = 10ms
const standard_length = 100ms # from start of tone 1 to start of tone 2

const SOA = 900ms # how long to wait to starting play stimulus 2
const response_delay = 300ms # how long to wait before asking for a response
const trial_spacing = 500ms # how long to wait between trials

const adapter = levitt_adapter(first_delta=2,down=3,up=1,
                               big=0.5,little=0.05,min_delta=0,
                               max_delta=(SOA/2-tone_length)/tone_length)

#======================================================================#

experiment = Experiment(
  skip=skip,
  columns = [
    :sid => sid,
    :condition => "example",
    :version => version,
    :standard => standard_length,
    :block
  ]
)

isresponse(e) = iskeydown(e,key"p") || iskeydown(e,key"q")

beep = attenuate(ramp(tone(tone_freq,tone_length)),atten_dB)
beep_beep(len) = [beep;silence(len - tone_length);beep]

function trial_stimuli(first_longer,delta)
  standard = beep_beep(standard_length)
  signal = beep_beep(standard_length*(1+delta))
  if first_longer
    [signal;silence(SOA-duration(signal));standard]
  else
    [standard;silence(SOA-duration(standard));signal]
  end
end

function duration2AFC(adapter;keys...)
  first_longer = rand(Bool)
  resp = response(adapter,key"q" => "first_longer",key"p" => "second_longer",
                  correct=(first_longer? "first_longer" : "second_longer");
                  keys...)

  stim() = trial_stimuli(first_longer,delta(adapter))
  [moment(trial_spacing,play,stim),show_cross(),
   moment(SOA + standard_length + response_delay,display,
          "Was the first [Q] or second sound [P] longer?"),
   resp,await_response(isresponse)]
end

setup(experiment) do
  addbreak(moment(record,"start"),
           moment(play,@> tone(1kHz,1s) ramp attenuate(atten_dB)),
           moment(1s))

  addbreak(instruct("""

      On each trial, you will hear two sounds. Indicate which of the two sounds
  you heard as longer. Hit 'Q' if the first sound was longer, and 'P' if the
  second was longer.

  """))

  for block in 1:n_blocks
    @addtrials let adapter = adapter
      for trial in 1:n_trials_per_block
        addtrial(duration2AFC(adapter,block=block))
      end
      # define this string during experiment setup
      # when we know what block we're on...
      heading = "Block $block of $n_blocks"

      function report_threshold()
        mean,sd = estimate(adapter)
        thresh = round(mean,3)*standard_length
        thresh_sd = round(sd,3)*standard_length

        # define this string during run time when we know
        # what the threshold estimate is.
        "$heading: Threshold $thresh (SD: $thresh_sd)\n"*
        "Hit spacebar to continue..."
      end

      addbreak(moment(display,report_threshold,clean_whitespace=false),
               await_response(iskeydown(key":space:")))
    end
  end
end

# play a test tone, to verify levels
run(experiment)


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

mutable struct IoK8sApiCoreV1LimitRangeItem <: SwaggerModel
    default::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: default
    defaultRequest::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: defaultRequest
    max::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: max
    maxLimitRequestRatio::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: maxLimitRequestRatio
    min::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: min
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiCoreV1LimitRangeItem(;default=nothing, defaultRequest=nothing, max=nothing, maxLimitRequestRatio=nothing, min=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1LimitRangeItem, Symbol("default"), default)
        setfield!(o, Symbol("default"), default)
        validate_property(IoK8sApiCoreV1LimitRangeItem, Symbol("defaultRequest"), defaultRequest)
        setfield!(o, Symbol("defaultRequest"), defaultRequest)
        validate_property(IoK8sApiCoreV1LimitRangeItem, Symbol("max"), max)
        setfield!(o, Symbol("max"), max)
        validate_property(IoK8sApiCoreV1LimitRangeItem, Symbol("maxLimitRequestRatio"), maxLimitRequestRatio)
        setfield!(o, Symbol("maxLimitRequestRatio"), maxLimitRequestRatio)
        validate_property(IoK8sApiCoreV1LimitRangeItem, Symbol("min"), min)
        setfield!(o, Symbol("min"), min)
        validate_property(IoK8sApiCoreV1LimitRangeItem, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiCoreV1LimitRangeItem

const _property_map_IoK8sApiCoreV1LimitRangeItem = Dict{Symbol,Symbol}(Symbol("default")=>Symbol("default"), Symbol("defaultRequest")=>Symbol("defaultRequest"), Symbol("max")=>Symbol("max"), Symbol("maxLimitRequestRatio")=>Symbol("maxLimitRequestRatio"), Symbol("min")=>Symbol("min"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiCoreV1LimitRangeItem = Dict{Symbol,String}(Symbol("default")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}", Symbol("defaultRequest")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}", Symbol("max")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}", Symbol("maxLimitRequestRatio")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}", Symbol("min")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1LimitRangeItem }) = collect(keys(_property_map_IoK8sApiCoreV1LimitRangeItem))
Swagger.property_type(::Type{ IoK8sApiCoreV1LimitRangeItem }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1LimitRangeItem[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1LimitRangeItem }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1LimitRangeItem[property_name]

function check_required(o::IoK8sApiCoreV1LimitRangeItem)
    true
end

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


text:
using StructJuMP

######
# DATA
######
NUMSTAGES = 3
TESTTIME = GEN = WIND = BUS = LIN = GEN = LOAD = 1:2
NODES = HORIZON = [1]
FIXEDGENCOST = THETASCALE = COSTSCALE = 1
Pmax = zeros(2,2)
np_cap = wind_share = X = V = max_ur = Pgen_init = gen_cost = ones(2)
rec_bus = snd_bus = bus_gen = bus_load = bus_wind = [1:2]
wind_total = [1 1]
node = NODES[1]
Pload = ones(2,2)
ref_bus = 1

m = StructuredModel()

@defStochasticVar(m, dummyV >= 0)

s0 = StructuredModel(parent=m)
@defStochasticVar(s0, 0 <= y[TESTTIME,GEN] <= 1)
@setObjective(s0, Min, FIXEDGENCOST*sum{y[t,j], t=TESTTIME,j=GEN})

for it in 1:NUMSTAGES
    st = StructuredModel(parent=s0)
    @defStochasticVar(st, 0 <= Pgen[TESTTIME,j=GEN] <= np_cap[j])
    @defStochasticVar(st, 0 <= Pwind[t=TESTTIME,j=WIND] <= wind_total[node,t]*wind_share[j])
    @defStochasticVar(st, -THETASCALE*π/2 <= theta[TESTTIME,j=BUS] <= THETASCALE*π/2)
    @defStochasticVar(st, -Pmax[i] <= P[TESTTIME,i=LIN] <= Pmax[i])
    for t in TESTTIME
        for j in BUS
            @addConstraint(st, ( sum{P[t,i], i=LIN; j==rec_bus[i]}
                                -sum{P[t,i], i=LIN; j==snd_bus[i]}
                                +sum{Pgen[t,i], i=GEN; j==bus_gen[i]}
                                -sum{Pload[t,i], i=LOAD; j==bus_load[i]}
                                +sum{wind_total[node,t]*wind_share[i]-Pwind[t,i], i=WIND; j==bus_wind[i]}
                                == 0))
        end
        @addConstraint(st, theta[t,ref_bus] == 0)
        for i in LIN
            @addConstraint(st, X[i]*P[t,i] - V[i]^2*(theta[t,snd_bus[i]]-theta[t,rec_bus[i]])/THETASCALE == 0)
        end
    end
    # ramp constraints
    for t in HORIZON
        for i in GEN
            @addConstraint(st, -max_ur[i] <= Pgen[t,i] - Pgen[t+1,i] <= max_ur[i]) # range constraint, beware...
        end
    end

    for i in GEN
        @addConstraint(st, -max_ur[i] <= Pgen[1,i] - Pgen_init[i] <= max_ur[i])
    end

    # linking constraints
    # gen[s,t,j] <= np_cap[j]*y[t,j]
    for t in TESTTIME
        for j in GEN
            @addConstraint(st, -np_cap[j] <= Pgen[t,j] - np_cap[j]*y[t,j] <= 0) # should wrap the call to y in an ancestor-type thing
        end
    end

    # need Exp here for earlier sml versions
    @setObjective(st, Min, sum{COSTSCALE*gen_cost[i]*Pgen[t,i], t=TESTTIME, i=GEN})

end


text:
using BioEnergeticFoodWebs
using Test

anyerrors = false

test_files = [
    "trophic_rank.jl",
    "checks.jl",
    "make_parameters.jl",
    "simulate.jl",
    "random.jl",
    "measures.jl",
    "rewiring/ADBM.jl",
    "rewiring/GilljamRewire.jl",
    "rewiring/StaniczenkoRewire.jl",
    "consumption.jl",
    "biological_rates.jl",
    "temperature_size.jl",
    "productivity.jl",
    "extinctions.jl",
    "mortality.jl"
]

test_n = 1
anyerrors = false

for current_test in test_files
    try
        include(current_test)
        println("\033[1m\033[32mPASSED\033[0m\t$(current_test)")
    catch e
        global anyerrors = true
        println("\033[1m\033[31mFAILED\033[0m\t$(current_test)")
        showerror(stdout, e, backtrace())
        println()
        throw("TEST FAILED")
    end
    global test_n += 1
end

if anyerrors
    throw("Test failed")
end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
# transact helper method

function transact(fa::FinancialAsset, qty::Int, fill::Float64)
    if fa.quantity + qty == 0
        fa.basis = 0
        fa.quantity = 0
    else
        fa.basis = ((fa.basis * fa.quantity) + (fill * qty)) / (fa.quantity + qty)
        fa.quantity += qty
    end
    fa
end


text:
function push_substring!(sr::StreamReassembler, data::AbstractString, index::Int, eof_flg::Bool=false)
    end_of_str = index + length(data)
    if eof_flg
        sr.eof_idx = end_of_str
        (end_of_str == sr.cur_idx) && return reach_eof(sr)
    end
    (index >= sr.capacity + sr.cur_idx || end_of_str <= sr.cur_idx || isempty(data)) && return

    start_point = max(sr.cur_idx - index, 0)
    end_point = min(sr.cur_idx + sr.capacity - index, length(data))
    abs_idx = start_point == 0 ? index : sr.cur_idx
    data = data[(1 + start_point):end_point]

    if !isempty(sr)
        prev_it = lb = searchsortedfirst(sr.unassembled_strings, abs_idx)
        if prev_it != startof(sr.unassembled_strings)
            prev_it = regress((sr.unassembled_strings, prev_it))
            k, v = deref((sr.unassembled_strings, prev_it))
            diff = k + length(v) - abs_idx
            if diff >= 0
                data = v * @view(data[1+min(end, diff):end])
                sr.unassembled_bytes -= length(v)
                abs_idx = k
            end
        end
        first_after = searchsortedafter(sr.unassembled_strings, abs_idx + length(data))
        if first_after != lb
            k, v = deref((sr.unassembled_strings, lb))
            adv = advance((sr.unassembled_strings, lb))
            while true
                k, v = deref((sr.unassembled_strings, lb))
                sr.unassembled_bytes -= length(v)
                delete!((sr.unassembled_strings, lb))
                lb = adv
                (adv == first_after) && break
                adv = advance((sr.unassembled_strings, lb))
            end
            diff = abs_idx + length(data) - k
            (diff >= 0) && (data *= @view(v[1+min(end, diff):end]))
        end
    end
    sr.unassembled_strings[abs_idx] = data
    sr.unassembled_bytes += length(data)

    if !isempty(sr)
        bg = startof(sr.unassembled_strings)
        fst = deref((sr.unassembled_strings, bg))
        (fst[1] != sr.cur_idx) && return
        s_to_write = fst[2]
        sr.unassembled_bytes -= length(s_to_write)
        sr.cur_idx += write!(sr.output, s_to_write)
        delete!((sr.unassembled_strings, bg))
        (sr.eof_idx == sr.cur_idx) && reach_eof(sr)
    end
    return
end

@inline stream_out(sr::StreamReassembler) = sr.output

unassembled_bytes(sr::StreamReassembler) = sr.unassembled_bytes

Base.isempty(sr::StreamReassembler) = isempty(sr.unassembled_strings)

function reach_eof(sr::StreamReassembler)::Nothing
    end_input!(sr.output)
    sr.unassembled_bytes = 0
    empty!(sr.unassembled_strings)
    return nothing
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
"""
    tens4checksymmetry(C4th)

If the fourth-order tensor of material elasticity has the full set of
symmetries, return true; otherwise false.
"""
function tens4checksymmetry(C4th)
	for I in 1:3
		for J in 1:3
			for K in 1:3
				for L in 1:3
					C4th[I, J, K, L] != C4th[K, L, I, J] && return false
					C4th[I, J, K, L] != C4th[J, I, K, L] && return false
					C4th[I, J, K, L] != C4th[I, J, L, K] && return false
				end
			end
		end
	end
	return true
end

"""
    tens4symmtto6x6t!(M::FMat{T}, ST::Array{T, 4}) where {T}

Convert a symmetric 4th-order tensor to a 6 x 6 matrix.

!!! Note
The order corresponds to the arrangement of the components of stress (or
strain) tensor, symmetric, three-dimensional, into a 6-component
vector.

# Example
```
J=tens4_ijkl(eye(3),eye(3))
produces the tracor:
T=rand(3);
sum(diag(T))*eye(3)
t= tens4_dot_2(J,T)
M= tens4_symm_to_6(ST)
```
"""
function tens4symmtto6x6t!(M::FMat{T}, ST::Array{T, 4}) where {T}
	# This corresponds to the arrangement of the components of stress (or
	# strain) tensor, symmetric, three-dimensional, into a 6-component
	# vector.
	ix=[1 1; 2 2; 3 3; 1 2; 1 3; 2 3];
	for  j in 1:6
		for  i in 1:6
			M[i,j] = ST[ix[i,1],ix[i,2],ix[j,1],ix[j,2]];
		end
	end
	return M
end

"""
    tens4symmt6x6tot!(ST::Array{T, 4}, M::FMat{T}) where {T}

Convert a symmetric 6 x 6 matrix to a symmetric 4th-order tensor.

!!! Note
The order corresponds to the arrangement of the components of stress (or
strain) tensor, symmetric, three-dimensional, into a 6-component
vector.
"""
function tens4symmt6x6tot!(ST::Array{T, 4}, M::FMat{T}) where {T}
	ix = [1 4 5; 4 2 6; 5 6 3]
	n = 3
	for i in 1:n
		for j in 1:n
			for k in 1:n
				for l in 1:n
					ST[i, j, k, l] = M[ix[i, j], ix[k, l]]
				end
			end
		end
	end
    return ST
end

"""
    tens4dot2!(R::Array{T, 2}, F::Array{T, 4}, S::Array{T, 2}) where {T}

Compute the double contraction of a 4th-order and a 2nd-order tensors.

!!! note
The double contraction  of two second-order sensors is defined as
`A:B = tr(A'*B) = A_ij B_ij`

The resulting second-order tensor is first zeroed out, and then the result is
accumulated.
"""
function tens4dot2!(R::Array{T, 2}, F::Array{T, 4}, S::Array{T, 2}) where {T}
	R .= zero(T)
	for  l in 1:3
		for  k in 1:3
			for  j in 1:3
				for  i in 1:3
					R[i,j] += F[i,j,k,l]*S[k,l];
				end
			end
		end
	end
	return R
end

"""
    tens4ijkl!(t::Array{T, 4}, A::FA, B::FB, op = :+) where {T, FA, FB}

Fill a 4th-order tensor as a dyadic product of two 2nd-order tensors.

The `i,j,k,l` component is given as `t[i,j,k,l]=A(i,j)*B(k,l)`.

!!! note
The tensor is accumulated to. It needs to be initialized to zero, if that is
desired as the initial state.

# Example
```
t = fill(0.0, 3, 3, 3, 3)
delta = (I, J) -> I == J ? 1.0 : 0.0
tens4ijkl!(t, delta, delta)
S = rand(3, 3)
@show tr(S) * I
tS = fill(0.0, 3, 3)
@show tens4dot2!(tS, t, S)
```
"""
function tens4ijkl!(t::Array{T, 4}, A::FA, B::FB) where {T, FA, FB}
	for  l in 1:3
		for  k in 1:3
			for  j in 1:3
				for  i in 1:3
					t[i,j,k,l] += A(i,j)*B(k,l);
				end
			end
		end
	end
	return t
end

"""
    tens4ikjl!(t::Array{T, 4}, A::FA, B::FB) where {T, FA, FB}

Fill a 4th-order tensor as a dyadic product of two 2nd-order tensors.

The `i,j,k,l` component is given as `t[i,j,k,l]=A(i,k)*B(j,l)`.

!!! note
The tensor is accumulated to. It needs to be initialized to zero, if that is
desired as the initial state.

# Example
```
t = fill(0.0, 3, 3, 3, 3)
delta = (I, J) -> I == J ? 1.0 : 0.0
tens4ikjl!(t, delta, delta)
S = rand(3, 3)
@show transpose(S)
tS = fill(0.0, 3, 3)
@show transpose(S) - tens4dot2!(tS, t, S)
```
"""
function tens4ikjl!(t::Array{T, 4}, A::FA, B::FB) where {T, FA, FB}
	for  l in 1:3
		for  k in 1:3
			for  j in 1:3
				for  i in 1:3
					t[i,j,k,l] += A(i,k)*B(j,l);
				end
			end
		end
	end
	return t
end

"""
    tens4iljk!(t::Array{T, 4}, A::FA, B::FB) where {T, FA, FB}

Fill a 4th-order tensor as a dyadic product of two 2nd-order tensors.

The `i,j,k,l` component is given as `t[i,j,k,l]=A(i,l)*B(j,k)`.

!!! note
The tensor is accumulated to. It needs to be initialized to zero, if that is
desired as the initial state.

# Example
```
t = fill(0.0, 3, 3, 3, 3)
delta = (I, J) -> I == J ? 1.0 : 0.0
tens4iljk!(t, delta, delta)
S = rand(3, 3)
tS = fill(0.0, 3, 3)
@show S - tens4dot2!(tS, t, S)
```
"""
function tens4iljk!(t::Array{T, 4}, A::FA, B::FB) where {T, FA, FB}
	for  l in 1:3
		for  k in 1:3
			for  j in 1:3
				for  i in 1:3
					t[i,j,k,l] += A(i,l)*B(j,k);
				end
			end
		end
	end
	return t
end

"""
    tens4identity!(t::Array{T, 4}) where {T}

Compute 4th-order identity tensor.

# Example

The product of the identity tensor with the second-order tensor `S` is
```
t = fill(0.0, 3, 3, 3, 3)
tens4identity!(t)
S = rand(3, 3)
tS = fill(0.0, 3, 3)
tens4dot2!(tS, t, S)
@show S - tS
```
"""
function tens4identity!(t::Array{T, 4}) where {T}
	delta = (I, J) -> I == J ? 1.0 : 0.0
	t .= zero(T)
	return tens4ikjl!(t, delta, delta)
end

"""
    tens4transposor!(t::Array{T, 4}) where {T}

Compute 4th-order transposor tensor.

# Example

The product of the transposor tensor with the second-order tensor `S` is
```
t = fill(0.0, 3, 3, 3, 3)
tens4transposor!(t)
S = rand(3, 3)
tS = fill(0.0, 3, 3)
tens4dot2!(tS, t, S)
@show S' - tS
```
"""
function tens4transposor!(t::Array{T, 4}) where {T}
	delta = (I, J) -> I == J ? 1.0 : 0.0
	t .= zero(T)
	return tens4iljk!(t, delta, delta)
end

"""
    tens4tracor!(t::Array{T, 4}) where {T}

Compute 4th-order tracor tensor.

Double contraction of a second order tensor with this fourth-order tensor
produces the spherical part of the second order tensor.

# Example

The product of the tracor tensor with the second-order tensor `S` is
```
t = fill(0.0, 3, 3, 3, 3)
tens4tracor!(t)
S = rand(3, 3)
tS = fill(0.0, 3, 3)
tens4dot2!(tS, t, S)
@show tr(S) * I - tS
```
"""
function tens4tracor!(t::Array{T, 4}) where {T}
	delta = (I, J) -> I == J ? 1.0 : 0.0
	t .= zero(T)
	return tens4ijkl!(t, delta, delta)
end

"""
    tens4symmetrizor!(t::Array{T, 4}) where {T}

Compute 4th-order symmetrizor tensor.

Double contraction of a second order tensor with this fourth-order tensor
produces the symmetric part of the second order tensor.

# Example

The product of the symmetrizor tensor with the second-order tensor `S` is
```
t = fill(0.0, 3, 3, 3, 3)
tens4symmetrizor!(t)
S = rand(3, 3)
tS = fill(0.0, 3, 3)
tens4dot2!(tS, t, S)
@show (S + S')/2 * I - tS
```
"""
function tens4symmetrizor!(t::Array{T, 4}) where {T}
	delta = (I, J) -> I == J ? 1.0 : 0.0
	t .= zero(T)
	tens4ikjl!(t, delta, delta) # identity
	tens4iljk!(t, delta, delta) # transposor
	t .*= 0.5
	return t

end

"""
    tens4skewor!(t::Array{T, 4}) where {T}

Compute 4th-order skewor tensor.

Double contraction of a second order tensor with this fourth-order tensor
produces the skew part of the second order tensor.

# Example

The product of the skewor tensor with the second-order tensor `S` is
```
t = fill(0.0, 3, 3, 3, 3)
tens4skewor!(t)
S = rand(3, 3)
tS = fill(0.0, 3, 3)
tens4dot2!(tS, t, S)
@show (S - S')/2 * I - tS
```
"""
function tens4skewor!(t::Array{T, 4}) where {T}
	delta = (I, J) -> I == J ? 1.0 : 0.0
	t .= zero(T)
	tens4iljk!(t, delta, delta) # transposor
	t .= -t # subtract that part
	tens4ikjl!(t, delta, delta) # identity
	t .*= 0.5
	return t
end

"""
    tens4deviator!(t::Array{T, 4}) where {T}

Compute 4th-order deviator tensor.

Double contraction of a second order tensor with this fourth-order tensor
produces the deviator part of the second order tensor.

# Example

The product of the deviator tensor with the second-order tensor `S` is
```
t = fill(0.0, 3, 3, 3, 3)
tens4deviator!(t)
S = rand(3, 3)
tS = fill(0.0, 3, 3)
tens4dot2!(tS, t, S)
@show tr((S - tr(S)/3*I) ), tr(tS)
```
"""
function tens4deviator!(t::Array{T, 4}) where {T}
	delta = (I, J) -> I == J ? 1.0 : 0.0
	t .= zero(T)
	tens4tracor!(t) # tracor
	t .= -(1.0/3).*t # subtract (1/3) that part
	tens4ikjl!(t, delta, delta) # identity
	t .*= 0.5
	return t
end


text:
using Test
import FlatBuffers

# generated code
include("MyGame/MyGame.jl")
import .MyGame
import .MyGame.Example
import .MyGame.Example2
import .MyGame.Example.Any_
import .MyGame.Example.Monster
import .MyGame.Example.TestSimpleTableWithEnum

# override the typeorder function, since we are working around
# circular type definitions using type parameters
FlatBuffers.typeorder(::Type{Any_}, i::Integer) = [Nothing, Monster{Any_}, TestSimpleTableWithEnum, Example2.Monster][i+1]

function loadmonsterfile(filename)
    open(joinpath(@__DIR__, filename), "r") do f Monster{Any_}(f) end
end

function checkmonster(monster)
    @test monster.hp == 80
    @test monster.mana == 150
    @test monster.name == "MyMonster"

    vec = monster.pos

    @test vec.x == 1.0
    @test vec.y == 2.0
    @test vec.z == 3.0
    @test vec.test1 == 3.0
    @test vec.test2 == MyGame.Example.ColorGreen
    @test vec.test3_a == 5
    @test vec.test3_b == 6

    monster2 = monster.test
    @test monster2.name == "Fred"

    @test length(monster.inventory) == 5
    @test sum(monster.inventory) == 10

    @test monster.vector_of_longs == [10 ^ (2*i) for i = 0:4]
    @test monster.vector_of_doubles == [-1.7976931348623157e+308, 0, 1.7976931348623157e+308]

    @test length(monster.test4) == 2

    (test0, test1) = monster.test4
    @test sum([test0.a, test0.b, test1.a, test1.b]) == 100

    @test monster.testarrayofstring == ["test1", "test2"]
    @test monster.testarrayoftables == []
    @test monster.testf == 3.14159f0
end

function checkpassthrough(monster)
    b = FlatBuffers.Builder(Monster{Any_})
    FlatBuffers.build!(b, monster)
    bytes = FlatBuffers.bytes(b)
    @test FlatBuffers.has_identifier(Monster{Any_}, bytes)
    newmonster = FlatBuffers.read(Monster{Any_}, bytes)
    checkmonster(newmonster)
end

function checkserialize(monster)
    io = IOBuffer()
    FlatBuffers.serialize(io, monster)
    bytes = take!(io)
    newmonster = FlatBuffers.deserialize(IOBuffer(bytes), Monster{Any_})
    checkmonster(newmonster)
end

@test FlatBuffers.root_type(Monster) == true
@test FlatBuffers.file_identifier(Monster) == "MONS"
@test FlatBuffers.file_extension(Monster) == "mon"

for testcase in ["test", "python_wire"]
    mon = loadmonsterfile("monsterdata_$testcase.mon")
    checkmonster(mon)
    checkpassthrough(mon)
    checkserialize(mon)
end

# test printing
mon = loadmonsterfile("monsterdata_test.mon")
b = FlatBuffers.Builder(Monster{Any_})
FlatBuffers.build!(b, mon)
io = IOBuffer()
show(io, b)
output = String(take!(io))
@test occursin("deprecated field", split(output, "\n")[9])


text:
using QuasiArrays, LinearAlgebra, Test

@testset "permutedims" begin
    A = QuasiArray(randn(5,2), (0:0.5:2, [1,4]))
    P = permutedims(A)
    @test P isa PermutedDimsQuasiArray
    @test size(P) == reverse(size(A))
    @test axes(P) == reverse(axes(A))
    @test P[4,0.5] == A[0.5,4]
    @test similar(P) isa QuasiArray
    P[4,0.5] = 2
    @test A[0.5,4] ≡ 2.0

    @testset "Diag of view" begin
        # this is used in ContinuumArrays
        D = Diagonal(view(A, 0.0, [1,4]))
        @test permutedims(D) ≡ D
    end

    @testset "vec" begin
        a = QuasiVector(randn(3),[0,2,3])
        @test permutedims(a) == a'
        @test permutedims(a .+ im) == transpose(a .+ im)
    end
end


text:
using BenchmarkTools,
    Particles,
    Distributions,
    ConjugatePriors,
    StaticArrays,
    Compat,
    Random

using ConjugatePriors: NormalInverseWishart
using Compat.LinearAlgebra

Random.seed!(100)

truth = MixtureModel([MvNormal(-2 .* ones(2), Matrix(1.0I,2,2)),
                      MvNormal(2 .* ones(2), Matrix(1.0I,2,2))])
y = [rand(truth) for _ in 1:100]

prior = NormalInverseWishart(zeros(2), 0.1, Matrix(1.0I,2,2), 3.0)

function f(y, prior)
    stateprior = ChineseRestaurantProcess(1.)
    ps = ChenLiuParticles(100, prior, stateprior)
    filter!(ps, y, false)
end

f(y, prior)
@btime(f($y, $prior))

sprior = NormalInverseWishart(SVector{2}(zeros(2)), 0.1, cholesky(SMatrix{2,2}(1.0I)), 3.0)

f(y, sprior)
@btime(f($y, $sprior))

sy = SVector{2}.(y)
f(sy, sprior)
ps3 = @btime f($sy, $sprior)


text:
using OTVisualizations
using Test

@testset "OTVisualizations.jl" begin
    # Write your tests here.
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to parse the table cells in text back-end.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    _parse_cell_text(cell::T; kwargs...)

Parse the table `cell` of type `T` and return a vector of `String` with the
parsed cell text, one component per line.
"""
function _parse_cell_text(
    cell::Any;
    autowrap::Bool = true,
    cell_data_type::DataType = Nothing,
    cell_first_line_only::Bool = false,
    column_width::Integer = -1,
    compact_printing::Bool = true,
    limit_printing::Bool = true,
    linebreaks::Bool = false,
    renderer::Union{Val{:print}, Val{:show}} = Val(:print),
    kwargs...
)
    isstring = cell_data_type <: AbstractString

    # Convert to string using the desired renderer.
    #
    # Due to the non-specialization of `data`, `cell` here is inferred as `Any`.
    # However, we know that the output of `_render_text` must be a vector of
    # String.
    cell_vstr::Vector{String} = _render_text(
        renderer, cell,
        compact_printing = compact_printing,
        isstring = isstring,
        limit_printing = limit_printing,
        linebreaks = linebreaks || cell_first_line_only
    )

    # Check if we must autowrap the text.
    autowrap && (cell_vstr = _str_autowrap(cell_vstr, column_width))

    # Check if the user only wants the first line.
    cell_first_line_only && (cell_vstr  = [first(cell_vstr)])

    return cell_vstr
end

function _parse_cell_text(
    cell::Markdown.MD;
    column_width::Integer = -1,
    linebreaks::Bool = false,
    has_color::Bool = true,
    kwargs...
)
    # The maximum size for Markdowns cells is 80.
    column_width ≤ 0 && (column_width = 80)

    # Render Markdown
    # ==========================================================================

    # First, we need to render the Markdown with all the colors.
    str = sprint(Markdown.term, cell, column_width; context = :color => true)

    # Now, we need to remove all ANSI escape sequences to count the printable
    # characters.
    #
    # This regex was obtained at:
    #
    #     https://stackoverflow.com/questions/14693701/how-can-i-remove-the-ansi-escape-sequences-from-a-string-in-python
    #
    str_nc = replace(str, r"\x1B(?:[@-Z\\-_]|\[[0-?]*[ -/]*[@-~])" => "")

    if !linebreaks
        if !has_color
            str_nc = replace(str_nc, "\n" => "\\n")
            return [str_nc]
        else
            str = replace(str, "\n" => "\\n")
            return [str]
        end
    else
        # Obtain the number of lines and the maximum number of used columns.
        tokens_nc = String.(split(str_nc, '\n'))

        if !has_color
            return tokens_nc
        else
            tokens = String.(split(str, '\n'))
            _reapply_ansi_format!(tokens)
            return tokens
        end
    end
end

function _parse_cell_text(cell::CustomTextCell; kwargs...)
    # Call the API function to reset all the fields in the custom text cell.
    reset!(cell)
    cell_vstr = parse_cell_text(cell; kwargs...)
    return cell_vstr
end

@inline _parse_cell_text(cell::Missing; kwargs...) = ["missing"]
@inline _parse_cell_text(cell::Nothing; kwargs...) = ["nothing"]
@inline _parse_cell_text(cell::UndefInitializer; kwargs...) = ["#undef"]

"""
    _process_cell_text(data::Any, i::Int, j::Int, l::Int, data_cell::Bool, data_str::String, data_len::Int, col_width::Int, crayon::Crayon, alignment::Symbol, highlighters::Ref{Any})

Process the cell by applying the right alignment and also verifying the
highlighters.
"""
function _process_cell_text(
    (@nospecialize data::Any),
    i::Int,
    j::Int,
    l::Int,
    data_cell::Bool,
    custom_cell::Bool,
    data_str::String,
    col_width::Int,
    crayon::Crayon,
    alignment::Symbol,
    (@nospecialize highlighters::Ref{Any})
)
    lstr = -1

    if data_cell
        # Check for highlighters.
        for h in highlighters.x
            if h.f(_getdata(data), i, j)
                crayon = h.fd(h, _getdata(data), i, j)::Crayon
                break
            end
        end

        # For Markdown cells, we will overwrite alignment and highlighters.
        if isassigned(data, i, j) && (data[i, j] isa Markdown.MD)
            alignment = :l
            crayon = Crayon()
            lstr = _printable_textwidth(data_str)
        end
    end

    if custom_cell
        # To align a custom text cell, we need to compute the alignment and
        # cropping data and apply it using the API functions.
        crop_chars, left_pad, right_pad = _str_compute_alignment_and_crop(
            data_str,
            alignment,
            col_width,
            -1
        )

        if crop_chars > 0
            apply_line_padding!(data[i, j], l, 0, 0)
            crop_line!(data[i, j], l, crop_chars + 1)
            append_suffix_to_line!(data[i, j], l, "…")
        else
            apply_line_padding!(data[i, j], l, left_pad, right_pad)
        end

        data_str = get_printable_cell_line(data[i, j], l)::String
    else
        # Align the string to be printed.
        data_str = _str_aligned(data_str, alignment, col_width, lstr)
    end

    return data_str, crayon
end


text:
#!/bin/env julia

# Playing with running the pricer on GPU
#
# Will port into main library once I'm happy with it

module CUDA

export mc_pricer, mc_pricer_pathdep

using CuArrays # Unfortunately this gives a warning - see https://github.com/JuliaPackaging/Requires.jl/issues/65
CuArrays.allowscalar(false)

using Statistics: mean, std

present_value(v,r,t,T) = v*CuArrays.CUDAnative.exp(-r*(T-t))

function rand_price(seed=randn();payoff=(S,K,t,T,r)->S,S_0=100,t=0,T=1,K=S_0,r=0.02,σ=0.05) # Doesn't work on GPU - doesn't like function as argument
    S = B(seed;T=T-t,B0=S_0,r=r,d=0,σ=σ)

    present_value(payoff(S,K,t,T,r),r,t,T)
end

function mc_pricer(payoff;rng=CURAND.randn,S_0=100,t=0,T=1.0,r=0.02,σ=0.05,trials=10_000)
  a = rng(Float32,trials)
  #a = B.(rng(Float32,trials),0.9) # ;B0=S_0,T=T-t,r=0.02,σ=0.05) # CuArrays doesn't like any extra args?

  # TODO: stuff gets converted to f64 here - which IIRC is much slower on GPUs
  Bbake(x) = S_0*CuArrays.CUDAnative.exp(r*T-0.5*σ^2*T+σ*CuArrays.CUDAnative.pow(T-t,0.5)*x) # sqrt(T-t) doesn't work
  prices = present_value.(payoff.(Bbake.(a)),r,t,T)
  (mean = mean(prices), sem = std(prices)/sqrt(trials))
end

# TODO:
# - this gives different results to JustJoshing.mc_pricer_pathdep - build a failing test
# - try to unify the interface with the CPU based one (tricky as we need to "reach in" and specify dims=2?)
function mc_pricer_pathdep(payoff;S_0=100,t=0:0.1:1,T=1.0f0,r=0.02f0,σ=0.05f0,trials=100_000)
  # Much faster than CPU - now check why!
  a = CuArrays.CURAND.randn(Float32,trials)

  Bbake(p,t)::Float32 = begin
    S_0*CuArrays.CUDAnative.exp(
        r*T -
        0.5f0*σ^2*T +
        p*σ*CuArrays.CUDAnative.sqrt(T-t)
    )
  end
  # If prices is not a cuarray we get a big slowdown
  # Need to ensure that our function operates on array for max
  prices::CuArray{Float32,2} = payoff(Bbake.(a,t'))

  (mean = mean(prices), sem = std(prices)/sqrt(trials))
end

end