text:
#!/usr/bin/env julia
# This script builds/pulls all static contents needed by storage server.
#
# Usage:
#  1. make sure you've added StorageMirrorServer.jl
#  2. generate/pull all tarballs: `julia gen_static_full.jl`
#  3. set a cron job to run step 2 regularly
#
# Note:
#   * Initialization would typically take days, depending on the network bandwidth and CPU
#   * set `JULIA_NUM_THREADS` to use multiple threads
#
# Disk space requirements for a complete storage (increases over time):
#   * `STATIC_DIR`: at least 500GB, would be better to have more than 3TB free space

using StorageMirrorServer
using Pkg

# This holds all the data you need to set up a storage server
# For example, my nginx service serves all files in `/mnt/mirrors` as static contents using autoindex
output_dir = ENV["LUG_path"]

# check https://status.julialang.org/ for available public storage servers
upstreams = [
    "https://kr.storage.juliahub.com",
    "https://us-east.storage.juliahub.com"
]

registries = [
    # (name, uuid, original_git_url)
    ("General", "23338594-aafe-5451-b93e-139f81909106", "https://github.com/JuliaRegistries/General")
]

# These are default parameter settings for StorageMirrorServer
# you can modify them accordingly to fit your settings
parameters = Dict(
    # if needed, you can pass custom http parameters
    :http_parameters => Dict{Symbol, Any}(
        :retry => true,
        :retries => 5,
        :readtimeout => 600
        # download data using proxy
        # it also respects `http_proxy`, `https_proxy` and `no_proxy` environment variables
        # :proxy => "http://localhost:1080"
    ),

    # whether to show the progress bar
    :show_progress => true,

    # for how long (hours) you want to skip resources in `/failed_resources.txt` until the next try
    :skip_duration => 12,
)

for reg in registries
    mirror_tarball(reg, upstreams, output_dir; parameters...)
end


text:
using Distributed: @everywhere, @spawnat, @distributed
@everywhere include("gru_dynamic_struct.jl")

const window_size = 50

create_model_definitions(layer) =
begin
    eval(Meta.parse("@everywhere " * "struct Model\n" * *(["l$i::Layer\n" for i in 1:length(layers) +1]...) * "end"))
    eval(Meta.parse("@everywhere " * "(model::Model)(io) =\n" * "begin\n" * *(["io = model.l$i(io)\n" for i in 1:length(layers) +1]...) * "end"))
end

soften = arr -> (begin
    new_arr = []
    soft = softmax(arr[1:Int(length(arr)/4)])
    for i in 1:Int(length(arr)/4)
        push!(new_arr, soft[i])
    end
    for i in Int(length(arr)/4):Int(length(arr)*3/4)
        push!(new_arr, sigm(arr[i]))
    end
new_arr
end)

prop(model, x) =
begin
    # for mfield in fieldnames(Model)
    #     layer = getfield(model, mfield)
    #     layer.state = zeros(1, length(layer.bs))
    # end
[soften(model(t)) for t in x]
end

cross_entropy(out, label) =
begin
    - sum(label .* log.(out))
end

train!(model, data, lr, batch_size) =
begin
    results = []
    for i in 1:trunc(Int, length(data)/batch_size)
        res = @distributed (vcat) for sequence in data[(i-1)*batch_size+1:i*batch_size]
            d = @diff (begin
                    sequence_parts = batchify(sequence, window_size)
                    part_losses = []
                    for part in sequence_parts
                        x = sequence[1:end-1]
                        y = sequence[2:end]
                        o = prop(model, x)
                        l = [cross_entropy(o_t, y_t) for (o_t, y_t) in zip(o, y)]
                        push!(part_losses, sum(l))
                    end
                sum(part_losses) end)
            grads = []
            for mfield in fieldnames(Model)
                layer = getfield(model, mfield)
                for lfield in fieldnames(Layer)
                    push!(grads, grad(d, getfield(layer, lfield)))
                end
            end
            @spawnat 1 print("/")
            grads, value(d)
        end
        push!(results, res)
    end

    loss = 0.0
    for (g,l) in results
        loss += l
        i = 0
        for mfield in fieldnames(Model)
            layer = getfield(model, mfield)
            for lfield in fieldnames(Layer)
                i +=1
                setfield!(layer, lfield, Param(getfield(layer, lfield) - g[i] .* lr))
            end
        end
    end
    @show loss
end

ext: jl


text:
module PDMatsExtras
using ChainRulesCore
using LinearAlgebra
using PDMats
using SparseArrays: sparse
using SuiteSparse
import Base: *, \

export PSDMat, WoodburyPDMat
export submat

include("psd_mat.jl")
include("woodbury_pd_mat.jl")
include("utils.jl")

end


text:
for (op, i) in ((:zero,:0), (:one,:1), (:two,:2), (:three,:3), (:four, :4))
  @eval begin
    function ($op){P}(::Type{ArbFloat{P}})
        z = initializer(ArbFloat{P})
        ccall(@libarb(arb_set_si), Void, (Ptr{ArbFloat}, Int), &z, $i)
        z
    end
    ($op)(::Type{ArbFloat}) = ($op)(ArbFloat{precision(ArbFloat)})
  end
end

function copy{P}(x::ArbFloat{P})
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_set), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}), &z, &x)
    z
end

function deepcopy{P}(x::ArbFloat{P})
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_set), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}), &z, &x)
    z
end

function copyradius{P}(target::ArbFloat{P}, source::ArbFloat{P})
    z = deepcopy(target)
    z.rad_exp = source.rad_exp
    z.rad_man = source.rad_man
    z
end

function deepcopyradius{P}(target::ArbFloat{P}, source::ArbFloat{P})
    target.rad_exp = source.rad_exp
    target.rad_man = source.rad_man
    target
end

function copymidpoint{P}(target::ArbFloat{P}, source::ArbFloat{P})
    z = deepcopy(target)
    z.mid_exp = source.mid_exp
    z.mid_size = source.mid_size
    z.mid_d1 = source.mid_d1
    z.mid_d2 = source.mid_d2
    z
end

"""
Rounds x to a number of bits equal to the accuracy of x (as indicated by its radius), plus a few guard bits.
The resulting ball is guaranteed to contain x, but is more economical if x has less than full accuracy.
(from arb_trim documentation)
"""
function trim{P}(x::ArbFloat{P})
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_trim), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}), &z, &x)
    z
end

"""
Rounds x to a clean estimate of x as a point value.
"""
function tidy{P}(x::ArbFloat{P})
    s = smartarbstring(x)
    ArbFloat{P}(s)
end

function decompose{P}(x::ArbFloat{P})
    # decompose x as num * 2^pow / den
    # num, pow, den = decompose(x)
    bfprec=precision(BigFloat)
    setprecision(BigFloat,P)
    bf = convert(BigFloat, x)
    n,p,d = decompose(bf)
    setprecision(BigFloat,bfprec)
    n,p,d
end

eps{P}(::Type{ArbFloat{P}}) = ldexp(1.0,-P) # for intertype workings
function eps{P}(x::ArbFloat{P})   # for intratype workings
    r = radius(x)
    if r == 0
       ldexp(1.0,-P)*x
    else
       r
    end
end

"""Similar to eps(x), epsilon(ArbFloat(x)) adjusts for the uncertainty as given by the radius.
   This function is limited to values within the range of Float64.
"""
function epsilon{P}(x::ArbFloat{P})
  if (radius(x) == 0)
     return eps(x)
  end
  midpoint_fr, midpoint_xp = frexp(eps(convert(Float64,midpoint(x))))
  radius_fr, radius_xp = frexp(eps(convert(Float64,radius(x))))
  fr = (midpoint_fr + radius_fr)*0.5
  xp = midpoint_xp + radius_xp
  if isodd(xp)
     fr = fr + 0.25
  end
  xp = xp >> 1
  ldexp(fr,xp)
end


text:
__precompile__(true)

module DistributedArrays

using Compat

if VERSION >= v"0.5.0-dev+4340"
    using Primes
    using Primes: factor
end

if VERSION < v"0.5.0-"
    typealias Future RemoteRef
    typealias RemoteChannel RemoteRef
end

importall Base
import Base.Callable
import Base.BLAS: axpy!

export (.+), (.-), (.*), (./), (.%), (.<<), (.>>), div, mod, rem, (&), (|), ($)
export DArray, SubDArray, SubOrDArray, @DArray
export dzeros, dones, dfill, drand, drandn, distribute, localpart, localindexes, ppeval, samedist
export close, darray_closeall

const registry=Dict{Tuple, Any}()
const refs=Set()  # Collection of darray identities created on this node

"""
    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...)
```
"""
type DArray{T,N,A} <: AbstractArray{T,N}
    identity::Tuple
    dims::NTuple{N,Int}
    pids::Array{Int,N}                          # pids[i]==p ⇒ processor p has piece i
    indexes::Array{NTuple{N,UnitRange{Int}},N}  # indexes held by piece i
    cuts::Vector{Vector{Int}}                   # cuts[d][i] = first index of chunk i in dimension d

    release::Bool

    function DArray(identity, dims, pids, indexes, cuts)
        # check invariants
        if dims != map(last, last(indexes))
            throw(ArgumentError("dimension of DArray (dim) and indexes do not match"))
        end
        release = (myid() == identity[1])

        global registry
        haskey(registry, (identity, :DARRAY)) && return registry[(identity, :DARRAY)]

        d = new(identity, dims, pids, indexes, cuts, release)
        if release
            push!(refs, identity)
            registry[(identity, :DARRAY)] = d

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

    DArray() = new()
end

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

## core constructors ##

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)
    identity = next_did()

    return construct_darray(identity, 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
function DArray(refs::Array{Future})
    dimdist = size(refs)
    identity = next_did()

    npids = [r.where for r in refs]
    nsizes = Array(Tuple, dimdist)
    @sync for i in 1:length(refs)
        let i=i
            @async nsizes[i] = remotecall_fetch(rr_localpart, npids[i], refs[i], identity)
        end
    end

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

    for i in 1:length(nindexes)
        subidx = ind2sub(dimdist, i)
        nindexes[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 nindexes]...)
    ncuts = Array{Int,1}[unshift!(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)]...)

    construct_darray(identity, refs, ndims, reshape(npids, dimdist), nindexes, ncuts)
end

macro DArray(ex::Expr)
    if ex.head !== :comprehension
        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}}})->($ex),
                tuple($(map(r->:(length($r)), ranges)...))) )
end

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

function construct_darray(identity, init, dims, pids, idxs, cuts)
    r=Channel(1)
    @sync begin
        for i = 1:length(pids)
            @async begin
                local typA
                if isa(init, Function)
                    typA=remotecall_fetch(construct_localparts, pids[i], init, identity, dims, pids, idxs, cuts)
                else
                    # constructing from an array of remote refs.
                    typA=remotecall_fetch(construct_localparts, pids[i], init[i], identity, dims, pids, idxs, cuts)
                end
                !isready(r) && put!(r, typA)
            end
        end
    end

    typA = take!(r)
    if myid() in pids
        d = registry[(identity, :DARRAY)]
    else
        d = DArray{eltype(typA),length(dims),typA}(identity, dims, pids, idxs, cuts)
    end
    d
end

function construct_localparts(init, identity, dims, pids, idxs, cuts)
    A = isa(init, Function) ? init(idxs[localpartindex(pids)]) : fetch(init)
    global registry
    registry[(identity, :LOCALPART)] = A
    typA = typeof(A)
    d = DArray{eltype(typA),length(dims),typA}(identity, dims, pids, idxs, cuts)
    registry[(identity, :DARRAY)] = d
    typA
end

let DID::Int = 1
    global next_did
    next_did() = (id = DID; DID += 1; (myid(), id))
end

function release_localpart(identity)
    global registry
    delete!(registry, (identity, :DARRAY))
    delete!(registry, (identity, :LOCALPART))
    nothing
end
release_localpart(d::DArray) = release_localpart(d.identity)

function close_by_identity(identity, pids)
#   @schedule println("Finalizer for : ", identity)
    global refs
    @sync begin
        for p in pids
            @async remotecall_fetch(release_localpart, p, identity)
        end
        if !(myid() in pids)
            release_localpart(identity)
        end
    end
    delete!(refs, identity)
    nothing
end

function close(d::DArray)
#    @schedule println("close : ", d.identity, ", object_id : ", object_id(d), ", myid : ", myid() )
    if (myid() == d.identity[1]) && d.release
        @schedule close_by_identity(d.identity, d.pids)
        d.release = false
    end
    nothing
end

function darray_closeall()
    global registry
    global refs
    crefs = copy(refs)
    for identity in crefs
        if identity[1] ==  myid() # sanity check
            haskey(registry, (identity, :DARRAY)) && close(registry[(identity, :DARRAY)])
            yield()
        end
    end
end

function rr_localpart(r::Future, identity)
    global registry
    lp = fetch(r)
    registry[(identity, :LOCALPART)] = lp
    return size(lp)
end

function Base.serialize(S::SerializationState, d::DArray)
    # Only send the ident for participating workers - we expect the DArray to exist in the
    # remote registry
    destpid = Base.worker_id_from_socket(S.io)
    Serializer.serialize_type(S, typeof(d))
    if (destpid in d.pids) || (destpid == d.identity[1])
        serialize(S, (true, d.identity))    # (identity_only, identity)
    else
        serialize(S, (false, d.identity))
        for n in [:dims, :pids, :indexes, :cuts]
            serialize(S, getfield(d, n))
        end
    end
end

function Base.deserialize{T<:DArray}(S::SerializationState, t::Type{T})
    what = deserialize(S)
    identity_only = what[1]
    identity = what[2]

    if identity_only
        global registry
        if haskey(registry, (identity, :DARRAY))
            return registry[(identity, :DARRAY)]
        else
            # access to fields will throw an error, at least the deserialization process will not
            # result in worker death
            d = T()
            d.identity = identity
            return d
        end
    else
        # We are not a participating worker, deser fields and instantiate locally.
        dims = deserialize(S)
        pids = deserialize(S)
        indexes = deserialize(S)
        cuts = deserialize(S)
        return T(identity, dims, pids, indexes, cuts)
    end
end

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

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

"""
    procs(d::DArray)

Get the vector of processes storing pieces of DArray `d`.
"""
Base.procs(d::DArray) = d.pids

chunktype{T,N,A}(d::DArray{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 = last(find(dims .== d))
        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 indexes for dividing sz into nc chunks
function defaultdist(sz::Int, nc::Int)
    if sz >= nc
        return round(Int, linspace(1, sz+1, nc+1))
    else
        return [[1:(sz+1);], zeros(Int, nc-sz);]
    end
end

# compute indexes 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}},chunks...)
    for cidx in CartesianRange(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.
"""
function localpart{T,N,A}(d::DArray{T,N,A})
    lpidx = localpartindex(d)
    if lpidx == 0
        return convert(A, Array(T, ntuple(zero, N)))::A
    end

    global registry
    return registry[(d.identity, :LOCALPART)]::A
end

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

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

"""
    localpart(A)

The identity when input is not distributed
"""
localpart(A) = A

"""
    localindexes(d)

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

# find which piece holds index (I...)
locate(d::DArray, I::Int...) =
    ntuple(i -> searchsortedlast(d.cuts[i], I[i]), ndims(d))

chunk{T,N,A}(d::DArray{T,N,A}, i...) = 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{T}(::Type{T}, dims::Dims, args...) = DArray(I->zeros(T,map(length,I)), dims, args...)
dzeros{T}(::Type{T}, d1::Integer, drest::Integer...) = 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{T}(::Type{T}, dims::Dims, args...) = DArray(I->ones(T,map(length,I)), dims, args...)
dones{T}(::Type{T}, d1::Integer, drest::Integer...) = 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{T}(::Type{T}, dims::Dims, args...) = DArray(I->rand(T,map(length,I)), dims, args...)
drand{T}(::Type{T}, d1::Integer, drest::Integer...) = drand(T, 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))

    owner = myid()
    rr = RemoteChannel()
    put!(rr, A)
    d = DArray(size(A), procs, dist) do I
        remotecall_fetch(() -> fetch(rr)[I...], owner)
    end
    return d
end

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

Base.convert{S,T,N}(::Type{Array{S,N}}, d::DArray{T,N}) = begin
    a = Array(S, size(d))
    @sync begin
        for i = 1:length(d.pids)
            @async a[d.indexes[i]...] = chunk(d, i)
        end
    end
    return a
end

Base.convert{S,T,N}(::Type{Array{S,N}}, s::SubDArray{T,N}) = begin
    I = s.indexes
    d = s.parent
    if isa(I,Tuple{Vararg{UnitRange{Int}}}) && S<:T && T<:S
        l = locate(d, map(first, I)...)
        if isequal(d.indexes[l...], I)
            # SubDArray corresponds to a chunk
            return chunk(d, l...)
        end
    end
    a = Array(S, size(s))
    a[[1:size(a,i) for i=1:N]...] = s
    return a
end

function Base.convert{T,N}(::Type{DArray}, SD::SubArray{T,N})
    D = SD.parent
    DArray(SD.dims, procs(D)) do I
        TR = typeof(SD.indexes[1])
        lindices = Array(TR, 0)
        for (i,r) in zip(I, SD.indexes)
            st = step(r)
            lrstart = first(r) + st*(first(i)-1)
            lrend = first(r) + st*(last(i)-1)
            if TR <: UnitRange
                push!(lindices, lrstart:lrend)
            else
                push!(lindices, lrstart:st:lrend)
            end
        end
        convert(Array, D[lindices...])
    end
end

Base.reshape{T,S<:Array}(A::DArray{T,1,S}, d::Dims) = begin
    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,sz)
        nr = size(B,1)
        sztail = size(B)[2:end]

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

            a = sub2ind(d, d1offs, globalidx...)

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

## indexing ##

getlocalindex(d::DArray, idx...) = localpart(d)[idx...]
function getindex_tuple{T}(d::DArray{T}, I::Tuple{Vararg{Int}})
    chidx = locate(d, I...)
    idxs = d.indexes[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

Base.getindex(d::DArray, i::Int) = getindex_tuple(d, ind2sub(size(d), i))
Base.getindex(d::DArray, i::Int...) = getindex_tuple(d, i)

Base.getindex(d::DArray) = d[1]
Base.getindex(d::DArray, I::Union{Int,UnitRange{Int},Colon,Vector{Int},StepRange{Int,Int}}...) = sub(d, I...)

Base.copy!(dest::SubOrDArray, src::SubOrDArray) = begin
    if !(dest.dims == src.dims &&
         procs(dest) == procs(src) &&
         dest.indexes == src.indexes &&
         dest.cuts == src.cuts)
        throw(DimensionMismatch("destination array doesn't fit to source array"))
    end
    @sync for p in procs(dest)
        @async remotecall_fetch((dest,src)->(copy!(localpart(dest), localpart(src)); nothing), p, dest, src)
    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.indexes[i]
        @async a[[I[j][K[j]] for j=1:n]...] = chunk(d, i)
    end
    return a
end

function Base.setindex!(a::Array, s::SubDArray,
        I::Union{UnitRange{Int},Colon,Vector{Int},StepRange{Int,Int}}...)
    n = length(I)
    d = s.parent
    J = s.indexes
    if length(J) < n
        a[I...] = convert(Array,s)
        return a
    end
    offs = [isa(J[i],Int) ? J[i]-1 : first(J[i])-1 for i=1:n]
    @sync for i = 1:length(d.pids)
        K_c = Any[d.indexes[i]...]
        K = [ intersect(J[j],K_c[j]) for j=1:n ]
        if !any(isempty, K)
            idxs = [ I[j][K[j]-offs[j]] for j=1:n ]
            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
                        sub(localpart(d), [K[j]-first(K_c[j])+1 for j=1:n]...)
                    end
            end
        end
    end
    return a
end

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

## higher-order functions ##

Base.map(f, d::DArray) = DArray(I->map(f, localpart(d)), d)

function Base.reduce(f, d::DArray)
    results=[]
    @sync begin
        for p in procs(d)
            @async push!(results, remotecall_fetch((f,d)->reduce(f, localpart(d)), p, f, d))
        end
    end
    reduce(f, results)
end

function _mapreduce(f, opt, d::DArray)
# TODO Change to an @async remotecall_fetch - will reduce one extra network hop -
# once bug in master is fixed.
    results=[]
    @sync begin
        for p in procs(d)
            @async push!(results, remotecall_fetch((f,opt,d)->mapreduce(f, opt, localpart(d)), p, f, opt, d))
        end
    end
    reduce(opt, results)
end
Base.mapreduce(f, opt::Union{typeof(@functorize(|)), typeof(@functorize(&))}, d::DArray) = _mapreduce(f, opt, d)
Base.mapreduce(f, opt::Function, d::DArray) = _mapreduce(f, opt, d)
Base.mapreduce(f, opt, d::DArray) = _mapreduce(f, opt, d)

Base.map!{F}(f::F, d::DArray) = begin
    @sync for p in procs(d)
        @async remotecall_fetch((f,d)->(map!(f, localpart(d)); nothing), p, f, d)
    end
    return d
end

# mapreducedim
Base.reducedim_initarray{R}(A::DArray, region, v0, ::Type{R}) = begin
    procsgrid = reshape(procs(A), size(A.indexes))
    gridsize = Base.reduced_dims(size(A.indexes), region)
    procsgrid = procsgrid[UnitRange{Int}[1:n for n = gridsize]...]
    return dfill(convert(R, v0), Base.reduced_dims(A, region), procsgrid, gridsize)
end
Base.reducedim_initarray{T}(A::DArray, region, v0::T) = Base.reducedim_initarray(A, region, v0, T)

Base.reducedim_initarray0{R}(A::DArray, region, v0, ::Type{R}) = begin
    procsgrid = reshape(procs(A), size(A.indexes))
    gridsize = Base.reduced_dims0(size(A.indexes), region)
    procsgrid = procsgrid[UnitRange{Int}[1:n for n = gridsize]...]
    return dfill(convert(R, v0), Base.reduced_dims0(A, region), procsgrid, gridsize)
end
Base.reducedim_initarray0{T}(A::DArray, region, v0::T) = Base.reducedim_initarray0(A, region, v0, T)

mapreducedim_within(f, op, A::DArray, region) = begin
    arraysize = [size(A)...]
    gridsize = [size(A.indexes)...]
    arraysize[[region...]] = gridsize[[region...]]
    return DArray(tuple(arraysize...), procs(A), tuple(gridsize...)) do I
        mapreducedim(f, op, localpart(A), region)
    end
end

function mapreducedim_between!(f, op, R::DArray, A::DArray, region)
    @sync for p in procs(R)
        @async remotecall_fetch(p, f, op, R, A, region) do f, op, R, A, region
            localind = [r for r = localindexes(A)]
            localind[[region...]] = [1:n for n = size(A)[[region...]]]
            B = convert(Array, A[localind...])
            Base.mapreducedim!(f, op, localpart(R), B)
            nothing
        end
    end
    return R
end

Base.mapreducedim!(f, op, R::DArray, A::DArray) = begin
    lsize = Base.check_reducedims(R,A)
    if isempty(A)
        return copy(R)
    end
    region = tuple(collect(1:ndims(A))[[size(R)...] .!= [size(A)...]]...)
    if isempty(region)
        return copy!(R, A)
    end
    B = mapreducedim_within(f, op, A, region)
    return mapreducedim_between!(identity, op, R, B, region)
end

Base.mapreducedim(f, op, R::DArray, A::DArray) = begin
    Base.mapreducedim!(f, op, Base.reducedim_initarray(A, region, v0), A)
end

function nnz(A::DArray)
    B = Array(Any, size(A.pids))
    @sync begin
        for i in eachindex(A.pids)
            @async B[i...] = remotecall_fetch(x -> nnz(localpart(x)), A.pids[i...], A)
        end
    end
    return reduce(+, B)
end

# LinAlg
Base.scale!(A::DArray, x::Number) = begin
    @sync for p in procs(A)
        @async remotecall_fetch((A,x)->(scale!(localpart(A), x); nothing), p, A, x)
    end
    return A
end

# TODO
# - avoid temporary darrays being created by sum, mean, std, etc when called along specific dimensions

# reduce like
for (fn, fr) in ((:sum, :+),
                 (:prod, :*),
                 (:maximum, :max),
                 (:minimum, :min),
                 (:any, :|),
                 (:all, :&))
    @eval (Base.$fn)(d::DArray) = reduce(@functorize($fr), d)
end

# mapreduce like
for (fn, fr1, fr2) in ((:maxabs, :abs, :max),
                       (:minabs, :abs, :min),
                       (:sumabs, :abs, :+),
                       (:sumabs2, :abs2, :+))
    @eval (Base.$fn)(d::DArray) = mapreduce(@functorize($fr1), @functorize($fr2), d)
end

# semi mapreduce
for (fn, fr) in ((:any, :|),
                 (:all, :&),
                 (:count, :+))
    @eval begin
        (Base.$fn)(f::typeof(@functorize(identity)), d::DArray) = mapreduce(f, @functorize($fr), d)
        (Base.$fn)(f::Base.Predicate, d::DArray) = mapreduce(f, @functorize($fr), d)
        # (Base.$fn)(f::Base.Func{1}, d::DArray) = mapreduce(f, @functorize $fr, d)
        (Base.$fn)(f::Callable, d::DArray) = mapreduce(f, @functorize($fr), d)
    end
end

# scalar ops
(+)(A::DArray{Bool}, x::Bool) = A .+ x
(+)(x::Bool, A::DArray{Bool}) = x .+ A
(-)(A::DArray{Bool}, x::Bool) = A .- x
(-)(x::Bool, A::DArray{Bool}) = x .- A
(+)(A::DArray, x::Number) = A .+ x
(+)(x::Number, A::DArray) = x .+ A
(-)(A::DArray, x::Number) = A .- x
(-)(x::Number, A::DArray) = x .- A

map_localparts(f::Callable, d::DArray) = DArray(i->f(localpart(d)), d)
map_localparts(f::Callable, d1::DArray, d2::DArray) = DArray(d1) do I
    f(localpart(d1), localpart(d2))
end
function map_localparts!(f::Callable, d::DArray)
    @sync for p in procs(d)
        @async remotecall_fetch((f,d)->(f(localpart(d)); nothing), p, f, d)
    end
    return d
end

# Here we assume all the DArrays have
# the same size and distribution
map_localparts(f::Callable, As::DArray...) = DArray(I->f(map(localpart, As)...), As[1])

for f in (:.+, :.-, :.*, :./, :.%, :.<<, :.>>, :div, :mod, :rem, :&, :|, :$)
    @eval begin
        ($f){T}(A::DArray{T}, B::Number) = map_localparts(r->($f)(r, B), A)
        ($f){T}(A::Number, B::DArray{T}) = map_localparts(r->($f)(A, r), B)
    end
end

function samedist(A::DArray, B::DArray)
    (size(A) == size(B)) || throw(DimensionMismatch())
    if (procs(A) != procs(B)) || (A.cuts != B.cuts)
        B = DArray(x->B[x...], A)
    end
    B
end

for f in (:+, :-, :div, :mod, :rem, :&, :|, :$)
    @eval begin
        function ($f){T}(A::DArray{T}, B::DArray{T})
            B = samedist(A, B)
            map_localparts($f, A, B)
        end
    end
end
for f in (:.+, :.-, :.*, :./, :.%, :.<<, :.>>)
    @eval begin
        function ($f){T}(A::DArray{T}, B::DArray{T})
            map_localparts($f, A, B)
        end
    end
end

function Base.ctranspose{T}(D::DArray{T,2})
    DArray(reverse(D.dims), procs(D)) do I
        lp = Array(T, map(length, I))
        rp = convert(Array, D[reverse(I)...])
        ctranspose!(lp, rp)
    end
end

function Base.transpose{T}(D::DArray{T,2})
    DArray(reverse(D.dims), procs(D)) do I
        lp = Array(T, map(length, I))
        rp = convert(Array, D[reverse(I)...])
        transpose!(lp, rp)
    end
end

for f in (:abs, :abs2, :acos, :acosd, :acosh, :acot, :acotd, :acoth,
          :acsc, :acscd, :acsch, :angle, :asec, :asecd, :asech, :asin,
          :asind, :asinh, :atan, :atand, :atanh, :big, :cbrt, :ceil, :cis,
          :complex, :cos, :cosc, :cosd, :cosh, :cospi, :cot, :cotd, :coth,
          :csc, :cscd, :csch, :dawson, :deg2rad, :digamma, :erf, :erfc,
          :erfcinv, :erfcx, :erfi, :erfinv, :exp, :exp10, :exp2, :expm1,
          :exponent, :float, :floor, :gamma, :imag, :invdigamma, :isfinite,
          :isinf, :isnan, :lfact, :lgamma, :log, :log10, :log1p, :log2, :rad2deg,
          :real, :sec, :secd, :sech, :sign, :sin, :sinc, :sind, :sinh, :sinpi,
          :sqrt, :tan, :tand, :tanh, :trigamma)
    @eval begin
        ($f)(A::DArray) = map($f, A)
    end
end

function mapslices{T,N}(f::Function, D::DArray{T,N}, dims::AbstractVector)
    #Ensure that the complete DArray is available on the specified dims on all processors
    for d in dims
        for idxs in D.indexes
            if length(idxs[d]) != size(D, d)
                throw(DimensionMismatch(string("dimension $d is distributed. ",
                    "mapslices requires dimension $d to be completely available on all processors.")))
            end
        end
    end

    refs = Future[remotecall((x,y,z)->mapslices(x,localpart(y),z), p, f, D, dims) for p in procs(D)]

    DArray(reshape(refs, size(procs(D))))
end

function _ppeval(f, A...; dim = map(ndims, A))
    if length(dim) != length(A)
        throw(ArgumentError("dim argument has wrong length. length(dim) = $(length(dim)) but should be $(length(A))"))
    end
    narg = length(A)
    dimlength = size(A[1], dim[1])
    for i = 2:narg
        if dim[i] > 0 && dimlength != size(A[i], dim[i])
            throw(ArgumentError("lengths of broadcast dimensions must be the same. size(A[1], $(dim[1])) = $dimlength but size(A[$i], $(dim[i])) = $(size(A[i], dim[i]))"))
        end
    end
    dims = []
    idx  = []
    args = []
    for i = 1:narg
        push!(dims, ndims(A[i]))
        push!(idx, Any[1:size(A[i], d) for d in 1:dims[i]])
        if dim[i] > 0
            idx[i][dim[i]] = 1
            push!(args, slice(A[i], idx[i]...))
        else
            push!(args, A[i])
        end
    end
    R1 = f(args...)
    ridx = Any[1:size(R1, d) for d in 1:ndims(R1)]
    push!(ridx, 1)
    Rsize = map(last, ridx)
    Rsize[end] = dimlength
    R = Array(eltype(R1), Rsize...)

    for i = 1:dimlength
        for j = 1:narg
            if dim[j] > 0
                idx[j][dim[j]] = i
                args[j] = slice(A[j], idx[j]...)
            else
                args[j] = A[j]
            end
        end
        ridx[end] = i
        R[ridx...] = f(args...)
    end

    return R
end

"""
     ppeval(f, D...; dim::NTuple)

Evaluates the callable argument `f` on slices of the elements of the `D` tuple.

#### Arguments
`f` can be any callable object that accepts sliced or broadcasted elements of `D`.
The result returned from `f` must be either an array or a scalar.

`D` has any number of elements and the alements can have any type. If an element
of `D` is a distributed array along the dimension specified by `dim`. If an
element of `D` is not distributed, the element is by default broadcasted and
applied on all evaluations of `f`.

`dim` is a tuple of integers specifying the dimension over which the elements
of `D` is slices. The length of the tuple must therefore be the same as the
number of arguments `D`. By default distributed arrays are slides along the
last dimension. If the value is less than or equal to zero the element are
broadcasted to all evaluations of `f`.

#### Result
`ppeval` returns a distributed array of dimension `p+1` where the first `p`
sizes correspond to the sizes of return values of `f`. The last dimention of
the return array from `ppeval` has the same length as the dimension over which
the input arrays are sliced.

#### Examples
```jl
addprocs(JULIA_CPU_CORES)

using DistributedArrays

A = drandn((10, 10, JULIA_CPU_CORES), workers(), [1, 1, JULIA_CPU_CORES])

ppeval(eigvals, A)

ppeval(eigvals, A, randn(10,10)) # broadcasting second argument

B = drandn((10, JULIA_CPU_CORES), workers(), [1, JULIA_CPU_CORES])

ppeval(*, A, B)
```
"""
function ppeval(f, D...; dim::NTuple = map(t -> isa(t, DArray) ? ndims(t) : 0, D))
    #Ensure that the complete DArray is available on the specified dims on all processors
    for i = 1:length(D)
        if isa(D[i], DArray)
            for idxs in D[i].indexes
                for d in setdiff(1:ndims(D[i]), dim[i])
                    if length(idxs[d]) != size(D[i], d)
                        throw(DimensionMismatch(string("dimension $d is distributed. ",
                            "ppeval requires dimension $d to be completely available on all processors.")))
                    end
                end
            end
        end
    end

    refs = Future[remotecall((x, y, z) -> _ppeval(x, map(localpart, y)...; dim = z), p, f, D, dim) for p in procs(D[1])]

    # The array of Futures has to be reshaped for the DArray constructor to work correctly.
    # This requires a fetch and the DArray is also fetching so it might be better to modify
    # the DArray constructor.
    sd = [size(D[1].pids)...]
    nd = remotecall_fetch((r)->ndims(fetch(r)), refs[1].where, refs[1])
    DArray(reshape(refs, tuple([sd[1:nd - 1], sd[end];]...)))
end

typealias DVector{T,A} DArray{T,1,A}
typealias DMatrix{T,A} DArray{T,2,A}

# Level 1

function axpy!(α, x::DVector, y::DVector)
    if length(x) != length(y)
        throw(DimensionMismatch("vectors must have same length"))
    end
    @sync for p in procs(y)
        @async remotecall_fetch(() -> (Base.axpy!(α, localpart(x), localpart(y)); nothing), p)
    end
    return y
end

function dot(x::DVector, y::DVector)
    if length(x) != length(y)
        throw(DimensionMismatch(""))
    end
    if (procs(x) != procs(y)) || (x.cuts != y.cuts)
        throw(ArgumentError("vectors don't have the same distribution. Not handled for efficiency reasons."))
    end

    results=Any[]
    @sync begin
        for i = eachindex(x.pids)
            @async push!(results, remotecall_fetch((x, y, i) -> dot(localpart(x), fetch(y, i)), x.pids[i], x, y, i))
        end
    end
    return reduce(@functorize(+), results)
end

function norm(x::DVector, p::Real = 2)
    results = []
    @sync begin
        for pp in procs(x)
            @async push!(results, remotecall_fetch(() -> norm(localpart(x), p), pp))
        end
    end
    return norm(results, p)
end

# Level 2
function add!(dest, src, scale = one(dest[1]))
    if length(dest) != length(src)
        throw(DimensionMismatch("source and destination arrays must have same number of elements"))
    end
    if scale == one(scale)
        @simd for i = eachindex(dest)
            @inbounds dest[i] += src[i]
        end
    else
        @simd for i = eachindex(dest)
            @inbounds dest[i] += scale*src[i]
        end
    end
    return dest
end

localtype{T,N,S}(A::DArray{T,N,S}) = S
localtype(A::AbstractArray) = typeof(A)

function A_mul_B!(α::Number, A::DMatrix, x::AbstractVector, β::Number, y::DVector)

    # error checks
    if size(A, 2) != length(x)
        throw(DimensionMismatch(""))
    end
    if y.cuts[1] != A.cuts[1]
        throw(ArgumentError("cuts of output vector must match cuts of first dimension of matrix"))
    end

    # Multiply on each tile of A
    R = Array(Future, size(A.pids)...)
    for j = 1:size(A.pids, 2)
        xj = x[A.cuts[2][j]:A.cuts[2][j + 1] - 1]
        for i = 1:size(A.pids, 1)
            R[i,j] = remotecall(procs(A)[i,j]) do
                localpart(A)*convert(localtype(x), xj)
            end
        end
    end

    # Scale y if necessary
    if β != one(β)
        @sync for p in y.pids
            if β != zero(β)
                @async remotecall_fetch(y -> (scale!(localpart(y), β); nothing), p, y)
            else
                @async remotecall_fetch(y -> (fill!(localpart(y), 0); nothing), p, y)
            end
        end
    end

    # Update y
    @sync for i = 1:size(R, 1)
        p = y.pids[i]
        for j = 1:size(R, 2)
            rij = R[i,j]
            @async remotecall_fetch(() -> (add!(localpart(y), fetch(rij), α); nothing), p)
        end
    end

    return y
end

function Ac_mul_B!(α::Number, A::DMatrix, x::AbstractVector, β::Number, y::DVector)

    # error checks
    if size(A, 1) != length(x)
        throw(DimensionMismatch(""))
    end
    if y.cuts[1] != A.cuts[2]
        throw(ArgumentError("cuts of output vector must match cuts of second dimension of matrix"))
    end

    # Multiply on each tile of A
    R = Array(Future, reverse(size(A.pids))...)
    for j = 1:size(A.pids, 1)
        xj = x[A.cuts[1][j]:A.cuts[1][j + 1] - 1]
        for i = 1:size(A.pids, 2)
            R[i,j] = remotecall(() -> localpart(A)'*convert(localtype(x), xj), procs(A)[j,i])
        end
    end

    # Scale y if necessary
    if β != one(β)
        @sync for p in y.pids
            if β != zero(β)
                @async remotecall_fetch(() -> (scale!(localpart(y), β); nothing), p)
            else
                @async remotecall_fetch(() -> (fill!(localpart(y), 0); nothing), p)
            end
        end
    end

    # Update y
    @sync for i = 1:size(R, 1)
        p = y.pids[i]
        for j = 1:size(R, 2)
            rij = R[i,j]
            @async remotecall_fetch(() -> (add!(localpart(y), fetch(rij), α); nothing), p)
        end
    end
    return y
end

# Level 3
function _matmatmul!(α::Number, A::DMatrix, B::AbstractMatrix, β::Number, C::DMatrix, tA)
    # error checks
    Ad1, Ad2 = (tA == 'N') ? (1,2) : (2,1)
    mA, nA = size(A, Ad1, Ad2)
    mB, nB = size(B)
    if mB != nA
        throw(DimensionMismatch("matrix A has dimensions ($mA, $nA), matrix B has dimensions ($mB, $nB)"))
    end
    if size(C,1) != mA || size(C,2) != nB
        throw(DimensionMismatch("result C has dimensions $(size(C)), needs ($mA, $nB)"))
    end
    if C.cuts[1] != A.cuts[Ad1]
        throw(ArgumentError("cuts of the first dimension of the output matrix must match cuts of dimension $Ad1 of the first input matrix"))
    end

    # Multiply on each tile of A
    if tA == 'N'
        R = Array(Future, size(procs(A))..., size(procs(C), 2))
    else
        R = Array(Future, reverse(size(procs(A)))..., size(procs(C), 2))
    end
    for j = 1:size(A.pids, Ad2)
        for k = 1:size(C.pids, 2)
            Acuts = A.cuts[Ad2]
            Ccuts = C.cuts[2]
            Bjk = B[Acuts[j]:Acuts[j + 1] - 1, Ccuts[k]:Ccuts[k + 1] - 1]
            for i = 1:size(A.pids, Ad1)
                p = (tA == 'N') ? procs(A)[i,j] : procs(A)[j,i]
                R[i,j,k] = remotecall(p) do
                    if tA == 'T'
                        return localpart(A).'*convert(localtype(B), Bjk)
                    elseif tA == 'C'
                        return localpart(A)'*convert(localtype(B), Bjk)
                    else
                        return localpart(A)*convert(localtype(B), Bjk)
                    end
                end
            end
        end
    end

    # Scale C if necessary
    if β != one(β)
        @sync for p in C.pids
            if β != zero(β)
                @async remotecall_fetch(() -> (scale!(localpart(C), β); nothing), p)
            else
                @async remotecall_fetch(() -> (fill!(localpart(C), 0); nothing), p)
            end
        end
    end

    # Update C
    @sync for i = 1:size(R, 1)
        for k = 1:size(C.pids, 2)
            p = C.pids[i,k]
            for j = 1:size(R, 2)
                rijk = R[i,j,k]
                @async remotecall_fetch(d -> (add!(localpart(d), fetch(rijk), α); nothing), p, C)
            end
        end
    end
    return C
end

A_mul_B!(α::Number, A::DMatrix, B::AbstractMatrix, β::Number, C::DMatrix) = _matmatmul!(α, A, B, β, C, 'N')
Ac_mul_B!(α::Number, A::DMatrix, B::AbstractMatrix, β::Number, C::DMatrix) = _matmatmul!(α, A, B, β, C, 'C')
At_mul_B!(α::Number, A::DMatrix, B::AbstractMatrix, β::Number, C::DMatrix) = _matmatmul!(α, A, B, β, C, 'T')
At_mul_B!(C::DMatrix, A::DMatrix, B::AbstractMatrix) = At_mul_B!(one(eltype(C)), A, B, zero(eltype(C)), C)

function (*)(A::DMatrix, x::AbstractVector)
    T = promote_type(Base.LinAlg.arithtype(eltype(A)), Base.LinAlg.arithtype(eltype(x)))
    y = DArray(I -> Array(T, map(length, I)), (size(A, 1),), procs(A)[:,1], (size(procs(A), 1),))
    return A_mul_B!(one(T), A, x, zero(T), y)
end
function (*)(A::DMatrix, B::AbstractMatrix)
    T = promote_type(Base.LinAlg.arithtype(eltype(A)), Base.LinAlg.arithtype(eltype(B)))
    C = DArray(I -> Array(T, map(length, I)), (size(A, 1), size(B, 2)), procs(A)[:,1:min(size(procs(A), 2), size(procs(B), 2))], (size(procs(A), 1), min(size(procs(A), 2), size(procs(B), 2))))
    return A_mul_B!(one(T), A, B, zero(T), C)
end

function Ac_mul_B(A::DMatrix, x::AbstractVector)
    T = promote_type(Base.LinAlg.arithtype(eltype(A)), Base.LinAlg.arithtype(eltype(x)))
    y = DArray(I -> Array(T, map(length, I)), (size(A, 2),), procs(A)[1,:], (size(procs(A), 2),))
    return Ac_mul_B!(one(T), A, x, zero(T), y)
end
function Ac_mul_B(A::DMatrix, B::AbstractMatrix)
    T = promote_type(Base.LinAlg.arithtype(eltype(A)), Base.LinAlg.arithtype(eltype(B)))
    C = DArray(I -> Array(T, map(length, I)), (size(A, 2), size(B, 2)), procs(A)[1:min(size(procs(A), 1), size(procs(B), 2)),:], (size(procs(A), 2), min(size(procs(A), 1), size(procs(B), 2))))
    return Ac_mul_B!(one(T), A, B, zero(T), C)
end

end # module


text:
using COSMO, Random, Test, Pkg
rng = Random.MersenneTwister(12345)

include("./UnitTests/COSMOTestUtils.jl")

@testset "COSMO native tests" begin

  include("./UnitTests/simple.jl")
  include("./UnitTests/sets.jl")
  include("./UnitTests/constraints.jl")
  include("./UnitTests/kktsolver.jl")
  include("./UnitTests/model.jl")
  include("./UnitTests/qp-box.jl")
  include("./UnitTests/socp-lasso.jl")
  include("./UnitTests/closestcorr.jl")
  include("./UnitTests/exp_cone.jl")
  include("./UnitTests/pow_cone.jl")
  include("./UnitTests/print.jl")
  include("./UnitTests/InfeasibilityTests/runTests.jl")
  include("./UnitTests/algebra.jl")
  include("./UnitTests/splitvector.jl")
  include("./UnitTests/interface.jl")
  include("./UnitTests/chordal_decomposition_triangle.jl")
  include("./UnitTests/psd_completion.jl")

    # optional unittests
  if in("Pardiso",keys(Pkg.installed()))
    include("./UnitTests/options_factory.jl")
  end

end


text:
type Results
    ptr::Ref{Ptr{Void}}
    ptr2::Ref{Ptr{Void}}
    dim0::Ref{DimT}
    dim1::Ref{DimT}
    dim2::Ref{DimT}
    dim3::Ref{DimT}
    dims::Vector{DimT}
    pointers::Vector{Ptr{Void}}
    dType::Ref{DType}
    bool::Ref{Bool}
    double::Ref{Float64}
    double2::Ref{Float64}
    uint::Ref{UInt32}
end

Results() = Results(
    Ref{Ptr{Void}}(C_NULL),
    Ref{Ptr{Void}}(C_NULL),
    Ref{DimT}(0),
    Ref{DimT}(0),
    Ref{DimT}(0),
    Ref{DimT}(0),
    Vector{DimT}(4),
    Vector{Ptr{Void}}(4),
    Ref{DType}(0),
    Ref{Bool}(false),
    Ref{Float64}(0.0),
    Ref{Float64}(0.0),
    Ref{UInt32}(0))

function collectDims(dims::Vector{DimT}, values)
    idx = 1
    for v in values
        dims[idx] = v
        idx += 1
    end
    for idx2 in idx:4
        dims[idx2] = 1
    end
    dims
end

function collectArrayPointers(ptrs::Vector{Ptr{Void}}, arrays)
    idx = 1
    for arr in arrays
        verifyAccess(arr)
        ptrs[idx] = arr.ptr
        idx += 1
    end
    return idx - 1, ptrs
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
using SWIP
using Base.Test

# Untemplated model
H ~ :H normal(0, 1)
X = normal(θ, 1)
Y = normal(X, Y)
X ⟂ (Y ¦ x => 2)

t = trace((H = 1.23, X = 1.34. Y = 1.23))

# Templated model
x ~ Real

# Unspecified Distributions
X ⟂ (Y ¦ x => x)
H ~ ud(0, 1)
X = ud(θ, 1)
Y = ud(X, Y)
@test_throws UninterpretedDensity sample(Y)


text:
function besselj_asymptotic_zero(nu::T, k::Integer) where {T <: AbstractFloat}
    T(π) * (k-1 + nu/2 + 3*one(T)/4)
end

function besselj_zero(nu, k)
    z = besselj_asymptotic_zero(nu, k)
    bf = (x) -> besselj(nu, x)
    return fzero(bf,z)
end

function besselroots(nu::AbstractFloat, n)
    [besselj_zero(nu, k) for k in 1:n]
end

asy_gaussjacobi(n::Integer, α::Number, β::Number) =
    asy_gaussjacobi(n, promote(float(α), float(β))...)

function asy_gaussjacobi(n::Integer, α::T, β::T) where {T <: AbstractFloat}
    x = zeros(T, n)
    w = zeros(T, n)
    asy_gaussjacobi!(x, w, α, β)
end

# The expansions are in terms of z = 1/(2n+α+β+1)
jac_z(n, α, β) = 1/(2n+α+β+1)

function asy_gaussjacobi!(x, w, α, β)
    @assert length(x) == length(w)

    n = length(x)
    k_left = ceil(Int, sqrt(n))
    k_right = n-k_left+1

    z = jac_z(n, α, β)
    jb = besselroots(β, k_left)
    δ = zero(eltype(x))
    for k in 1:k_left
        x_left, w_left, δ_left = asy_jacobi_leftendpoint(n, k, α, β, jb[k], z)
        x_bulk, w_bulk, δ_bulk = asy_jacobi_bulk(n, k, α, β, z)
        if δ_left < δ_bulk
            x[k], w[k] = x_left, w_left
            δ = max(δ, δ_left)
        else
            x[k], w[k] = x_bulk, w_bulk
            δ = max(δ, δ_bulk)
        end
    end
    ja = besselroots(α, k_left)
    for k in k_right:n
        x_right, w_right, δ_right = asy_jacobi_rightendpoint(n, k, α, β, ja[n-k+1], z)
        x_bulk, w_bulk, δ_bulk = asy_jacobi_bulk(n, k, α, β, z)
        if δ_right < δ_bulk
            x[k], w[k] = x_right, w_right
            δ = max(δ, δ_right)
        else
            x[k], w[k] = x_bulk, w_bulk
            δ = max(δ, δ_bulk)
        end
    end
    for k in k_left+1:k_right-1
        x[k], w[k], δ_bulk = asy_jacobi_bulk(n, k, α, β, z)
        δ = max(δ, δ_bulk)
    end
    x, w, δ
end

two(T) = 2*one(T)

function asy_jacobi_leftendpoint(n::Int, k::Int, α, β, jbk, z = jac_z(n, α, β))
    T = promote_type(typeof(α), typeof(β))
    x = zero(T)
    w = zero(T)

    # Add terms one by one and continue until the next term is larger than the previous one
    x3 = 2*jbk^2*z^2
    δ = abs(x3)
    if δ > eps(T) && δ < 1
        x += x3

        x5 = - 2/3*(jbk^4 - (3α^2 + β^2 - 1)*jbk^2)*z^4
        w5 = 1/3*(3α^2 + β^2 - 2*jbk^2 - 1)*z^2
        δ = abs(x5)
        if δ > eps(T) && δ < abs(x3)
            x += x5
            w += w5

            x7 = two(T)/45*(2*jbk^6 - 3*(5α^2 + 3β^2 - 2)*jbk^4 + (45α^4 + 7β^4 + 20*(3α^2 - 1)*β^2 - 60α^2 + 13)*jbk^2)*z^6
            w7 = one(T)/45*(45α^4 + 7β^4 + 6*jbk^4 + 20*(3α^2- 1)*β^2
                - 6*(5α^2 + 3β^2 - 2)*jbk^2 - 60α^2 + 13)*z^4
            δ = abs(x7)
            if δ > eps(T) && δ < abs(x5)
                x += x7
                w += w7

                x9 = - two(T)/2835*(9*jbk^8 - 18*(7α^2 + 5β^2 - 3)*jbk^6
                        + (328β^4 + (1512α^2 - 575)*β^2 + 567α^2 -113)*jbk^4
                        - (2835α^6 + 247β^6 + 1407*(3α^2 - 1)*β^4 - 8505α^4 + 21*(405α^4 - 600α^2 + 133)*β^2 +8379α^2 - 1633)*jbk^2)*z^8
                w9 = one(T)/2835*(2835α^6 + 247β^6 - 36*jbk^6 + 1407*(3α^2 - 1)*β^4
                        + 54*(7α^2 + 5β^2 - 3)*jbk^4 - 8505α^4
                        + 21*(405α^4 - 600α^2 + 133)*β^2 - 2*(328β^4 +
                        (1512α^2 - 575)*β^2 + 567α^2 - 113)*jbk^2 + 8379α^2 -1633)*z^6
                δ = abs(x9)
                if δ > eps(T) && δ < abs(x7)
                    x += x9
                    w += w9

                    x11 = two(T)/42525*(6*jbk^10 - 15*(9α^2 + 7β^2 - 4)*jbk^8 + (6615α^4
                            + 769β^4 + 2*(1620α^2 - 589)*β^2 - 12150α^2 + 2668)*jbk^6 +
                            3*(9450α^6 - 999β^6 - 23*(400α^2 - 147)*β^4 - 40635α^4 -
                            (2835α^4 + 1850α^2 - 294)*β^2 + 50650α^2 - 10236)*jbk^4 +
                            (42525α^8 + 2327β^8 + 22340*(3α^2 - 1)*β^6 - 226800α^6 + 168*(1530α^4 -
                            2415α^2 + 542)*β^4 + 517860α^4 + 20*(11340α^6 - 38745α^4 + 42399α^2
                             - 8488)*β^2 - 509280α^2 + 98717)*jbk^2)*z^10
                     if abs(x11) > eps(T) && abs(x11) < abs(x9)
                         x += x11
                     end
                end
            end
        end
    end

    x += -1
    jbm = besselj(β-1, jbk)
    w = 8*z^2/jbm^2 * (1-x)^α*(1+x)^β * (1+w)
    x, w, δ
end

function asy_jacobi_rightendpoint(n, k, α, β, jak, z = jac_z(n, α, β))
    x, w, δ = asy_jacobi_leftendpoint(n, n-k, β, α, jak, z)
    -x, w, δ
end

function asy_jacobi_bulk(n, k, α, β, z = jac_z(n, α, β))

    t = cos(pi * (4n-4k+2α+3)/(4n+2α+2β+2))

    T = eltype(α)
    x = zero(T)
    w = zero(T)

    x3 = one(T)/2*(2α^2 - 2β^2 + (2α^2 + 2β^2 - 1)*t)*z^2
    w3 = -(1 -2α^2 -2β^2)*z^2

    δ = abs(x3)
    if δ > eps(T) && δ < abs(t)
        x += x3
        w += w3

        x5 = - one(T)/24*(32α^4 - 32β^4 - (16α^4 + 16β^4 + 4*(12α^2 -5)*β^2 - 20α^2 + 5)*t^3
            - 24*(α^2 - β^2)*t^2 - 40α^2 + 40β^2+ 3*(16α^4 + 16β^4 + 4*(4α^2 - 7)*β^2
            - 28α^2 + 11)*t)*z^4/(t^2- 1)
        w5 = one(T)/12*((16α^4 + 16β^4 + 4*(12α^2 - 5)*β^2 - 20α^2 + 5)*t^4
            + 48α^4 + 48β^4 + 12*(4α^2 - 7)*β^2 - 6*(4*(4α^2 + 1)*β^2
            + 4α^2 - 3)*t^2 - 84α^2 + 64*(α^4 - β^4 - 2α^2 + 2β^2)*t +33)/(t^4 - 2*t^2 + 1)*z^4
        δ = abs(x5)
        if δ > eps(T) && δ < abs(x3)
            x += x5
            w += w5

            x7 = one(T)/240*(576α^6 - 576β^6 + (96α^6 + 96β^6 + 80*(8α^2 - 3)*β^4
                - 240α^4 + 2*(320α^4 - 440α^2 + 101)*β^2 + 202α^2 - 39)*t^5
                -320*(α^2 - 6)*β^4 + 240*(α^2 - β^2)*t^4 - 1920α^4
                -10*(32*(5α^2 + 3)*β^4 + 96α^4 + 2*(80α^4 - 152α^2 - 97)*β^2
                - 194α^2 + 99)*t^3 + 16*(20α^4 - 127)*β^2 + 160*(6α^6 - 6β^6 +
                2*(α^2 + 15)*β^4 - 30α^4 - (2α^4 + 41)*β^2 + 41α^2)*t^2 +
                2032α^2 + 15*(96α^6 + 96β^6 + 16*(4α^2 - 23)*β^4 - 368α^4
                + 2*(32α^4 - 72α^2 + 223)*β^2 + 446α^2 - 173)*t)*z^6/(t^4 - 2*t^2 +1)
            w7 = one(T)/120*((96α^6 + 96β^6 + 80*(8α^2 - 3)*β^4 - 240α^4 + 2*(320α^4 - 440α^2 +101)*β^2
                + 202α^2 - 39)*t^6 - 1440α^6 - 1440β^6 - 240*(4α^2 - 23)*β^4 - 5*(96α^6+ 96β^6
                + 16*(20α^2 - 27)*β^4 - 432α^4 + 2*(160α^4 - 136α^2 + 295)*β^2 + 590α^2 - 237)*t^4
                + 5520α^4 - 640*(3α^6 - 3β^6 + (α^2 + 15)*β^4 - 15α^4
                -(α^4 + 22)*β^2 + 22α^2)*t^3 - 30*(32α^4 - 72α^2 + 223)*β^2 - 15*(288α^6 +288β^6
                - 16*(8α^2 + 81)*β^4 - 1296α^4 - 2*(64α^4 - 88α^2 - 863)*β^2 +1726α^2 - 717)*t^2
                - 6690α^2 - 128*(33α^6 - 33β^6 - 5*(α^2 - 27)*β^4 - 135α^4 +(5α^4 - 166)*β^2 + 166α^2)*t
                + 2595)*z^6/(t^6 - 3*t^4 + 3*t^2 - 1)
            δ = abs(x7)
            if δ > eps(T) && δ < abs(x5)
                x += x7
                w += w7

                x9 = -one(T)/40320*(219648α^8 - 219648β^8 - (9728α^8 + 9728β^8 + 896*(138α^2 - 49)*β^6
                    - 43904α^6 + 224*(1160α^4 - 1720α^2 + 389)*β^4 + 87136α^4
                    + 8*(15456α^6 -48160α^4 + 49364α^2 - 9785)*β^2 - 78280α^2 + 14921)*t^7 - 10752*(14α^2 - 127)*β^6
                    -40320*(α^2 - β^2)*t^6 - 1365504α^6 + 21*(2048α^8 + 2048β^8 + 128*(146α^2 -123)*β^6
                    - 15744α^6 + 32*(1320α^4 - 1760α^2 + 2023)*β^4 + 64736α^4 + 8*(2336α^6- 7040α^4
                    + 3644α^2 - 11275)*β^2 - 90200α^2 + 37111)*t^5 + 75264*(5α^2 - 49)*β^4
                    + 4480*(44α^8 - 44β^8 + 8*(3α^2 + 55)*β^6 - 440α^6 - 24*(7α^2 + 72)*β^4
                    +1728α^4 - (24α^6 - 168α^4 - 2405)*β^2 - 2405α^2)*t^4 + 3687936α^4
                    + 105*(6656α^8 + 6656β^8 - 128*(50α^2 + 443)*β^6 - 56704α^6 - 32*(296α^4 -696α^2 - 6027)*β^4
                    + 192864α^4 - 8*(800α^6 - 2784α^4 + 3580α^2 + 30285)*β^2 -242280α^2 + 99933)*t^3
                    + 384*(392α^6 - 980α^4 + 10527)*β^2 + 2688*(424α^8 - 424β^8+ 4*(4α^2 + 783)*β^6 - 3132α^6
                    + 4*(35α^2 - 2407)*β^4 + 9628α^4 - (16α^6 +140α^4 - 11429)*β^2 - 11429α^2)*t^2
                    - 4042368α^2 + 35*(23552α^8 + 23552β^8 +128*(90α^2 - 1231)*β^6 - 157568α^6
                    + 32*(328α^4 - 1376α^2 + 14095)*β^4 +451040α^4 + 8*(1440α^6 - 5504α^4 + 9964α^2 - 65439)*β^2
                    - 523512α^2 +206379)*t)*z^8/(t^2 - 1)^3
                if abs(x9) > eps(T) && abs(x9) < abs(x7)
                    x += x9
                end
            end
        end
    end

    x += t
    w = pi * sqrt(1-t^2) * z * (1-x)^α*(1+x)^β * (2+w)
    x, w, δ
end


text:
# using Random, Statistics, DelimitedFiles, Arpack, Distributions
# using Images, CSV, LinearAlgebra, SparseArrays, Plots

using Plots, Images, DelimitedFiles, Statistics

begin
  function init2(arg)            #THIS FUNCTION IMPORTS THE SKELETONIZED PIC AND
                              #RETURNS THE GRAYSCALE PIC
    if arg==0
      estrin="/Users/javier/Desktop/Retinopathy/DRIVE/test/1st_manual/01_manual1_Sket.tif"
    else
      estrin=string(arg)
    end
    pic=load(estrin)
    GPic=Gray.(pic)
    return GPic
  end

  function addZeros(a)
    A=channelview(a)
    if size(A)[1]- size(A)[2]>0
      A=hcat(A,fill(0.0, (size(A)[1], size(A)[1]- size(A)[2]) ))
    elseif size(A)[1]- size(A)[2]<0
      A=vcat(A,fill(0.0, ( -size(A)[1]+ size(A)[2], size(A)[2]) ))
    end
    return A
  end

  function checkIfNeighborIsNode(nod,x,y)
    arj= [x y].==nod
    r = arj[:,1]'*arj[:,2]

    if r==1
      return true
    elseif r==0
      return false
    else
      @warn "Your code is doing some nasty things... Burn it and throw it in a river!!"
      return true
    end
  end
  println(pwd())
  function removeOverlapNodes(a)
    newA=copy(a)              #SKELETONIZED WITHOUT NODES (THIS MATRIX IS BINARY)
    nod=[0 0]               #POSITION OF NODES AND ENDPOINTS
    numberOfSegPerNode=[0]    #NUMBER OF LINKS FROM GIVEN NODE (COUPLED WITH PREVIOUS ARRAY)
    adjMat=zeros(size(a)[1],size(a)[2])       #MISLEADING NAME. BINARY MATRIX WHICH CONTAINS ONLY THE NODES
    endpoints=zeros(size(a)[1],size(a)[2])    #BINARY MATRIX WHICH CONTAINS ONLY THE ENDPOINTS
    ep=[0 0]          #POSITIONS OF ENDPOINTS

    for i=1:size(a)[1], j=1:size(a)[2]
      counter=0;
      isitaNeighbor=false
      if a[i,j]==1.0
        if a[i-1,j]==1.0
          counter +=1
          isitaNeighbor=checkIfNeighborIsNode(nod,i-1,j)
        end
        if a[i+1,j]==1.0
          counter +=1
          isitaNeighbor = isitaNeighbor ? true : checkIfNeighborIsNode(nod,i+1,j)
        end
        if a[i,j-1]==1.0
          counter +=1
          isitaNeighbor = isitaNeighbor ? true : checkIfNeighborIsNode(nod,i,j-1)
        end
        if a[i,j+1]==1.0
          counter +=1
          isitaNeighbor = isitaNeighbor ? true : checkIfNeighborIsNode(nod,i,j+1)
        end
        if a[i+1,j+1]==1.0
          counter +=1
          isitaNeighbor = isitaNeighbor ? true : checkIfNeighborIsNode(nod,i+1,j+1)
        end
        if a[i-1,j+1]==1.0
          counter +=1
          isitaNeighbor = isitaNeighbor ? true : checkIfNeighborIsNode(nod,i-1,j+1)
        end
        if a[i+1,j-1]==1.0
          counter +=1
          isitaNeighbor = isitaNeighbor ? true : checkIfNeighborIsNode(nod,i+1,j-1)
        end
        if a[i-1,j-1]==1.0
          counter +=1
          isitaNeighbor = isitaNeighbor ? true : checkIfNeighborIsNode(nod,i-1,j-1)
        end
      end

      if isitaNeighbor
        #println("It is a neighbor node!")
      else
        if counter>2 || counter==1
          #println("removing pixel")
          for k=-1:1, l=-1:1
            newA[i+k,j+l]=0.0
          end
          adjMat[i,j]=1.0
          nod=vcat(nod,[i j])
          append!(numberOfSegPerNode,counter)
          if counter==1
            #adjMat[i,j]=1.0
            endpoints[i,j]=1.0
            newA[i,j]=0.0
            ep=vcat(ep,[i j])
            #println("1 neighbor")
          end
        end
      end
    end
    return newA, nod[2:size(nod)[1],:], adjMat, ep[2:size(ep,1),:], numberOfSegPerNode[2:size(numberOfSegPerNode)[1]]
  end

  function algSearch(a, nod, index, newNode, ASeg; tMAX=3000)
    initialX, initialY=nod[newNode,1], nod[newNode,2]
    oldvalues=[initialX initialY]
    initialnod=[nod[newNode,1] nod[newNode,2]]
    firstneighborlist=[0 0]
    outputNodes=[0 0]
    diamMean=[0.0]
    traj=[initialX initialY]
    qq=false
    for ind=1:index
      diamArray=[0]

      for tim=1:tMAX
          if tim>1
              newx, newy, qq = secondNeiIsNode(nod, initialX, initialY, oldvalues[1], oldvalues[2],initialnod)
              if qq
                initialX=newx
                initialY=newy
                traj=vcat(traj,[initialX initialY])
                #println("if qq " * string([initialX initialY]) * "  " * string([nod[newNode,1] nod[newNode,2]]))
                break
              end
          end

          for k=-1:1, l=-1:1
            if [k,l] !=[0,0] && a[initialX+k, initialY+l]==1 && [initialX+k initialY+l]!=oldvalues && [initialX+k initialY+l]!=initialnod && NotInthelis([initialX+k initialY+l], firstneighborlist)
              #if ind==index
                oldvalues=[initialX initialY]
                initialX=initialX+k
                initialY=initialY+l
                traj=vcat(traj,[initialX initialY])

                if tim<=100 #tim>1 && tim%2==0 &&
                  diam=measureDiam(ASeg, oldvalues[1], oldvalues[2], initialX, initialY)
                  if diam <=10
                    #global diamMat[oldvalues[1]-Int(floor(diam/2)*(oldvalues[2]-initialY)): oldvalues[1]+ Int(floor(diam/2)*(oldvalues[2]-initialY)), oldvalues[2]- Int(floor(diam/2)*(initialX-oldvalues[1])): oldvalues[2]+ Int(floor(diam/2)*(initialX-oldvalues[1])) ].=1
                    append!(diamArray,diam)
                  end
                end

                #println("if qq NOT " * string([initialX initialY]) * " " * string(firstneighborlist))
                if tim==1
                  firstneighborlist=vcat(firstneighborlist,[initialX initialY])
                end
                if tim>1 && [initialX,initialY] !=initialnod && abs(initialX-initialnod[1])<2 && abs(initialY-initialnod[2])<2
                  firstneighborlist=firstneighborlist[1:size(firstneighborlist,1)-1,:]
                  firstneighborlist=vcat(firstneighborlist,[initialX initialY])
                end
                break
            end
          end
          if tim==tMAX
           @warn "Reach tMax= $(tMAX)!! $(ind) at node $(newNode)"
           initialX=initialnod[1]
           initialY=initialnod[2]
           traj=vcat(traj,[0 0])
          end
      end
      tmpList=[initialX initialY] .== outputNodes

      if tmpList[:,1]'*tmpList[:,2]==0 && [initialX initialY]!=initialnod
        outputNodes=vcat(outputNodes,[initialX initialY])
      end
      initialX, initialY=nod[newNode,1], nod[newNode,2]
      oldvalues=[initialX initialY]
      traj=vcat(traj, [initialX initialY])

      popfirst!(diamArray)
      if length(diamArray)==0
        length(diamMean)>2  ? append!(diamMean,diamMean[size(diamMean,1)]) : append!(diamMean,1)
      else
        append!(diamMean,mean(diamArray))
      end
      #println(size(diamMean))
      length(diamArray) < 100 ? append!(diamArray,zeros(100-length(diamArray))) : 1

      #append!(diamStat,diamArray)

    end
    return outputNodes[2:size(outputNodes,1),:], diamMean[2:size(diamMean,1)], traj
  end

  function NotInthelis(l, firstneighborlist)
    activate=true;
    for i=1:size(firstneighborlist)[1]
      if l ==firstneighborlist[i,:]'
        activate=false
        break
      end
    end
    return activate
  end

  function secondNeiIsNode(nod, x, y, oldx, oldy, initnode)
    r=0
    newx=0
    newy=0
    for j=-1:1, k=-1:1
      if [x+j,y+k]!=[oldx,oldy] && [x+j y+k]!= initnode
        arj= [x+j y+k].==nod
        r = arj[:,1]'*arj[:,2]

        if r==1
          newx=x+j
          newy=y+k
          break
        end
      end
    end

    if r==1
      return newx, newy, true
    else
      return x,y,false
    end
  end

  function NotANode(l,nod, firsttime, initNod)     #IF ACTIVATE IS TRUE, RETURN FALSE.
                                          # BASICALLY, IF l IS A NODE, THEN STOP THERE.
    activate=false;
    for i=1:size(nod)[1]
      if l==nod[i,:] && l!=initNod
        activate=true;
      end
    end

    if firsttime==1
      return true
    else
      if activate
        return false
      else
        return true
      end
    end
  end

  function NotANode2(nod, cNode, index)
    activate=false
    for i=1:size(cNode)[1]
      if nod[index]==cNode[i,:]
        activate=true;
      end
    end

      if activate
        return false
      else
        return true
      end
  end

  function getIndex(nod, x, y)
    aa=nod .==[x y]
    b=aa[:,1].*aa[:,2]
    ind =findall(x-> x==true, b)
    return ind[1]
  end

  function IterationOnAlgSearch(A, nod, numberOfSegPerNode, ASeg)
    #array1=zeros(size(A)[1],size(A)[2])
    #array1[nod[1,1],nod[1,2]]=1.0
    capturedNodes=nod[1,:]'
    adjMat=zeros(size(nod)[1],size(nod)[1])
    adjMatD=zeros(size(nod)[1],size(nod)[1])
    outputNodes, diamMean, traj=algSearch(A,nod,numberOfSegPerNode[1],1, ASeg)
    for j=1:size(outputNodes,1)
      vere, indie = calDist(nod, outputNodes, 1, j, traj)
      if vere
        outputNodes[j,:]=traj[indie,:]
        nod=vcat(nod,outputNodes[j,:]')
        append!(numberOfSegPerNode,2)
        adjMat=vcat(adjMat,zeros(1,size(adjMat,2)))
        adjMat=hcat(adjMat,zeros(size(adjMat,1),1))
        adjMatD=vcat(adjMatD,zeros(1,size(adjMatD,2)))
        adjMatD=hcat(adjMatD,zeros(size(adjMatD,1),1))
      end
    end
    #global trajTest=traj
    for i=1:size(outputNodes,1)
      x, y=outputNodes[i,1], outputNodes[i,2]
      #println([x,y])
      #array1[x,y]=1.0
      ind=getIndex(nod,x,y)
      #println(ind)
      adjMat[1,ind]=1.0
      adjMat[ind,1]=1.0
      adjMatD[1,ind]=diamMean[i]
      adjMatD[ind,1]=diamMean[i]
      #diamMean[i]==0 ? println("diam zero") : 1
    end

    nodesCp=copy(nod)
    for k=2:size(nodesCp)[1]
        outputNodes, diamMean, traj=algSearch(A,nod,numberOfSegPerNode[k],k,ASeg)
        for j=1:size(outputNodes,1)
          vere, indie = calDist(nod, outputNodes, k, j, traj)
          if vere
            outputNodes[j,:]=traj[indie,:]
            nod=vcat(nod,outputNodes[j,:]')
            append!(numberOfSegPerNode,2)
            adjMat=vcat(adjMat,zeros(1,size(adjMat,2)))
            adjMat=hcat(adjMat,zeros(size(adjMat,1),1))
            adjMatD=vcat(adjMatD,zeros(1,size(adjMatD,2)))
            adjMatD=hcat(adjMatD,zeros(size(adjMatD,1),1))
          end
        end
        #global trajTest=vcat(trajTest, traj)
        for i=1:size(outputNodes,1)
          x, y=outputNodes[i,1], outputNodes[i,2]
          #println([x,y])
          #array1[x,y]=1.0
          ind = getIndex(nod,x,y)
          #println(ind)
          adjMat[k,ind]=1
          adjMat[ind,k]=1
          adjMatD[k,ind]=diamMean[i]
          adjMatD[ind,k]=diamMean[i]
          #diamMean[i]==0 ? println("diam zero 2") : 1
        end
    end

    for k=size(nodesCp,1)+1:size(nod,1)
        outputNodes, diamMean, traj=algSearch(A,nod,numberOfSegPerNode[k],k,ASeg)
        # for j=1:size(outputNodes,1)
        #   vere, indie = calDist(nod, outputNodes, k, j, traj)
        #   if vere
        #     outputNodes[j,:]=traj[indie,:]
        #     nod=vcat(nod,outputNodes[j,:]')
        #     append!(numberOfSegPerNode,2)
        #     adjMat=vcat(adjMat,zeros(1,size(adjMat,2)))
        #     adjMat=hcat(adjMat,zeros(size(adjMat,1),1))
        #     adjMatD=vcat(adjMatD,zeros(1,size(adjMatD,2)))
        #     adjMatD=hcat(adjMatD,zeros(size(adjMatD,1),1))
        #   end
        # end
        #global trajTest=vcat(trajTest, traj)
        for i=1:size(outputNodes,1)
          x, y=outputNodes[i,1], outputNodes[i,2]
          #println([x,y])
          #array1[x,y]=1.0
          ind = getIndex(nod,x,y)
          #println(ind)
          adjMat[k,ind]=1
          adjMat[ind,k]=1
          adjMatD[k,ind]=diamMean[i]
          adjMatD[ind,k]=diamMean[i]
          #diamMean[i]==0 ? println("diam zero 3") : 1
        end
    end

    return nod, capturedNodes, adjMat, adjMatD
  end

  function measureDiam(aSeg, xold, yold, x, y)
    dr=[0 -1; 1 0]*[x-xold;y-yold]
    dmax=minimum(size(aSeg))
    d=0
    for i=1:dmax
      if aSeg[xold+dr[1]*i, yold+dr[2]*i]==0
        d=i
        break
      end
      if i==dmax
        @warn "Something is wrong when calculating the diameter. You're getting a diameter larger than the image"
      end
    end
    for i=1:dmax
      if aSeg[xold-dr[1]*i, yold-dr[2]*i]==0
        d=d+i-1
        break
      end
      if i==dmax
        @warn "something is wrong when calculating the diameter. You're getting a diameter larger than the image"
      end
    end

    return d
  end

  function calDist(nod, outnod, i, j, traj)
    d= sqrt( (nod[i,1]-outnod[j,1])^2 + (nod[i,2]-outnod[j,2])^2 )
    d2=0
    c=0;
    index=0
    thrs=1.1
    for k=1:size(traj,1)
      if traj[k,:] == nod[i,:]
        c=0
        index=k
      end
      if traj[k,:] == outnod[j,:]
        d2=c
        index = index +Int(floor((k - index)/2))
        break
      end
      if sum(abs.(traj[k+1,:]-traj[k,:]))==1
        c=c+1
      elseif sum(abs.(traj[k+1,:]-traj[k,:]))>1
        c=c+sqrt(2)
      end
    end
    #println(string(i)*"   "* string(d) *"   " * string(d2) * "  " * string(nod[i,:]) * "   " * string(outnod[j,:]))
    if d2/d > thrs && d>10
      return true, index
    else
      return false, index
    end
  end

  function arrayCleanup(aMat, aMatD, nod, ep)
    n=[0,0]
    newaMat, newaMatD, newNod, ep = [], [],[],[]
    for i=1:size(aMat,1)
      if sum(aMat[i,:])==2
        count=1
        for j=1:size(aMat,2)
          if aMat[i,j]==1
            n[count]=j
            count=count+1
          end
        end

        if aMat[n[1],n[2]]==1 && sqrt( (nod[n[1],1]-nod[n[2],1])^2 + (nod[n[1],2]-nod[n[2],2])^2) < 10
          newaMat = vcat(aMat[1:i-1,:], aMat[i+1:size(aMat,1),:])
          newaMat = hcat(newaMat[:,1:i-1], newaMat[:,i+1:size(aMat,2)])
          newaMatD = vcat(aMatD[1:i-1,:], aMatD[i+1:size(aMat,1),:])
          newaMatD = hcat(newaMatD[:,1:i-1], newaMatD[:,i+1:size(aMat,2)])
          newNod=vcat(nod[1:i-1,:], nod[i+1:size(nod,1),:])
          for k=1:size(ep,1)
            ep[k]>i ? ep[k]=ep[k]-1 : 1
          end
          break
        end
      end
    end
    return newaMat, newaMatD, newNod, ep
  end
  ############################################DATA MANIPULATION
  function overImposedPic(a,nodes, ep)
    Aep=zeros(size(a,1)+4,size(a,2)+4)
    Anod=zeros(size(a,1)+4,size(a,2)+4)

    for i=1:size(ep,1)
      Aep[ep[i,1]-2+2:ep[i,1]+2+2, ep[i,2]-2+2:ep[i,2]+2+2] .=1
    end

    for i=1:size(nodes,1)
      Anod[nodes[i,1]-2+2:nodes[i,1]+2+2, nodes[i,2]-2+2:nodes[i,2]+2+2] .=1
    end

    return colorview(RGB, a, Anod[3:end-2,3:end-2], Aep[3:end-2,3:end-2])
  end

  function readAndDir(arg2, arg)
    list=readdir("./" * string(arg2) * "/")
    if (isdir(string(arg))) == false
       mkdir(string(arg))
     end
     return list, arg
  end

  function saveFig(pathh, fig)
    pathname= string(pathh)
    @info fig, pathname
    savefig(fig,pathname)
  end

  function saveData(pathh, dat)
    pathname=string(pathh)
    writedlm(pathname,dat)
  end
end
######################################################

function walkTheWalk(thepath, readDir, writeDir)
  println("./"* string(readDir) * "/" * string(thepath) * "_Sket.tif")
    A=addZeros(init2("./"* string(readDir) * "/" * string(thepath) * "_Sket.tif"))
    ASeg=addZeros(init2("./"* string(readDir) * "/" * string(thepath) * ".gif"))
    #global diamMat=zeros(size(A,1),size(A,2))
    newA, nodes, arS, endpoints, numberOfSegPerNode=@time removeOverlapNodes(A )
    #colorview(RGB,A,arS,zeroarray)
    nodes, capNodes, adjMat, adjMatD=@time IterationOnAlgSearch(A,nodes,numberOfSegPerNode, ASeg)
    for i=1:size(adjMat,1)
      adjMat[i,i]=0
      adjMatD[i,i]=0
    end
    f(c)=round(c,digits=2)
    for i=1:size(adjMatD,1), j=1:size(adjMatD,2)
      if adjMat[i,j]>0 && adjMatD[i,j]==0
      end
    end

    positionEndpoints=[0]
    for i=1:size(nodes,1), j=1:size(endpoints,1)
      if endpoints[j,:]==nodes[i,:]
        append!(positionEndpoints,i-1)
      end
    end
    popfirst!(positionEndpoints)

    #adjMat,adjMatD,nodes,positionEndpoints=arrayCleanup(adjMat,adjMatD,nodes,positionEndpoints)

    loc2=[0 0];
    sstep=1;
    for i=1:sstep:size(A,1), j=1:sstep:size(A,1)
      if A[i,j]==1
        loc2=vcat(loc2,[j i])
      end
    end
    loc2=loc2[2:size(loc2)[1],:]

    d = scatter(loc2[:,1],-loc2[:,2].+size(A,1), markershapes = :square, markerstrokewidth=0, frame=:box, size=(500,500), legend=false)
    for i=1:size(nodes)[1], j=1:size(nodes)[1]
      if adjMat[j,i]==1
        d = plot!([nodes[j,2],nodes[i,2]],[size(A,1)-nodes[j,1],size(A,1)-nodes[i,1]],
              markershapes = :circle, lw=2, ms=5, legend=:none, c=:black, markerstrokewidth=0)
      end
    end
    fig=plot(d)
    saveFig(writeDir * "/" * thepath*"_network.png", fig)
    save(writeDir * "/" * thepath*"_Pic.png",overImposedPic(A,nodes, endpoints))
    # writedlm(thepath[1:10]*"_adj.csv",Array{Int,2}(adjMat))
    # writedlm(thepath[1:10]*"_nodes.csv",hcat(nodes,[10.0 for i=1:size(nodes,1)]))
    # writedlm(thepath[1:10]*"_endpoints.csv",Array{Int,1}(positionEndpoints))
    saveData(writeDir * "/" * thepath*"_adj.csv",Array{Int,2}(adjMat))
    saveData(writeDir * "/" * thepath*"_adjD.csv",Array{Float16,2}(f.(adjMatD)))
    saveData(writeDir * "/" * thepath*"_nodes.csv",hcat(nodes,[10.0 for i=1:size(nodes,1)]))
    saveData(writeDir * "/" * thepath*"_endpoints.csv",Array{Int,1}(positionEndpoints))
    return d, adjMat, positionEndpoints, nodes, adjMatD, endpoints
end

function theWholeProcess(readDir,writeDir)
  list, outDir=readAndDir( readDir, writeDir)
  # list=setdiff([list[i][1:10] for i=1:size(list,1)])
  list = setdiff([size(split(list[i],"_Sket"),1) == 2 ? split(list[i],"_Sket")[1] : split(list[i],".")[1] for i in 1:size(list,1)])

  for i in list
    d, adjMat, endpoints, nodes, adjMatD, ep=walkTheWalk(i, readDir,writeDir)

  end
  return 1

end

###########################################################
######################################################
#######################################################
#######################################################
######################################################
#######################################################
cd("/Users/javier/Desktop/SyntheticVasculature/Data")

thepath="01_manual1"
thepath = "01_dr_AVmanual"

@time d, adjMat, endpoints, nodes, adjMatD, ep=walkTheWalk(thepath, "sket", "networks")

# This function will go over all the skeletonized segmented .tif images
# and compute the Adj Matrix, end points
theWholeProcess("sket", "networks")

max_stars_repo_name: jquetzalcoatl/SyntheticVasculature

max_issues_repo_name: jquetzalcoatl/SyntheticVasculature

max_forks_repo_name: jquetzalcoatl/SyntheticVasculature


text:
"""
$SIGNATURES

Recursively process a conditional block from an opening HTML_COND_OPEN to a
{{end}}.
"""
function process_html_cond(hs::AS, qblocks::Vector{AbstractBlock},
                           i::Int)::Tuple{String,Int,Int}
    if i == length(qblocks)
        throw(HTMLBlockError("Could not close the conditional block " *
                             "starting with '$(qblocks[i].ss)'."))
    end

    init_idx      = i
    else_idx      = 0
    elseif_idx    = Vector{Int}()
    accept_elseif = true # false as soon as we see an else block
    content       = ""

    # inbalance keeps track of whether we've managed to find a
    # matching {{end}}. It increases if it sees other opening {{if..}}
    # and decreases if it sees a {{end}}
    inbalance = 1
    probe = qblocks[i] # just for allocation
    while i < length(qblocks) && inbalance > 0
        i    += 1
        probe = qblocks[i]
        # if we have a {{elseif ...}} or {{else}} with the right level,
        # keep track of them
        if inbalance == 1
            if probe isa HElseIf
                accept_elseif || throw(HTMLBlockError("Saw a {{elseif ...}} " *
                                            "after a {{else}} block."))
                push!(elseif_idx, i)
            elseif probe isa HElse
                else_idx      = i
                accept_elseif = false
            end
        end
        inbalance += hbalance(probe)
    end
    # we've exhausted the candidate qblocks and not found an appropriate {{end}}
    if inbalance > 0
        throw(HTMLBlockError("Could not close the conditional block " *
                             "starting with '$(qblocks[init_idx].ss)'."))
    end

    # we've found the closing {{end}} at index `i`
    β_close = probe
    i_close = i

    # We now have a complete conditional block with possibly
    # several else-if blocks and possibly an else block.
    # Check which one applies and discard the rest

    k   = 0 # index of the first verified condition
    lag = 0 # aux var to distinguish 1st block case

    # initial block may be {{if..}} or {{isdef..}} or...
    # if its not just a {{if...}}, need to act appropriately
    # and if the condition is verified then k=1
    βi = qblocks[init_idx]
    if βi isa Union{HIsDef,HIsNotDef,HIsPage,HIsNotPage}
        lag = 1
        if βi isa HIsDef
            k = Int(haskey(LOCAL_VARS, βi.vname))
        elseif βi isa HIsNotDef
            k = Int(!haskey(LOCAL_VARS, βi.vname))
        else
            # HIsPage//HIsNotPage
            rpath = splitext(unixify(locvar("fd_rpath")))[1]
            if FD_ENV[:STRUCTURE] < v"0.2"
                # current path is relative to /src/ for instance
                # /src/pages/blah.md -> pages/blah
                # if starts with `pages/`, replaces by `pub/`:
                # pages/blah => pub/blah
                rpath = replace(rpath, Regex("^pages") => "pub")
            end
            # compare with β.pages
            inpage = any(p -> match_url(rpath, p), βi.pages)

            if βi isa HIsPage
                k = Int(inpage)
            else
                k = Int(!inpage)
            end
        end
    end

    # If we've not yet found a verified condition, keep looking
    # Now all cond blocks ahead are {{if ...}} {{elseif  ...}}
    if iszero(k)
        if lag == 0
            conds = [βi.vname, (qblocks[c].vname for c in elseif_idx)...]
        else
            conds = [qblocks[c].vname for c in elseif_idx]
        end
        known = [haskey(LOCAL_VARS, c) for c in conds]

        if !all(known)
            idx = findfirst(.!known)
            throw(HTMLBlockError("At least one of the condition variable " *
                                 "could not be found: check '$(conds[idx])'."))
        end
        # check that all these variables are bool
        bools = [isa(locvar(c), Bool) for c  in conds]
        if !all(bools)
            idx = findfirst(.!bools)
            throw(HTMLBlockError("At least one of the condition variable is " *
                                 " not a Bool: check '$(conds[idx])'."))
        end
        # which one is verified?
        u = findfirst(c -> locvar(c), conds)
        k = isnothing(u) ? 0 : u + lag
    end

    if iszero(k)
        # if we still haven't found a verified condition, use the else
        # if there is one
        if !iszero(else_idx)
            # use elseblock, the content is from it to the β_close
            head = nextind(hs, to(qblocks[else_idx]))
            tail = prevind(hs, from(β_close))
            # get the stack of blocks in there
            action_qblocks = qblocks[else_idx+1:i_close-1]
            # process
            content = process_html_qblocks(hs, action_qblocks, head, tail)
        end
    else
        # determine the span of blocks
        oidx = 0 # opening conditional
        cidx = 0 # closing conditional
        if k == 1
            oidx = init_idx
        else
            oidx = elseif_idx[k-1]
        end
        if k ≤ length(elseif_idx)
            cidx = elseif_idx[k]
        else
            cidx = ifelse(iszero(else_idx), i_close, else_idx)
        end
        # Get the content
        head = nextind(hs, to(qblocks[oidx]))
        tail = prevind(hs, from(qblocks[cidx]))
        # Get the stack of blocks in there
        action_qblocks = qblocks[oidx+1:cidx-1]
        # process
        content = process_html_qblocks(hs, action_qblocks, head, tail)
    end
    # move the head after the final {{end}}
    head = nextind(hs, to(β_close))

    return content, head, i_close
end

"""
$SIGNATURES

Process a for block.
"""
function process_html_for(hs::AS, qblocks::Vector{AbstractBlock},
                          i::Int)::Tuple{String,Int,Int}
    # check that the iterable is known
    β_open = qblocks[i]
    vname  = β_open.vname
    iname  = β_open.iname
    if !haskey(LOCAL_VARS, iname)
        throw(HTMLBlockError("The iterable '$iname' is not recognised. " *
                             "Please make sure it's defined."))
    end

    # try to close the for loop
    if i == length(qblocks)
        throw(HTMLBlockError("Could not close the conditional block " *
                             "starting with '$(qblocks[i].ss)'."))
    end

    init_idx = i
    content  = ""

    # inbalance keeps track of whether we've managed to find a
    # matching {{end}}. It increases if it sees other opening {{if..}}
    # and decreases if it sees a {{end}}
    inbalance = 1
    while i < length(qblocks) && inbalance > 0
        i         += 1
        inbalance += hbalance(qblocks[i])
    end
    # we've exhausted the candidate qblocks and not found an appropriate {{end}}
    if inbalance > 0
        throw(HTMLBlockError("Could not close the conditional block " *
                             "starting with '$(qblocks[init_idx].ss)'."))
    end
    # we've found the closing {{end}} and index `i`
    β_close = qblocks[i]
    i_close = i

    isempty(locvar(iname)) && @goto final_step

    # is vname a single variable or multiple variables?
    # --> {{for v in iterate}}
    # --> {{for (v1, v2) in iterate }}
    vnames = [vname]
    if startswith(vname, "(")
        vnames = strip.(split(vname[2:end-1], ","))
    end
    # check that the first element of the iterate has the same length
    el1 = first(locvar(iname))
    length(vnames) in (1, length(el1)) ||
        throw(HTMLBlockError("In a {{for ...}}, the first element of" *
                "the iterate has length $(length(el1)) but tried to unpack" *
                "it as $(length(vnames)) variables."))

    # so now basically we have to simply copy-paste the content replacing
    # variable `vname` when it appears in a html block {{...}}
    # users should try not to be dumb about this... if vname or iname
    # corresponds to something they shouldn't touch, they'll crash things.

    # content of the for block
    inner = subs(hs, nextind(hs, to(β_open)), prevind(hs, from(β_close)))
    isempty(strip(inner)) && @goto final_step
    content = ""
    if length(vnames) == 1
        for value in locvar(iname)
            # at the moment we only consider {{fill ...}}
            content *= replace(inner,
                        Regex("({{\\s*fill\\s+$vname\\s*}})") => "$value")
        end
    else
        for value in locvar(iname)
            temp = inner
            for (vname, value) in zip(vnames, value)
                temp = replace(temp,
                        Regex("({{\\s*fill\\s+$vname\\s*}})") => "$value")
            end
            content *= temp
        end
    end
    @label final_step
    head = nextind(hs, to(β_close))
    return convert_html(content), head, i_close
end


text:
# Once this stabilizes, migrate to a Graphics layer? Only if that supports text, which seems unlikely.
using Colors
if VERSION < v"0.4.0-dev+3275"
    using Base.Graphics
else
    using Graphics
end

abstract AbstractAnnotation

# Use this type when you want your annotation to be linked to particular data-coordinates
# (for example, to highlight a particular data point)
# devicebb(data) returns a BoundingBox in device coordinates, userbb(data) returns one in user coordinates.
# It's up to the draw() function to decide how to exploit these (most commonly, with set_coords())
type AnchoredAnnotation{T} <: AbstractAnnotation
    devicebb::Function
    userbb::Function
    data::T
    valid::Bool
end
AnchoredAnnotation{T}(devicebb::Function, userbb::Function, data::T) =
    AnchoredAnnotation{T}(devicebb, userbb, data, true)

# Use this type when you want your annotation to appear at a specific point on the screen, regardless
# of zoom/resize state (e.g., a scale bar)
type FloatingAnnotation{T} <: AbstractAnnotation
    devicebb::Function
    data::T
end
# FloatingAnnotation{T}(devicebb::Function, data::T) = AnchoredAnnotation{T}(devicebb, data)

type AnnotationScalebarFixed{T}
    width::T   # Probably has units
    height::T
    getsize::Function   # syntax w,h = getsize(width,height)
    centerx::Float64
    centery::Float64
    color::Color
end
# AnnotationScalebar{T}(width::T, height::T, getsize::Function, centerx::Real, centery::Real, color::Color = RGB(1,1,1)) = AnnotationScalebar{T}(width, height, getsize, @compat(Float64(centerx)), @compat(Float64(centery)), color)

## Text annotations

type AnnotationText
    x::Float64
    y::Float64
    z::Float64
    t::Float64
    string::AbstractString
    color::Color
    fontfamily::String
    fontoptions::String
    fontsize::Integer
    fontdesc::String
    angle::Float64
    halign::String
    valign::String
    markup::Bool
    scale::Bool
end

function AnnotationText(x::Real, y::Real, str::AbstractString;
                        z = NaN, t = NaN,
                        color = RGB(0,0,0), angle = 0.0, fontfamily = "sans", fontsize = 10,
                        fontoptions = "",  halign = "center", valign = "center",
                        markup = false, scale=true)
    AnnotationText(@compat(Float64(x)), @compat(Float64(y)), @compat(Float64(z)),
                   @compat(Float64(t)), str, color, fontfamily, fontoptions,
                   fontsize, fontdescription(fontfamily, fontoptions, fontsize),
                   @compat(Float64(angle)), halign, valign, markup, scale)
end

fontdescription(fontfamily, fontoptions, fontsize) =
    string(fontfamily, " ", fontoptions, " ", fontsize)

## Point annotations

type AnnotationPoints{T}
    pts::T
    z::Float64
    t::Float64
    size::Float64
    shape::Char
    color::Color
    linewidth::Float64
    linecolor::Color
    scale::Bool
end

function AnnotationPoints{R<:Real}(xys::Vector{@compat(Tuple{R,R})}=@compat(Tuple{Float64,Float64})[];
                                   z = NaN, t = NaN, size=10.0, shape::Char='x', color = RGB(1,1,1),
                                   linewidth=1.0, linecolor=color, scale::Bool=false)
    AnnotationPoints{typeof(xys)}(xys, z, t, float(size), shape, to_colorant(color),
                                  float(linewidth), to_colorant(linecolor), scale)
end

function AnnotationPoints{R<:Real}(xys::Matrix{R}; z = NaN, t = NaN, size=10.0,
                                   shape::Char='x', color = RGB(1,1,1), linewidth=1.0,
                                   linecolor=color, scale::Bool=false)
    AnnotationPoints{Matrix{R}}(xys, z, t, float(size),
                                shape, to_colorant(color), float(linewidth),
                                to_colorant(linecolor), scale)
end

function AnnotationPoint(xy::(@compat Tuple{Real,Real}); z = NaN, t = NaN, size=10.0,
                         shape::Char='x', color = RGB(1,1,1), linewidth=1.0,
                         linecolor=color, scale::Bool=false)
    AnnotationPoints{@compat(Tuple{Float64,Float64})}((@compat(Float64(xy[1])),
                                                       @compat(Float64(xy[2]))),
                                                      z, t, float(size), shape,
                                                      to_colorant(color), float(linewidth),
                                                      to_colorant(linecolor), scale)
end

AnnotationPoint(x::Real, y::Real; args...) =
    AnnotationPoint((@compat(Float64(x)), @compat(Float64(y))); args...)

## Line annotations

@compat type AnnotationLines{R<:Union{Real, Tuple{Real, Real}}, T}
    lines::T
    z::Float64
    t::Float64
    linecolor::Color
    linewidth::Float64
    coordinate_order::Vector{Int}

    function AnnotationLines(lines::T, z, t, linecolor, linewidth, coord_order_str)
        ord = sortperm(coord_order_str.data)
        @assert coord_order_str[ord] == "xxyy"
        new(lines, z, t, linecolor, linewidth, ord)
    end
end

function AnnotationLines{R<:Real}(lines::Vector{(@compat Tuple{(@compat Tuple{R, R}),
                                                               (@compat Tuple{R, R})})}=
                                  (@compat Tuple{(@compat Tuple{Float64, Float64}),
                                                 (@compat Tuple{Float64, Float64})})[];
                                  z = NaN, t = NaN, color=RGB(1,1,1),
                                  linewidth=1.0, coord_order="xyxy")
    AnnotationLines{R,Vector{(@compat Tuple{(@compat Tuple{R, R}),
                                            (@compat Tuple{R, R})})}}(lines, z, t, color,
                                                                      linewidth, coord_order)
end

function AnnotationLines{R<:Real}(lines::Matrix{R}; z = NaN, t = NaN, color=RGB(1,1,1),
                                  linewidth=1.0, coord_order="xyxy")
    AnnotationLines{R, Matrix{R}}(lines, z, t, color, linewidth, coord_order)
end

function AnnotationLine{R<:Real}(line::(@compat Tuple{(@compat Tuple{R,R}),(@compat Tuple{R,R})});
                                 z = NaN, t = NaN, color=RGB(1,1,1), linewidth=1.0)
    AnnotationLines{R,(@compat Tuple{(@compat Tuple{R,R}),
                                     (@compat Tuple{R,R})})}(line, z, t, color, linewidth, "xyxy")
end

AnnotationLine(pt1::(@compat Tuple{Real,Real}), pt2::(@compat Tuple{Real,Real}); args...) =
    AnnotationLine((pt1, pt2); args...)

function AnnotationLine(c1::Real, c2::Real, c3::Real, c4::Real; coord_order="xyxy", args...)
    ord = sortperm(coord_order.data)
    @assert coord_order[ord] == "xxyy"
    (x1,x2,y1,y2) = [c1,c2,c3,c4][ord]
    AnnotationLine((@compat(Float64(x1)), @compat(Float64(y1))),
                   (@compat(Float64(x2)), @compat(Float64(y2))); args...)
end

## Box annotations

type AnnotationBox
    left::Float64
    top::Float64
    right::Float64
    bottom::Float64
    z::Float64
    t::Float64
    linecolor::Color
    linewidth::Float64
end

function AnnotationBox(c1::Real, c2::Real, c3::Real, c4::Real; z = NaN, t = NaN,
                       color=RGB(1,1,1), linewidth=1.0, coord_order="xyxy")
    ord = sortperm(coord_order.data)
    @assert coord_order[ord] == "xxyy"
    (x1, x2, y1, y2) = [c1, c2, c3, c4][ord]
    (x1, x2) = minmax(x1, x2)
    (y1, y2) = minmax(y1, y2)
    AnnotationBox(x1, y1, x2, y2, z, t, color, linewidth)
end

function AnnotationBox(pt1::(@compat Tuple{Real, Real}), pt2::(@compat Tuple{Real, Real});
                       coord_order="xyxy", args...)
    AnnotationBox(pt1..., pt2...; coord_order=coord_order, args...)
end

AnnotationBox(bb::BoundingBox; args...) = AnnotationBox(bb.xmin, bb.ymin, bb.xmax, bb.ymax; args...)

##############

setvalid!(ann::AnchoredAnnotation, z, t) = (ann.valid = annotation_isvalid(ann.data, z, t))

function annotation_isvalid(dat::@compat(Union{AnnotationText,
                                               AnnotationPoints,
                                               AnnotationLines,
                                               AnnotationBox}), z, t)
    (isnan(dat.z) || round(dat.z) == z) && (isnan(dat.t) || round(dat.t) == t)
end

annotation_isvalid(x, z, t) = true

function setvalid!(ann::FloatingAnnotation, z, t)
end

function draw(c::Canvas, ann::AnchoredAnnotation)
    if ann.valid
        ctx = getgc(c)
        Graphics.save(ctx)
        data = ann.data
        set_coords(ctx, ann.devicebb(data), ann.userbb(data))
        scale_x = width(ann.userbb(data))/width(ann.devicebb(data))
        scale_y = height(ann.userbb(data))/height(ann.devicebb(data))
        draw_anchored(ctx, data, scale_x, scale_y)
        restore(ctx)
    end
end

function draw{T}(c::Canvas, ann::FloatingAnnotation{AnnotationScalebarFixed{T}})
    ctx = getgc(c)
    Graphics.save(ctx)
    data = ann.data
    set_coords(ctx, ann.devicebb(data), BoundingBox(0,1,0,1))
    set_source(ctx, data.color)
    w, h = data.getsize(data.width, data.height)
    bb = BoundingBox(-w/2+data.centerx, w/2+data.centerx, -h/2+data.centery, h/2+data.centery)
    rectangle(ctx, bb)
    fill(ctx)
    restore(ctx)
end

function draw_anchored(ctx::CairoContext, data::AnnotationText, scale_x, scale_y)
    set_source(ctx, data.color)
    if data.scale && scale_x != 1
        fontdesc = fontdescription(data.fontfamily, data.fontoptions, round(Int,data.fontsize/scale_x))
    else
        fontdesc = data.fontdesc
    end
    Cairo.set_font_face(ctx, fontdesc)
    Cairo.text(ctx, data.x-0.5, data.y-0.5, data.string, halign = data.halign, valign = data.valign,
               angle = data.angle, markup = data.markup)
end

function draw_anchored(ctx::CairoContext, data::AnnotationPoints, scale_x, scale_y)
    set_line_width(ctx, data.linewidth)
    set_source(ctx, data.linecolor)
    if data.scale
        sz_x = sz_y = data.size
    else
        sz_x =  scale_x * data.size
        sz_y =  scale_y * data.size
    end
    draw_pts(ctx, data.pts, sz_x, sz_y, data.shape, data.color, data.linecolor)
end

draw_pts(ctx::CairoContext, pt::NTuple{2}, args...) = draw_pt(ctx, pt, args...)

function draw_pts{R<:Real}(ctx::CairoContext, pts::Vector{(@compat Tuple{R,R})}, args...)
    for pt in pts
        draw_pt(ctx, pt, args...)
    end
end

function draw_pts(ctx::CairoContext, pts::Matrix, args...)
    @assert size(pts,1) == 2
    for i in 1:size(pts,2)
        pt = pts[:,i]
        draw_pt(ctx, pt, args...)
    end
end

function draw_pt(ctx::CairoContext, pt, sz_x, sz_y, shape::Char, color::Color, linecolor::Color)
    x::Float64,y::Float64 = pt
    hsz_x = sz_x/2
    hsz_y = sz_y/2

    if (shape == '.') | (shape == 'o')
        move_to(ctx, x, y)
        if sz_x == sz_y
            new_sub_path(ctx)
            circle(ctx, x, y, sz_x)
        else
            # draw an ellipse
            Graphics.save(ctx)
            translate(ctx, x, y)
            scale(ctx, sz_x, sz_y)
            new_sub_path(ctx)
            circle(ctx, 0, 0, 1);
            restore(ctx)
        end
        if shape == '.'
            set_source(ctx, color)
            fill_preserve(ctx)
            set_source(ctx, linecolor)
        end
    elseif (shape == 'x') | (shape == '*') | (shape == '+')
        if (shape == 'x') | (shape == '*')
            move_to(ctx, x-hsz_x, y-hsz_y)
            line_to(ctx, x+hsz_x, y+hsz_y)
            move_to(ctx, x-hsz_x, y+hsz_y)
            line_to(ctx, x+hsz_x, y-hsz_y)
        end
        if (shape == '+') | (shape == '*')
            move_to(ctx, x-hsz_x, y)
            line_to(ctx, x+hsz_x, y)
            move_to(ctx, x, y-hsz_y)
            line_to(ctx, x, y+hsz_y)
        end
    end

    stroke(ctx)
end

function draw_anchored(ctx::CairoContext, data::AnnotationLines, args...)
    set_line_width(ctx, data.linewidth)
    set_source(ctx, data.linecolor)
    draw_lines(ctx, data.lines, data.coordinate_order)
end

draw_lines(ctx::CairoContext, line::(@compat Tuple{(@compat Tuple{Real, Real}),
                                                   (@compat Tuple{Real, Real})}), _) =
    draw_line(ctx, line)

function draw_lines{R<:Real}(ctx::CairoContext,
                             lines::Vector{(@compat Tuple{(@compat Tuple{R, R}),
                                                          (@compat Tuple{R, R})})}, _)
    for line in lines
        draw_line(ctx, line)
    end
end

function draw_lines{R<:Real}(ctx::CairoContext, lines::Matrix{R}, coordinate_order)
    for i in 1:size(lines, 2)
        pt = tuple(lines[coordinate_order,i]...)
        draw_line(ctx, pt)
    end
end

function draw_lines{R<:(@compat Tuple{Real,Real})}(ctx::CairoContext, lines::Matrix{R}, _)
    for i in 1:size(lines, 2)
        pt = tuple(lines[:,i]...)
        draw_line(ctx, pt)
    end
end

function draw_line(ctx::CairoContext, line::(@compat Tuple{(@compat Tuple{Real, Real}),
                                                           (@compat Tuple{Real, Real})}))
    (x1, y1), (x2, y2) = line
    move_to(ctx, x1, y1)
    line_to(ctx, x2, y2)
    stroke(ctx)
end

function draw_line(ctx::CairoContext, line::(@compat Tuple{Real, Real, Real, Real}))
    x1, y1, x2, y2 = line
    move_to(ctx, x1, y1)
    line_to(ctx, x2, y2)
    stroke(ctx)
end

## Box

function draw_anchored(ctx::CairoContext, data::AnnotationBox, args...)
    set_line_width(ctx, data.linewidth)
    set_source(ctx, data.linecolor)
    draw_box(ctx, data.top, data.bottom, data.left, data.right)
end

function draw_box(ctx::CairoContext, top, bottom, left, right)
    move_to(ctx, left, top)
    line_to(ctx, right, top)
    line_to(ctx, right, bottom)
    line_to(ctx, left, bottom)
    line_to(ctx, left, top)
    stroke(ctx)
end

to_colorant(c::Colorant) = c
to_colorant(str::AbstractString) = parse(Colorant, str)