text:
using Test
using ProbabilisticCircuits
using DataFrames
using LogicCircuits

@testset "ensemble tests with SamplePSDD" begin

    # Set up a logic constraint ϕ as a BDD and scope size n. Sample m PSDDs.
    function case(ϕ::Bdd, n::Integer, strategy::Symbol; m::Integer = 20, atol::Real = 1e-2)::Ensemble{StructProbCircuit}
        # All possible valuations (including impossible ones).
        M = all_valuations(collect(1:n))
        # Get only possible worlds.
        W = M[findall(ϕ.(eachrow(M))),:]
        # Assign random probabilities for each world in W.
        R = rand(1:20, size(W, 1))
        # Construct a dataset that maps the distribution of R (world W[i] repeats R[i] times).
        D = DataFrame(vcat([repeat(W[i,:], 1, R[i])' for i ∈ 1:size(W, 1)]...), :auto)
        # Learn ensemble of PSDDs from ϕ and D.
        E = ensemble_sample_psdd(m, ϕ, 16, D; strategy, pseudocount = 0.0, maxiter = 100,
                                 verbose = false, vtree_bias = 0.8)
        T = DataFrame(M, :auto)
        # Test consistency.
        @test (EVI(E, T) .> -Inf) == ϕ.(eachrow(M))
        # Test probabilities.
        evi = exp.(EVI(E, T))
        @test isapprox(evi[findall(>(0), evi)], (R/sum(R)); atol)
        return E
    end

    # some tests commented out to speed up tests
    Es = Vector{Ensemble{StructProbCircuit}}()
    for strategy ∈ [:likelihood, :uniform, :em, :stacking]
        case((1 ∧ 2) ∨ (3 ∧ ¬4) ∨ (¬1 ∧ 5), 5, strategy)
        # case((1 → 3) ∧ (5 → ¬2), 5, strategy)
        # case((1 ∧ 2 ∧ 3) ∨ (4 ∧ 5), 5, strategy)
        # case(exactly(3, collect(1:5)), 5, strategy)
        # case(atleast(3, collect(1:5)), 5, strategy)
        case(atmost(3, collect(1:5)), 5, strategy)
    end

    tmp = mktempdir()
    @testset "Saving and loading ensembles" begin
        for (i, E) ∈ enumerate(Es)
            @test_nowarn save_as_ensemble("$tmp/$i.esbl", E; quiet = true)
        end
    end
    Rs = Vector{Ensemble{StructProbCircuit}}()
    T = DataFrame(all_valuations(1:5), :auto)
    @testset "Loading ensembles" begin
        for i ∈ 1:length(Es)
            E = load_as_ensemble("$tmp/$i.esbl"; quiet = true)
            @test EVI(E, T) ≈ EVI(Es[i], T)
        end
    end
end


text:
mutable struct ImageReplay{R, SB, PRE, POST} <: AbstractReplay
    replay::R
    state_buffer::SB

    image_preproc::PRE
    image_postproc::POST
    function ImageReplay(replay::R,
                         state_buffer::SB,
                         image_preproc::PRE=identity,
                         image_postproc::POST=identity) where {R, SB, PRE, POST}
        new{R, SB, PRE, POST}(replay, state_buffer, image_preproc, image_postproc)
    end
end

function get_image(er::ImageReplay, s)
    er.state_buffer[s] |> er.image_postproc
end

function get_raw_image(er::ImageReplay, s)
    er.state_buffer[s]
end

function proc_state(er::ImageReplay, s)
    s |> er.image_preproc |> er.image_postproc
end

Base.length(er::ImageReplay) = length(er.replay)
function Base.getindex(er::ImageReplay, idx)
    experience = er.replay[idx]
    get_exp = (k, exp) -> begin
        if k == :s
            get_image(er, exp[k])
        elseif k == :sp
            get_image(er, exp[k])
        else
            exp[k]
        end
    end
    ks = keys(experience)
    nt_exp = (; zip(ks, [get_exp(k, experience) for k ∈ ks])...)
end
Base.view(er::ImageReplay, idx) = Base.view(er, idx)

function start_statebuffer!(er::ImageReplay, s)
    push!(er.state_buffer, er.image_preproc(s))
end

function Base.push!(er::ImageReplay, experience)
    s_idx = laststate(er.state_buffer)
    push!(er.state_buffer, er.image_preproc(experience.sp))
    sp_idx = laststate(er.state_buffer)
    ks = keys(experience)
    get_exp = (k, exp) -> begin
        if k == :s
            s_idx
        elseif k == :sp
            sp_idx
        else
            exp[k]
        end
    end
    nt_exp = (; zip(ks, [get_exp(k, experience) for k ∈ ks])...)
    push!(er.replay, nt_exp)
end

function get_expr(er::ImageReplay{<:AbstractSequenceReplay}, k, expr)
    if k == :s
        get_state(e) = get_image(er, e.s)
        Flux.batchseq([[get_state.(seq); [get_image(er, seq[end].sp)]] for seq in expr], zero(get_state(expr[1][1])))
    elseif k == :sp
        [get_raw_image(er, seq[end].sp) for seq ∈ expr]
    elseif k == :am1
        l = maximum(length.(expr)) + 1
        temp = [rpad([[seqi_j.am1[] for seqi_j ∈ seq]; [seq[end].a[]]], l, one(expr[1][1].a[])) for seq in expr]
        [[temp[b][t] for b ∈ 1:length(temp)] for t ∈ 1:length(temp[1])]
    elseif length(expr[1][1][k]) > 1
        temp = [[seqi_j[k] for seqi_j ∈ seq] for seq in expr]
    else
        temp = [[seqi_j[k][] for seqi_j ∈ seq] for seq in expr]
    end
end

# lets just go ahead and deal w/ batching here.
function make_experience(er::ImageReplay{<:AbstractSequenceReplay}, expr)
    ks = keys(expr[1][1])
    data = [get_expr(er, k, expr) for k ∈ ks]
    length.(expr), (; zip(ks, data)...)
end

function make_experience(er::ImageReplay, experience)
    throw("Not Implemented!")
end

sample(er::ImageReplay, batch_size, args...) = sample(Random.GLOBAL_RNG, er, batch_size, args...)
function sample(rng::Random.AbstractRNG, er::ImageReplay, batch_size, args...)
    exp_s_idx, experience = sample(rng, er.replay, batch_size, args...)
    exp_s_idx, make_experience(er, experience)
end

ext: jl


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

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

"""
    next_did()

Produces an incrementing ID that will be used for DArrays.
"""
next_did

"""
    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}}

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

## core constructors ##

function 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

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 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
# The refs must have the same layout as the parts distributed.
# i.e.
#    size(refs) must specify the distribution of dimensions across processors
#    prod(size(refs)) must equal number of parts
# FIXME : Empty parts are currently not supported.
function DArray(refs)
    dimdist = size(refs)
    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)]...)

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

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

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

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(ref, identity)
    global registry
    lp = fetch(ref)
    registry[(identity, :LOCALPART)] = lp
    return size(lp)
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))
    np = prod(dist)
    procs_used = procs[1:np]
    idxs, _ = chunk_idxs([size(A)...], dist)

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

"""
    distribute(A, DA)

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

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

    s = verified_destination_serializer(procs(DA), size(DA.indexes)) do pididx
        A[DA.indexes[pididx]...]
    end
    return DArray(I->localpart(s), DA)
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(size(SD), 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}}...) = view(d, I...)

Base.copy!(dest::SubOrDArray, src::SubOrDArray) = begin
    asyncmap(procs(dest)) do p
        remotecall_fetch(p) do
            localpart(dest)[:] = src[localindexes(dest)...]
        end
    end
    return dest
end

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

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

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

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

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

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

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

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

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


text:
using SCS
const SCS_PRIMAL_FACTORY = MOI.OptimizerWithAttributes(SCS.Optimizer, MOI.Silent() => true)
const SCS_DUAL_FACTORY = dual_optimizer(SCS_PRIMAL_FACTORY)
const SCS_PRIMAL_OPT = MOI.instantiate(SCS_PRIMAL_FACTORY)
const SCS_DUAL_OPT = MOI.instantiate(SCS_DUAL_FACTORY)

push!(primal_power_cone_factory, SCS_PRIMAL_FACTORY)
push!(dual_power_cone_factory, SCS_DUAL_FACTORY)
push!(dual_power_cone_optimizer, SCS_DUAL_OPT)
push!(primal_power_cone_optimizer, SCS_PRIMAL_OPT)

@testset "SCS Exponential Cone Problems" begin
    list_of_exp_problems = [
        exp1_test,
        exp2_test
    ]
    test_strong_duality(list_of_exp_problems, SCS_PRIMAL_FACTORY)
end

@testset "SCS Power Cone Problems" begin
    list_of_pow_problems = [
        pow1_test,
        pow2_test
    ]
    test_strong_duality(list_of_pow_problems, SCS_PRIMAL_FACTORY; atol = 1e-3, rtol = 1e-3)
end


text:
@polyvar x y
p = x^4*y^2 + x^2*y^4 - 3*x^2*y^2 + 1
@constraint(model, p >= 0)
optimize!(model)
# Solution status is `Infeasible`
@show termination_status(model)

max_stars_repo_name: sums-of-squares/tools

max_issues_repo_name: sums-of-squares/tools

max_forks_repo_name: sums-of-squares/tools


text:
include("ring.jl")
using .rings

ar = rings.aring(1.0e20,10,9,2,2)

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
println(Base.stderr, "Starting notebook kernel server")

pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", "packages"))
using VSCodeServer
popfirst!(LOAD_PATH)

println(Base.stderr, "Core notebook support loaded")

using InteractiveUtils

let
    outputchannel_logger = Base.CoreLogging.SimpleLogger(Base.stderr)

    Base.with_logger(outputchannel_logger) do
        @info "Processing command line arguments..."
    end

    args = [popfirst!(Base.ARGS) for _ in 1:2]

    conn_pipeline, telemetry_pipeline = args[1:2]

    Base.with_logger(outputchannel_logger) do
        @info "Command line arguments processed"
    end

    ccall(:jl_exit_on_sigint, Nothing, (Cint,), false)

    Base.with_logger(outputchannel_logger) do
        @info "Handing things off to VSCodeServer.serve_notebook"
    end

    VSCodeServer.serve_notebook(conn_pipeline, outputchannel_logger, crashreporting_pipename=telemetry_pipeline)
end

max_stars_repo_name: TanookiToad/julia-vscode

max_issues_repo_name: TanookiToad/julia-vscode

max_forks_repo_name: TanookiToad/julia-vscode


text:
function get_date(obj::TimeStamped)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_observed_value(obj::ComparableMeasurement)
    return jcall(obj, "getObservedValue", Vector{jdouble}, ())
end


text:
function SeisNMO(in;dt=0.001,offset=1000.,tnmo=[0.],vnmo=[1500.],max_stretch=1000)

	nt,nx = size(in)
	if length(offset) < size(in,2)
		offset = offset[1]*fill!(similar(in[1,:]),one(eltype(in[1,:])))
	end

	# interpolate tau,v pairs to match sampling of input data
	if (length(vnmo) == 1)
		tnmo = convert(Float64,tnmo[1]);
		vnmo = convert(Float64,vnmo[1]);
		ti = collect(0:1:nt-1)*dt
		vi = ones(1,nt)*vnmo

	else
		tnmo = (convert(Array{Float64,1},vec(tnmo)),)
		vnmo = convert(Array{Float64,1},vec(vnmo))
		ti = collect(0:1:nt-1)*dt
#g = InterpIrregular(tnmo, vnmo, BCnan, InterpLinear)
		g = interpolate(tnmo, vnmo, Gridded(Linear()))
		ge = extrapolate(g, Line())
		vi = ge(ti)
	end
	out = zeros(typeof(in[1,1]),size(in))
	M = zeros(nt,1)
	for it = 1:nt
		for ix = 1:nx
			time = sqrt(ti[it]^2 + (offset[ix]/vi[it]).^2)
			stretch = (time-ti[it])/(ti[it]+1e-10)
			if (stretch<max_stretch/100)
				its = round(Int,time/dt)+1
				it1 = round(Int,floor(time/dt))+1
				it2 = it1+1
				a = its-it1
				if (it2 <= nt)
					out[it,ix] = (1-a)*in[it1,ix]+a*in[it2,ix]
				end
			end
		end
	end
	return out
end

"""
    SeisNMO(in,out,parameters; <keyword arguments>)

# Arguments
- `in::String`: Input file - Seis format.
- `out::String`: Output file - Seis format.
- `parameters` : list of the keyword arguments for the function SeisMute.

# Keyword arguments
- `group="gather"` : Options are all, some or gather
- `key=["imx","imy"]` : Defines type of gather
- `itrace=1` : Initial trace number
- `ntrace=10000` : Total number of traces to process at once

# Example
```
julia> param = Dict(:tmute=>0.0, :vmute=>10000, :taper=>0.05,:dt=>0.01)
julia> SeisMute(filein,fileout, param,group="some")
```

"""
function SeisNMO(in::String,out::String,parameters;group="gather",key=["imx","imy"],itrace=1,ntrace=10000)

	dt = get(parameters,:dt,0.001)
	offset = get(parameters,:offset,1000)
	tnmo = get(parameters,:tnmo,[0.])
	vnmo = get(parameters,:vnmo,[1500.])
	max_stretch = get(parameters,:max_stretch,1000)

	if (group=="all")
		headers = SeisMain.SeisReadHeaders(in);
		offset = SeisMain.ExtractHeader(headers,"h")
		dt = headers[1].d1
		parameters = Dict(:offset=>offset,:tnmo=>tnmo,:vnmo=>vnmo,:max_stretch=>max_stretch,:dt=>dt)

		SeisProcessFile(in,out,[SeisNMO],[parameters];group=group)
	else
		itrace_in = itrace

		nx = SeisMain.GetNumTraces(in)
		while itrace_in <= nx
			h = SeisMain.SeisReadHeaders(in,group=group,key=key,itrace=itrace_in,ntrace=ntrace)
			offset = SeisMain.ExtractHeader(h,"h")
			dt = h[1].d1

			parameters = Dict(:offset=>offset,:tnmo=>tnmo,:vnmo=>vnmo,:max_stretch=>max_stretch,:dt=>dt)
			#SeisProcessFile(in,out,[SeisNMO],[parameters];group=group,key=key,itrace=itrace_in,ntrace=ntrace)
			d1,h1,e1 = SeisMain.SeisRead(in,group=group,key=key,itrace=itrace_in,ntrace=ntrace)

			d2 = SeisNMO(d1;parameters...)
			d1 = copy(d2)

			SeisMain.SeisWrite(out,d1,h1,e1,itrace=itrace_in)
			itrace_in += ntrace
		end
	end
end


text:
using OrdinaryDiffEq
using LinearAlgebra
using NLsolve
using Test

p_inv = [500.0
    0.084
    4.69
    2.0
    400.0
    20.0
    0.2
    1000.0
    0.59
    736.0
    0.0
    0.0
    0.2
    1.27
    14.3
    0.0
    50.0
    0.2
    0.08
    0.003
    0.074
    0.2
    0.01 # Matt
    1.0  # V_ref
    1.0  # ω ref
    0.7  # Pref
    0    # Q_ref
    0.5  #Xtrans
    0.0  # Rtrans
    1.01 #Vm
    0.0] # Vθ

function vsm(dx, x, p, t)
    #PARAMETERS
    ω_lp = p[1]
    kp_pll = p[2]
    ki_pll = p[3]
    Ta = p[4]
    kd = p[5]
    kω = p[6]
    kq = p[7]
    ωf = p[8]
    kpv = p[9]
    kiv = p[10]
    kffv = p[11]
    rv = p[12]
    lv = p[13]
    kpc = p[14]
    kic = p[15]
    kffi = p[16]
    ωad = p[17]
    kad = p[18]
    lf = p[19]
    rf = p[20]
    cf = p[21]
    lg = p[22]
    rg = p[23]
    Vref = p[24]
    ωref = p[25]
    Pref = p[26]
    Qref = p[27]
    Xtrans = p[28]
    Rtrans = p[29]
    Vm = p[30]
    Vθ = p[31]

    #STATE INDEX AND STATES
    i__vi_filter, vi_filter = 1, x[1]
    i__γd_ic, γd_ic = 2, x[2]
    i__vq_pll, vq_pll = 3, x[3]
    i__γq_ic, γq_ic = 4, x[4]
    i__ir_filter, ir_filter = 5, x[5]
    i__ξd_ic, ξd_ic = 6, x[6]
    i__ϕd_ic, ϕd_ic = 7, x[7]
    i__ε_pll, ε_pll = 8, x[8]
    i__ir_cnv, ir_cnv = 9, x[9]
    i__vr_filter, vr_filter = 10, x[10]
    i__ω_oc, ω_oc = 11, x[11]
    i__ξq_ic, ξq_ic = 12, x[12]
    i__vd_pll, vd_pll = 13, x[13]
    i__q_oc, q_oc = 14, x[14]
    i__ϕq_ic, ϕq_ic = 15, x[15]
    i__θ_pll, θ_pll = 16, x[16]
    i__θ_oc, θ_oc = 17, x[17]
    i__ii_cnv, ii_cnv = 18, x[18]
    i__ii_filter, ii_filter = 19, x[19]
    i__ir_source, ir_source = 20, x[20]
    i__ii_source, ii_source = 21, x[21]

    ω_base = 60.0 * 2 * pi
    ω_sys = 1.0

    #PLL
    δω_pll = kp_pll * atan(vq_pll / vd_pll) + ki_pll * ε_pll
    ω_pll = δω_pll + ω_sys
    vd_filt_pll = sin(θ_pll + pi / 2) * vr_filter - cos(θ_pll + pi / 2) * vi_filter
    vq_filt_pll = cos(θ_pll + pi / 2) * vr_filter + sin(θ_pll + pi / 2) * vi_filter

    dx[i__vd_pll] = ω_lp * (vd_filt_pll - vd_pll)                                 #docs:(1a)
    dx[i__vq_pll] = ω_lp * (vq_filt_pll - vq_pll)                                 #docs:(1b)
    dx[i__ε_pll] = atan(vq_pll / vd_pll)                                         #docs:(1c)
    dx[i__θ_pll] = ω_base * δω_pll                                             #docs:(1d)

    #OUTER LOOP CONTROL
    pe = vr_filter * ir_filter + vi_filter * ii_filter                              #docs:(2d)
    qe = vi_filter * ir_filter - vr_filter * ii_filter                              #docs:(2e)
    v_ref_olc = Vref + kq * (Qref - q_oc)

    dx[i__ω_oc] = (Pref - pe - kd * (ω_oc - ω_pll) - kω * (ω_oc - ωref)) / Ta        #docs:(2a)
    dx[i__θ_oc] = ω_base * (ω_oc - ω_sys)                                            #docs:(2b)
    dx[i__q_oc] = ωf * (qe - q_oc)                                               #docs:(2c)

    #INNER LOOP CONTROL
    #reference transormations
    vd_filt_olc = sin(θ_oc + pi / 2) * vr_filter - cos(θ_oc + pi / 2) * vi_filter
    vq_filt_olc = cos(θ_oc + pi / 2) * vr_filter + sin(θ_oc + pi / 2) * vi_filter
    id_filt_olc = sin(θ_oc + pi / 2) * ir_filter - cos(θ_oc + pi / 2) * ii_filter
    iq_filt_olc = cos(θ_oc + pi / 2) * ir_filter + sin(θ_oc + pi / 2) * ii_filter
    id_cnv_olc = sin(θ_oc + pi / 2) * ir_cnv - cos(θ_oc + pi / 2) * ii_cnv
    iq_cnv_olc = cos(θ_oc + pi / 2) * ir_cnv + sin(θ_oc + pi / 2) * ii_cnv

    #Voltage control equations
    Vd_filter_ref = v_ref_olc - rv * id_filt_olc + ω_oc * lv * iq_filt_olc      #docs:(3g)
    Vq_filter_ref = -rv * iq_filt_olc - ω_oc * lv * id_filt_olc                 #docs:(3h)
    dx[i__ξd_ic] = Vd_filter_ref - vd_filt_olc                                 #docs:(3a)
    dx[i__ξq_ic] = Vq_filter_ref - vq_filt_olc                                 #docs:(3b)

    #current control equations
    Id_cnv_ref = (
        kpv * (Vd_filter_ref - vd_filt_olc) + kiv * ξd_ic -                           #docs:(3i)
        cf * ω_oc * vq_filt_olc + kffi * id_filt_olc
    )
    Iq_cnv_ref = (
        kpv * (Vq_filter_ref - vq_filt_olc) +
        kiv * ξq_ic +                           #docs:(3j)
        cf * ω_oc * vd_filt_olc +
        kffi * iq_filt_olc
    )
    dx[i__γd_ic] = Id_cnv_ref - id_cnv_olc                                      #docs:(3c)
    dx[i__γq_ic] = Iq_cnv_ref - iq_cnv_olc                                      #docs:(3d)

    #active damping equations
    Vd_cnv_ref = (
        kpc * (Id_cnv_ref - id_cnv_olc) + kic * γd_ic - lf * ω_oc * iq_cnv_olc +          #docs:(3k)
        kffv * vd_filt_olc - kad * (vd_filt_olc - ϕd_ic)
    )
    Vq_cnv_ref = (
        kpc * (Iq_cnv_ref - iq_cnv_olc) +
        kic * γq_ic +
        lf * ω_oc * id_cnv_olc +          #docs:(3l)
        kffv * vq_filt_olc - kad * (vq_filt_olc - ϕq_ic)
    )
    dx[i__ϕd_ic] = ωad * (vd_filt_olc - ϕd_ic)                                 #docs:(3e)
    dx[i__ϕq_ic] = ωad * (vq_filt_olc - ϕq_ic)                                  #docs:(3f)

    #LCL FILTER
    #reference transformations
    Vr_cnv = sin(θ_oc + pi / 2) * Vd_cnv_ref + cos(θ_oc + pi / 2) * Vq_cnv_ref
    Vi_cnv = -cos(θ_oc + pi / 2) * Vd_cnv_ref + sin(θ_oc + pi / 2) * Vq_cnv_ref

    Vr_pcc = Vm * cos(Vθ)
    Vi_pcc = Vm * sin(Vθ)

    dx[i__ir_cnv] =
        (ω_base / lf) *                                                #docs:(5a)
        (Vr_cnv - vr_filter - rf * ir_cnv + ω_sys * lf * ii_cnv)
    dx[i__ii_cnv] =
        (ω_base / lf) *                                                #docs:(5b)
        (Vi_cnv - vi_filter - rf * ii_cnv - ω_sys * lf * ir_cnv)
    dx[i__vr_filter] =
        (ω_base / cf) *                                             #docs:(5c)
        (ir_cnv - ir_filter + ω_sys * cf * vi_filter)
    dx[i__vi_filter] =
        (ω_base / cf) *                                             #docs:(5d)
        (ii_cnv - ii_filter - ω_sys * cf * vr_filter)
    dx[i__ir_filter] =
        (ω_base / lg) *                                             #docs:(5e)
        (vr_filter - Vr_pcc - rg * ir_filter + ω_sys * lg * ii_filter)
    dx[i__ii_filter] =
        +(ω_base / lg) *                                            #docs:(5f)
        (vi_filter - Vi_pcc - rg * ii_filter - ω_sys * lg * ir_filter)

    # Network interface algebraic equations 0 = I - YV
    line_currents = ((Vr_cnv + Vi_cnv * 1im) - (Vr_pcc + Vi_pcc * 1im)) / (Rtrans + Xtrans * 1im)
    dx[i__ii_source] = ii_source - imag(line_currents)
    dx[i__ir_source] = ir_source - real(line_currents)
    return
end

u0 = [0.0,
    0.01,
    0.01,
    0.01,
    0.01,
    0.01,
    0.01,
    0.01,
    1.0,
    1.0,
    0.0,
    0.01,
    0.01,
    0.01,
    0.01,
    0.01,
    0.01,
    0.01,
    0.0,
    1.0,
    0.0,
]
init_f! = (dx, x) -> vsm(dx, x, p_inv, 0)
res = nlsolve(init_f!, u0)

M = diagm(ones(21))
#Last 2 equations are algebraic
M[20, 20] = M[21, 21] = 0.0

f = ODEFunction(vsm, mass_matrix=M)

condition(u, t, integrator) = t == 1.0
affect!(integrator) = integrator.p[28] += 0.2
cb = DiscreteCallback(condition, affect!)

prob = ODEProblem(f, deepcopy(res.zero), (0, 20.0), deepcopy(p_inv))
refsol = solve(prob, Rodas4(), saveat = 0.1, callback=cb, tstops=[1.0], reltol=1e-12, abstol=1e-18)

for solver in (Rodas4, Rodas4P, Rodas5, Rodas5P, FBDF, QNDF, Rosenbrock23)
    @show solver
    prob = ODEProblem(f, deepcopy(res.zero), (0, 20.0), deepcopy(p_inv))
    sol = solve(prob, solver(), saveat=0.1, callback=cb, tstops=[1.0], reltol=1e-12, abstol=1e-16)
    @test sol.retcode == :Success
    @test sol.t[end] == 20.0
    @test maximum(sol - refsol) < 1e-11
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()
VERSION < v"0.4.0" && (AbstractString=String)

module RawFile

using IniFile
using ProgressMeter

import Base: size, ndims, read, get, setindex!, getindex

export Rawfile, read_raw, write_raw

type Rawfile
  filename::AbstractString
  extRaw::AbstractString
  extHeader::AbstractString
  dtype
  size

  function Rawfile(filename::AbstractString; extRaw=".dat", extHeader=".srm")
    s = split(filename,".")
    if length(s) == 1
      fileWithoutExt = filename
    else
      pop!(s)
      fileWithoutExt = *(s...)
    end
    f = new(fileWithoutExt, extRaw, extHeader)
    if filesize(f.filename*f.extHeader) > 0
      read_srm(f)
    end
    f
  end

  function Rawfile(filename::AbstractString, dtype, size; extRaw=".dat", extHeader=".srm")
    f = Rawfile(filename, extRaw = extRaw, extHeader = extHeader)
    f.size = size
    f.dtype = dtype
    f
  end
end

function read_raw(filename::AbstractString; extRaw=".dat", extHeader=".srm")
  f = Rawfile(filename, extRaw = extRaw, extHeader = extHeader)
  read_srm(f)
  f[]
end

function write_raw(filename::AbstractString, x; extRaw=".dat", extHeader=".srm")
  f = Rawfile(filename, extRaw = extRaw, extHeader = extHeader)
  f[] = x
end

function read_srm(f::Rawfile)
  ini = Inifile()
  read(ini, f.filename*f.extHeader)
  f.size = int(split(get(ini,"size"),","))
  # The following is due to type aliases
  dtypestr = get(ini,"datatype")
  if dtypestr == "Complex{Float64}"
    f.dtype = Complex128
  elseif dtypestr == "Complex{Float32}"
    f.dtype = Complex64
  else
    f.dtype = eval(symbol( get(ini,"datatype") ))
  end
end

function write_srm(f::Rawfile)
  ini = Inifile()
  set(ini, "size", string(f.size)[2:end-1] )
  set(ini, "datatype", string(f.dtype) )
  open(f.filename*f.extHeader,"w") do fd
    write(fd, ini)
  end
end

ndims(f::Rawfile) = length(f.size)
size(f::Rawfile) = f.size
size(f::Rawfile, dir::Integer) = f.size[dir]

function setindex!(f::Rawfile,x)
  f.size = size(x)
  f.dtype = eltype(x)
  open(f.filename*f.extRaw,"w") do fd
    write(fd, x)
  end
  write_srm(f)
end

function getindex(f::Rawfile)
  x = open(f.filename*f.extRaw,"r") do fd
    read(fd, f.dtype, f.size...)
  end
end

function getindex(f::Rawfile,::Colon,::Colon)
   f[]
end

function getindex(f::Rawfile,x::Colon,args...)
  f[1:size(f, 1),args...]
end

function getindex(f::Rawfile,x::UnitRange,y)
  filename = f.filename*f.extRaw
  fd = open(filename,"r")

  matrix = zeros(f.dtype, (length(x),length(y)))
  p = Progress(length(y), 1, "Loading data from $filename...")
  for l=1:length(y)
    seek(fd, ((y[l]-1)*f.size[1] + x[1] - 1 )*sizeof(f.dtype))
    matrix[:,l] = read(fd, f.dtype, length(x))
    next!(p)
  end

  close(fd)
  matrix
end

function getindex(f::Rawfile,x::UnitRange, y, z)
  filename = f.filename*f.extRaw
  fd = open(filename,"r")

  data = zeros(f.dtype, (length(x),length(y),length(z)))
  p = Progress(length(z), 1, "Loading data from $filename...")
  for r=1:length(z)
    for l=1:length(y)
      seek(fd, (((z[r]-1)*f.size[2] + (y[l]-1))*f.size[1] + x[1] - 1 )*sizeof(f.dtype))
      data[:,l,r] = read(fd, f.dtype, length(x))
    end
    next!(p)
  end

  close(fd)
  data
end

function setindex!(f::Rawfile, A,::Colon,::Colon)
  f[] = A
end

function setindex!(f::Rawfile,x::Colon,args...)
  setindex!(f,1:size(f, 1),args...)
end

function setindex!(f::Rawfile, A, x::UnitRange, y)
  if f.dtype != eltype(A); error("wrong datatype") end
  if length(x) != size(A,1); error("dimension missmatch") end

  open(f.filename*f.extRaw,"a+") do fd
    for l=1:length(y)
      seek(fd, ((y[l]-1)*f.size[1] + x[1] - 1 )*sizeof(f.dtype))
      write(fd, A[x,l])
    end
  end
end

end # module


text:
export tuple_with_left, tuple_with_right

import Base: show

tuple_with_left(value::T) where T = TupleWithLeftOperator{T}(value)

struct TupleWithLeftOperator{T} <: InferableOperator
    value :: T
end

function on_call!(::Type{L}, ::Type{Tuple{T, L}}, operator::TupleWithLeftOperator{T}, source) where T where L
    return proxy(Tuple{T, L}, source, TupleWithLeftProxy{T}(operator.value))
end

operator_right(operator::TupleWithLeftOperator{T}, ::Type{L}) where { T, L } = Tuple{T, L}

struct TupleWithLeftProxy{T} <: ActorProxy
    value :: T
end

actor_proxy!(::Type{Tuple{T, L}}, proxy::TupleWithLeftProxy{T}, actor::A) where { T, L, A } = TupleWithLeftActor{T, L, A}(proxy.value, actor)

struct TupleWithLeftActor{T, L, A} <: Actor{L}
    value :: T
    actor :: A
end

on_next!(actor::TupleWithLeftActor{T, L}, data::L) where { T, L } = next!(actor.actor, (actor.value, data))
on_error!(actor::TupleWithLeftActor, err)                         = error!(actor.actor, err)
on_complete!(actor::TupleWithLeftActor)                           = complete!(actor.actor)

Base.show(io::IO, ::TupleWithLeftOperator)                   = print(io, "TupleWithLeftOperator()")
Base.show(io::IO, ::TupleWithLeftProxy{T})    where { T    } = print(io, "TupleWithLeftProxy()")
Base.show(io::IO, ::TupleWithLeftActor{T, L}) where { T, L } = print(io, "TupleWithLeftActor($L -> Tuple{$T, $L})")

tuple_with_right(value::T) where T = TupleWithRightOperator{T}(value)

struct TupleWithRightOperator{T} <: InferableOperator
    value :: T
end

function on_call!(::Type{L}, ::Type{Tuple{L, T}}, operator::TupleWithRightOperator{T}, source) where T where L
    return proxy(Tuple{L, T}, source, TupleWithRightProxy{T}(operator.value))
end

operator_right(operator::TupleWithRightOperator{T}, ::Type{L}) where { T, L } = Tuple{L, T}

struct TupleWithRightProxy{T} <: ActorProxy
    value :: T
end

actor_proxy!(::Type{Tuple{L, T}}, proxy::TupleWithRightProxy{T}, actor::A) where { T, L, A } = TupleWithRightActor{T, L, A}(proxy.value, actor)

struct TupleWithRightActor{T, L, A} <: Actor{L}
    value :: T
    actor :: A
end

on_next!(actor::TupleWithRightActor{T, L}, data::L) where { T, L } = next!(actor.actor, (data, actor.value))
on_error!(actor::TupleWithRightActor, err)                         = error!(actor.actor, err)
on_complete!(actor::TupleWithRightActor)                           = complete!(actor.actor)

Base.show(io::IO, ::TupleWithRightOperator)                   = print(io, "TupleWithRightOperator()")
Base.show(io::IO, ::TupleWithRightProxy{T})    where { T    } = print(io, "TupleWithRightProxy()")
Base.show(io::IO, ::TupleWithRightActor{T, L}) where { T, L } = print(io, "TupleWithRightActor($L -> Tuple{$L, $T})")


text:
#### Deprecate on 0.5 (to be removed on 0.6)

function dim(d::MultivariateDistribution)
    Base.depwarn("dim(d::MultivariateDistribution) is deprecated. Please use length(d).", :dim)
    return length(d)
end

function binaryentropy(d::UnivariateDistribution)
    Base.depwarn("binaryentropy is deprecated. Please use entropy(d, 2).", :binaryentropy)
    return entropy(d) / log(2)
end

@Base.deprecate logpmf logpdf
@Base.deprecate logpmf! logpmf!
@Base.deprecate pmf pdf

#### Deprecate on 0.6 (to be removed on 0.7)

@Base.deprecate expected_logdet meanlogdet

function probs(d::DiscreteUnivariateDistribution)
    Base.depwarn("probs(d::$(typeof(d))) is deprecated. Please use pdf(d) instead.", :probs)
    return probs(d)
end

function Binomial(n::Real, p::Real)
    Base.depwarn("Binomial(n::Real, p) is deprecated. Please use Binomial(n::Integer, p) instead.", :Binomial)
    Binomial(Int(n), p)
end

function Binomial(n::Real)
    Base.depwarn("Binomial(n::Real) is deprecated. Please use Binomial(n::Integer) instead.", :Binomial)
    Binomial(Int(n))
end

function BetaBinomial(n::Real, α::Real, β::Real)
    Base.depwarn("BetaBinomial(n::Real, α, β) is deprecated. Please use BetaBinomial(n::Integer, α, β) instead.", :BetaBinomial)
    BetaBinomial(Int(n), α, β)
end


text:
module Fields

export SpectralElementField, assemble, data

using CUDA
using Printf

using ClimateMachine.Mesh.Grids: DiscontinuousSpectralElementGrid
using ClimateMachine.MPIStateArrays: MPIStateArray
using ..Domains: AbstractDomain

#####
##### SpectralElementField
#####

struct SpectralElementField{E, D, G, A, V, R}
    elements::E
    domain::D
    grid::G
    data::A
    realdata::V
    realelems::R # Required to build realdata from data
end

"""
    SpectralElementField(domain::RectangularDomain, state::MPIStateArray, variable_index::Int)

Returns a Cartesian `view` into `state.realdata[:, variable_index, :]`,
assuming that `state.realdata` lives on `RectangularDomain`.

`SpectralElementField.elements` is a three-dimensional array of `RectangularElements`.
"""
function SpectralElementField(
    domain::AbstractDomain,
    grid::DiscontinuousSpectralElementGrid,
    state::MPIStateArray,
    variable_index::Int,
)

    data = view(state.data, :, variable_index, :)
    realdata = view(state.realdata, :, variable_index, :)
    realelems = state.realelems

    return SpectralElementField(domain, grid, data, realdata, realelems)
end

const SEF = SpectralElementField

Base.@propagate_inbounds Base.getindex(field::SEF, i, j, k) =
    field.elements[i, j, k]
Base.size(field::SEF) = size(field.elements)
Base.eltype(field::SEF) = eltype(field.realdata)

Base.maximum(f, field::SEF) = maximum([maximum(f, el) for el in field.elements])
Base.minimum(f, field::SEF) = minimum([minimum(f, el) for el in field.elements])

Base.maximum(field::SEF) = maximum([maximum(el) for el in field.elements])
Base.minimum(field::SEF) = minimum([minimum(el) for el in field.elements])

Base.show(io::IO, field::SEF{E}) where {E} =
    print(io, "SpectralElementField{$(E.name.wrapper)}")

#####
##### Domain-specific stuff
#####

include("rectangular_element.jl")
include("rectangular_spectral_element_fields.jl")

#####
##### CPU and GPU-friendly element assembly
#####

"""
    assemble(u::SpectralElementField)

Assemble `u.elements` into a single `element::eltype(u)`, averaging shared nodes.
"""
assemble(u::SpectralElementField) = assemble(u.elements)

"""
    assemble_data(u::SpectralElementField{<:CuArray})

Assemble the data in `u.elements` into a single `Array`, averaging shared nodes.
"""
function assemble(u::SpectralElementField{<:CuArray})

    domain = u.domain
    index = u.variable_index

    cpudata = Array(u.realdata)

    u_cpu = SpectralElementField(cpudata, domain)

    return assemble(u_cpu)
end

end # module

max_forks_repo_name: jm-c/CLIMA


text:
##############################################################################
# Validation
##############################################################################

testloss(R::PointResults) = R.testloss
trainloss(R::PointResults) = R.trainloss
thetaopt(R::PointResults) = R.theta
lambda(R::PointResults) = R.lambda

testloss(R::FoldResults,i) = testloss(R.results[i])
trainloss(R::FoldResults,i) = trainloss(R.results[i])
thetaopt(R::FoldResults,i) = thetaopt(R.results[i])
lambda(R::FoldResults,i) = lambda(R.results[i])

lambdapath(R::RegPathResults) = [r.lambda for r in R.results]
testlosspath(R::RegPathResults) = [r.testloss for r in R.results]
trainlosspath(R::RegPathResults) = [r.trainloss for r in R.results]

testloss(R::RegPathResults) = testloss(R.results[R.imin])
trainloss(R::RegPathResults) = trainloss(R.results[R.imin])
thetaopt(R::RegPathResults) = thetaopt(R.results[R.imin])
lambdaopt(R::RegPathResults) = lambda(R.results[R.imin])

function thetapath(R::RegPathResults)
    r = length(R.results)
    d = length(R.results[1].theta)
    T = zeros(d,r)
    for i=1:r
        T[:,i] = R.results[i].theta
    end
    return T'
end

testloss(M::Model,i) = testloss(M.D.results,i)
trainloss(M::Model,i) = trainloss(M.D.results,i)
thetaopt(M::Model,i) = thetaopt(M.D.results,i)
lambda(M::Model,i) = lambda(M.D.results,i)

"`testloss(M)` returns the average testing loss"
testloss(M::Model) = testloss(M.D.results)
"`trainloss(M)` returns the average training loss"
trainloss(M::Model) = trainloss(M.D.results)
"`thetaopt(M)` returns the optimal `theta` chosen during training"
thetaopt(M::Model) = thetaopt(M.D.results)
"`lambda(M)` returns the `lambda` used during training. If multiple `lambda` values were specified,
returns the optimal one"
lambda(M::Model) = lambda(M.D.results)

"`lambdaopt(M)` returns the optimal lambda computed during training of M"
lambdaopt(M::Model) = lambdaopt(M.D.results)

"`lambdapath(M)` returns the vector of lambdas used"
lambdapath(M::Model) = lambdapath(M.D.results)
"`testlosspath(M)` returns the vector of test losses at various lambdas"
testlosspath(M::Model) = testlosspath(M.D.results)
"`trainlosspath(M)` returns the vector of train losses at various lambdas"
trainlosspath(M::Model) = trainlosspath(M.D.results)
"`thetapath(M)` returns the vector of optimal thetas at various lambdas"
thetapath(M::Model) = thetapath(M.D.results)
##############################################################################


text:
#=
A mask is a BitArray

Masking is boolean.
Not a byte, so no degree of masking.

A mask is usually associated with an image.

It can represent:
- selection, something a user specifies
- transparency: is the corresponding image pixel visible,
- validity: is the corresponding image pixel part of a non-rectangular shape.
=#

include("boundingBox.jl")

function nand(x, y)
    return ~(x & y)
end

#=
Boolean subtraction
(clamped subtraction on 0,1)

Truth table

    0 1
0   0 0
1   1 0
=#
function subtractBool(x, y)
    if x
        return ~y
    else
        return false
    end
end

#=
Subtract masks

Case: not intersect
    all false
Case: intersect but second mask not contained
    the region outside the intersection but in mask1 or mask2
Case: mask2 contained by mask1
    the region in mask1 but not in mask2
=#

function subtractMasks(mask1, mask2)
    # FAIL, changes to Int64
    # result = .-(mask1, mask2)

    result = subtractBool.(mask1, mask2)
    @debug "subtractMasks" result
    @assert typeof(result) <: BitArray
    return result
end

#=
Create a new  mask which is a frisket around a given mask.

The frisket may not be a donut if the given mask
is closer than frisketDepth to the edge of the array.

The frisket is larger than the given mask unless the given mask is entire,
in which case the frisket is empty.
=#
function frisketMaskAroundMask(image, mask, frisketDepth)::BitArray

    # get bounding box of the given mask
    boundingBox = BoundingBox(mask)

    # Method of BoundingBox
    newBoundingBox = expandOrShrink(boundingBox, frisketDepth)

    # clamp expanded box to image
    constrainedBoundingBox = clampBoundingBoxToImage(newBoundingBox, image)

    # BitArray from bounding box, similar to mask
    moreMask = bitMaskFromInShapeOf(constrainedBoundingBox, mask)

    # Punch out (nand) the original mask (the selection in the target) to make a frisket (mask with a hole.)
    result = subtractMasks(moreMask, mask)

    @debug "frisket" result typeof(result)
    @assert typeof(result) <: BitArray
    return result
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
include("rrc_types.jl")

function getFileContent(fName::String)
  return unsafe_string(ccall(dlsym(rrlib, :getFileContent), cdecl, Ptr{UInt8}, (Ptr{UInt8},), fName))
end

function createText(text::String)
  return unsafe_string(ccall(dlsym(rrlib, :createText), cdecl, Ptr{UInt8}, (Ptr{UInt8},), text))
end

function createTextMemory(count::Int64)
  return unsafe_string(ccall(dlsym(rrlib, :createTextMemory), cdecl, Ptr{UInt8}, (Int64,), count))
end

function createRRList()
  return ccall(dlsym(rrlib, :createRRList), cdecl, Ptr{RRList}, ())
end

function freeRRList(theList)
  ccall(dlsym(rrlib, :freeRRList), cdecl, Nothing, (Ptr{RRList},), theList)
end

function getListLength(myList)
  return ccall(dlsym(rrlib, :getListLength), cdecl, Int64, (Ptr{RRList},), myList)
end

function createIntegerItem(value::Int64)
  return ccall(dlsym(rrlib, :createIntegerItem ), cdecl, Ptr{RRListItem}, (Int64,), value)
end

function createDoubleItem(value::Float64)
  return ccall(dlsym(rrlib, :createDoubleItem  ), cdecl, Ptr{RRListItem}, (Float64,), value)
end

function createStringItem(value::String)
  return ccall(dlsym(rrlib, :createStringItem), cdecl, Ptr{RRListItem}, (Ptr{UInt8},), value)
end

function createListItem(value::Ptr{RRList})
  return ccall(dlsym(rrlib, :createListItem), cdecl, Ptr{RRListItem}, (Ptr{RRList},), value)
end

function addItem(list::Ptr{RRList}, item::Ptr{RRListItem})
  return ccall(dlsym(rrlib, :addItem), cdecl, Int64, (Ptr{RRList}, Ptr{RRListItem}), list, item)
end

function getListItem(list::Ptr{RRList}, index::Int64)
  return ccall(dlsym(rrlib, :getListItem), cdecl, Ptr{RRListItem}, (Ptr{RRList}, Int64), list, index)
end

function isListItemInteger(item::Ptr{RRListItem})
  return ccall(dlsym(rrlib, :isListItemInteger ), cdecl, Bool, (Ptr{RRListItem},), item)
end

function isListItemDouble(item::Ptr{RRListItem})
  return ccall(dlsym(rrlib, :isListItemDouble ), cdecl, Bool, (Ptr{RRListItem},), item)
end

function isListItemString(item::Ptr{RRListItem})
  return ccall(dlsym(rrlib, :isListItemString), cdecl, Bool, (Ptr{RRListItem},), item)
end

function isListItemList(item::Ptr{RRListItem})
  return ccall(dlsym(rrlib, :isListItemList), cdecl, Bool, (Ptr{RRListItem},), item)
end

function isListItem(item::Ptr{RRListItem}, itemType)
  return ccall(dlsym(rrlib, :isListItem), cdecl, Bool, (Ptr{RRListItem}, Ptr{Nothing}), item, itemType)
end

function getIntegerListItem(item::Ptr{RRListItem})
  value = Array{Int64}(undef,1)
  status = ccall(dlsym(rrlib, :getIntegerListItem), cdecl, Bool, (Ptr{RRListItem}, Ptr{Int64}), item, value)
  if status == false
    error(getLastError())
  end
  return value[1]
end

function getDoubleListItem(item::Ptr{RRListItem})
  value = Array{Float64}(undef,1)
  status = ccall(dlsym(rrlib, :getDoubleListItem ), cdecl, Bool, (Ptr{RRListItem}, Ptr{Float64}), item, value)
  if status == false
    error(getLastError())
  end
  return value[1]
end

function getStringListItem(item::Ptr{RRListItem})
  return bytestring(ccall(dlsym(rrlib, :getStringListItem), cdecl, Ptr{UInt8}, (Ptr{RRListItem},), item))
end

function getList(item::Ptr{RRListItem})
  return ccall(dlsym(rrlib, :getList), cdecl, Ptr{RRList}, (Ptr{RRListItem},), item)
end

function freeRRCData(handle::Ptr{RRCData})
  status = ccall(dlsym(rrlib, :freeRRCData), cdecl, Bool, (Ptr{RRCData},), handle)
  if status == false
    (error(getLastError()))
  end
end

function freeText(text::String)
  status = ccall(dlsym(rrlib, :freeText), cdecl, Bool, (Ptr{UInt8},), text)
  if status == false
    (error(getLastError()))
  end
end

function freeStringArray(sl::Ptr{RRStringArray})
  status =  ccall(dlsym(rrlib, :freeStringArray), cdecl, Bool, (Ptr{RRStringArray},), sl)
  if status == false
    (error(getLastError()))
  end
end

function freeVector(vector::Ptr{RRVector})
  status = ccall(dlsym(rrlib, :freeVector), cdecl, Bool, (Ptr{RRVector},), vector)
  if status == false
    (error(getLastError()))
  end
end

function freeMatrix(matrix::Ptr{RRDoubleMatrix})
  status = ccall(dlsym(rrlib, :freeMatrix), cdecl, Bool, (Ptr{RRDoubleMatrix},), matrix)
  if status == false
    (error(getLastError()))
  end
end

function getVectorLength(vector::Ptr{RRVector})
  return ccall(dlsym(rrlib, :getVectorLength), cdecl, Int64, (Ptr{RRVector},), vector)
end

function createVector(size::Int64)
  return ccall(dlsym(rrlib, :createVector), cdecl, Ptr{RRVector}, (Int64,), size)
end

function getVectorElement(vector::Ptr{RRVector}, index::Int64)
  value = Array{Float64}(undef,1)
  status = ccall(dlsym(rrlib, :getVectorElement), cdecl, Bool, (Ptr{RRVector}, Int64, Float64), vector, index, value)
  if status == false
    error(getLastError())
  end
  return value[1]
end

function setVectorElement(vector::Ptr{RRVector}, index::Int64, value::Float64)
  status = ccall(dlsym(rrlib, :setVectorElement), cdecl, Bool, (Ptr{RRVector}, Int64, Float64), vector, index, value)
  if status == false
    error(getLastError())
  end
end

function createRRMatrix(r::Int64, c::Int64)
  return ccall(dlsym(rrlib, :createRRMatrix), cdecl, Ptr{RRDoubleMatrix}, (Int64, Int64), r, c)
end

function getMatrixNumRows(m::Ptr{RRDoubleMatrix})
  return ccall(dlsym(rrlib, :getMatrixNumRows), cdecl, Int64, (Ptr{RRDoubleMatrix},), m)
end

function getMatrixNumCols(m::Ptr{RRDoubleMatrix})
  return ccall(dlsym(rrlib, :getMatrixNumCols), cdecl, Int64, (Ptr{RRDoubleMatrix},), m)
end

function getMatrixElement(m::Ptr{RRDoubleMatrix}, r::Int64, c::Int64)
  value = Array{Float64}(undef,1)
  status = ccall(dlsym(rrlib, :getMatrixElement), cdecl, Bool, (Ptr{RRDoubleMatrix}, Int64, Int64, Ptr{Cdouble}), m, r, c, value)
  if status == false
    error(getLastError())
  end
  return value[1]
end

function setMatrixElement(m::Ptr{RRDoubleMatrix}, r::Int64, c::Int64, value::Float64)
  status = ccall(dlsym(rrlib, :setMatrixElement), cdecl, Bool, (Ptr{RRDoubleMatrix}, Int64, Int64, Float64), m, r, c, value)
  if status == false
    error(getLastError())
  end
end

function getComplexMatrixElement(m::Ptr{RRComplexMatrix}, r::Int64, c::Int64)
  value = Ptr{RRComplex}
  status = ccall(dlsym(rrlib, :getComplexMatrixElement), cdecl, Bool, (Ptr{RRComplexMatrix}, Int64, Int64, Ptr{RRComplex}), m, r, c, value)
  if status == false
    error(getLastError())
  end
  return value
end

function setComplexMatrixElement(m::Ptr{RRComplexMatrix}, r::Int64, c::Int64, value::Float64)
  status = ccall(dlsym(rrlib, :setComplexMatrixElement), cdecl, Bool, (Ptr{RRComplexMatrix}, Int64, Int64, Ptr{RRComplex}), m, r, c, value)
  if status == false
    error(getLastError())
  end
end

function getRRDataNumRows(rrData::Ptr{RRCData})
  return ccall(dlsym(rrlib, :getRRDataNumRows), cdecl, Int64, (Ptr{RRCData},), rrData)
end

function getRRDataNumCols(rrData::Ptr{RRCData})
  return ccall(dlsym(rrlib, :getRRDataNumCols), cdecl, Int64, (Ptr{RRCData},), rrData)
end

function getRRCDataElement(rrData::Ptr{RRCData}, r::Int64, c::Int64)
  value = Array{Float64}(undef,1)
  status = ccall(dlsym(rrlib, :getRRCDataElement), cdecl, Bool, (Ptr{RRCData}, Int64, Int64, Ptr{Float64}), rrData, r, c, value)
  if status == false
    error(getLastError())
  end
  return value[1]
end

function getRRDataColumnLabel(rrData::Ptr{RRCData}, column::Int64)
  return unsafe_string(ccall(dlsym(rrlib, :getRRDataColumnLabel), cdecl, Ptr{UInt8}, (Ptr{RRCData}, Int64), rrData, column))
end

function vectorToString(vecHandle::Ptr{RRVector})
  return unsafe_string(ccall(dlsym(rrlib, :vectorToString), cdecl, Ptr{UInt8}, (Ptr{RRVector},), vecHandle))
end

function complexVectorToString(vecHandle::Ptr{RRComplexVector})
  return unsafe_string(ccall(dlsym(rrlib, :complexVectorToString), cdecl, Ptr{UInt8}, (Ptr{RRComplexVector},), vecHandle))
end

function rrCDataToString(rrData::Ptr{RRCData})
  return unsafe_string(ccall(dlsym(rrlib, :rrCDataToString), cdecl, Ptr{UInt8}, (Ptr{RRCData},), rrData))
end

function matrixToString(matrixHandle::Ptr{RRDoubleMatrix})
  return unsafe_string(ccall(dlsym(rrlib, :matrixToString), cdecl, Ptr{UInt8}, (Ptr{RRDoubleMatrix},), matrixHandle))
end

function complexMatrixToString(matrixHandle::Ptr{RRComplexMatrix})
  return unsafe_string(ccall(dlsym(rrlib, :complexMatrixToString), cdecl, Ptr{UInt8}, (Ptr{RRComplexMatrix},), matrixHandle))
end

function getNumberOfStringElements(list::Ptr{RRStringArray})
  return ccall(dlsym(rrlib, :getNumberOfStringElements), cdecl, Int64, (Ptr{RRStringArray},), list)
end

function getStringElement(list::Ptr{RRStringArray}, index::Int64)
  return unsafe_string(ccall(dlsym(rrlib, :getStringElement), cdecl, Ptr{UInt8}, (Ptr{RRStringArray}, Int64), list, index))
end

function stringArrayToString(list::Ptr{RRStringArray})
  return unsafe_string(ccall(dlsym(rrlib, :stringArrayToString), cdecl, Ptr{UInt8}, (Ptr{RRStringArray},), list))
end

function listToString(list::Ptr{RRList})
  return unsafe_string(ccall(dlsym(rrlib, :listToString), cdecl, Ptr{UInt8}, (Ptr{RRList},), list))
end

function writeRRData(rr::Ptr{Nothing}, fileNameAndPath::String)
  status = ccall(dlsym(rrlib, :writeRRData), cdecl, Bool, (Ptr{Nothing}, Ptr{UInt8}), rr, fileNameAndPath)
  if status == false
    error(getLastError())
  end
end

function compileSource(rr::Ptr{Nothing}, fName::String)
  status = ccall(dlsym(rrlib, :compileSource), cdecl, Bool, (Ptr{Nothing}, Ptr{UInt8}), rr, fName)
  if status == false
    error(getLastError())
  end
end


text:
export read_resp_from_sacpz

# Set the response to acceleration
const flag = 2

function read_resp_from_sacpz(input_sacpz::String, sampling_rate::Float64, N::Int)
    read_zeros = false
    read_poles = false
    zeros = Array{Complex{Float64}}(undef, 1)
    poles = Array{Complex{Float64}}(undef, 1)
    zeros_cnt, poles_cnt = 1, 1
    constant = .0

    lines = readlines(input_sacpz)
    for line in lines
        contents = split(line, " ")
        if contents[1] == "CONSTANT"
            constant = parse(Float64, contents[2])
        elseif contents[1] == "ZEROS"
            zeros_size = parse(Int64, contents[2])
            zeros = Base.zeros(Complex{Float64}, zeros_size)
            read_zeros = true
            read_poles = false
        elseif contents[1] == "POLES"
            poles_size = parse(Int64, contents[2])
            poles = Base.zeros(Complex{Float64}, poles_size)
            read_zeros = false
            read_poles = true
        else
            # Seems that some SACPZ file will let Julia pase as follows:
            # "<real>", "", "<imag>"
            num = parse(Float64, contents[1]) + parse(Float64, contents[3]) * im
            if read_zeros == true
                zeros[zeros_cnt] = num
                zeros_cnt += 1
            elseif read_poles == true
                poles[poles_cnt] = num
                poles_cnt += 1
            end
        end
    end

    return create_resp(poles, zeros, constant, sampling_rate, N)
end

function create_resp(poles::Vector{Complex{Float64}},
                     zeros::Vector{Complex{Float64}},
                     constant::Float64,
                     sampling_rate::Float64,
                     N::Int)
    freqs = Vector{Float64}(undef, N)
    delta = 1. / sampling_rate
    total_duration = delta * N
    for i = 1:(div(N, 2) + 1)
        freqs[i] = i / total_duration
    end
    for i = (div(N, 2) + 1):N
        freqs[i] = -(N - i) / total_duration
    end

    instrument_response = Vector{Complex{Float64}}(undef, N)
    for i = 1:N
        numerator = 1. + 0. * im
        denominator = 1. + 0. * im
        omega = 2 * π * freqs[i]
        i_omega = 0. + omega * im

        for j = 1:(size(zeros, 1) - flag)
            numerator *= (i_omega - zeros[j])
        end
        for j = 1:size(poles, 1)
            denominator *= (i_omega - poles[j])
        end

        instrument_response[i] = constant * numerator / denominator
    end

    return instrument_response
end


text:
####
#### Subroutines for GLE
####

"""
    add_sub_marker!(f, m)

Internal function to add an appropriate subroutine to the GLE script to deal with markers that
must have a different color than the line they are associated with. For instance if you want a
blue line with red markers, you need to define a specfici subroutine for red-markers otherwise
both line and markers are going to be of the same color.
Note: these subroutines all start with the name `mk_....` (see [`str`](@ref)).
"""
function add_sub_marker!(f::Figure, m::MarkerStyle)
    str(m) ∈ keys(f.subroutines) && return nothing
    f.subroutines[str(m)] = """
        sub _$(str(m)) size mdata
        	gsave
            set color $(col2str(m.color))
            marker $(m.marker) size
        	grestore
        end sub
        define marker $(str(m)) _$(str(m))
        """
    return nothing
end

const GLE_DRAW_SUB = Dict{String,String}()

###############################################################
####
#### Boxplot subroutine
#=
NOTE: the outliers are drawn separately (see the boxplot function)

 ===  ----  y: whigh            x: p-wwidth/2 <> p+wwidth/2
  |
 +=+  ----  y: q75              x: p-bwidth/2 <> p+bwidth/2
 | |
 | |
 +=+  ----  y: q25
  |
 ===  ----  y: wlow

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

const boxplot_box_lstyle = """set lstyle blstyle lwidth blwidth color blcolor\$"""
const boxplot_med_lstyle = """set lstyle medlstyle lwidth medlwidth color medcolor\$"""
const boxplot_mean_mstyle = """
    gsave
    \tset color mmcol\$
    \tmarker mmarker\$ mmsize
    \tgrestore"""
const boxplot_core_vertical = """
    ! -----------------------------------------
    \t! SET LINE STYLE THEN DRAW BOX AND WHISKERS
    \t! -----------------------------------------
    \tgsave
    \tset cap round
    \t$boxplot_box_lstyle
    \t! LOWER WHISKER
    \tamove xg(p-wwidth/2) yg(wlow)
    \taline xg(p+wwidth/2) yg(wlow)
    \t! CONNECTION LOWER WHISKER-BOX
    \tamove xg(p) yg(wlow)
    \t! BOX
    \taline xg(p) yg(q25)
    \tamove xg(p-bwidth/2) yg(q25)
    \tbox xg(p+bwidth/2)-xg(p-bwidth/2) yg(q75)-yg(q25)
    \t! CONNECTION BOX-UPPER WHISKER
    \tamove xg(p) yg(q75)
    \taline xg(p) yg(whigh)
    \t! UPPER WHISKER
    \tamove xg(p-wwidth/2) yg(whigh)
    \taline xg(p+wwidth/2) yg(whigh)
    \tgrestore
    \t! ------------------------------------
    \t! SET LINE STYLE THEN DRAW MEDIAN LINE
    \t! ------------------------------------
    \tgsave
    \t$boxplot_med_lstyle
    \tamove xg(p-bwidth/2) yg(q50)
    \taline xg(p+bwidth/2) yg(q50)
    \tgrestore
    \t! --------------------------------------
    \t! SET MARKER STYLE THEN DRAW MEAN MARKER
    \t! --------------------------------------
    \tif (mshow > 0) then
    \t    amove xg(p) yg(mean)
    \t$boxplot_mean_mstyle
    \tend if"""
const boxplot_core_horizontal = """
    ! -----------------------------------------
    \t! SET LINE STYLE THEN DRAW BOX AND WHISKERS
    \t! -----------------------------------------
    \tgsave
    \tset cap round
    \t$boxplot_box_lstyle
    \tamove xg(wlow) yg(p-wwidth/2)
    \taline xg(wlow) yg(p+wwidth/2)
    \tamove xg(wlow) yg(p)
    \taline xg(q25) yg(p)
    \tamove xg(q25) yg(p-bwidth/2)
    \tbox xg(q75)-xg(q25) yg(p+bwidth/2)-yg(p-bwidth/2)
    \tamove xg(q75) yg(p)
    \taline xg(whigh) yg(p)
    \tamove xg(whigh) yg(p-wwidth/2)
    \taline xg(whigh) yg(p+wwidth/2)
    \tgrestore
    \t! ------------------------------------
    \t! SET LINE STYLE THEN DRAW MEDIAN LINE
    \t! ------------------------------------
    \tgsave
    \t$boxplot_med_lstyle
    \tamove xg(q50) yg(p-bwidth/2)
    \taline xg(q50) yg(p+bwidth/2)
    \tgrestore
    \t! --------------------------------------
    \t! SET MARKER STYLE THEN DRAW MEAN MARKER
    \t! --------------------------------------
    \tif (mshow > 0) then
    \t    amove xg(mean) yg(p)
    \t    $boxplot_mean_mstyle
    \tend if"""
const boxplot_args = ("p wlow q25 q50 q75 whigh mean ", # NOTE don't forget space at the end
                      "bwidth wwidth blstyle blwidth blcolor\$ ",
                      "medlstyle medlwidth medcolor\$ ",
                      "mshow mmarker\$ mmsize mmcol\$ ")
GLE_DRAW_SUB["bp_vert"] = """
    sub bp_vert $(prod(boxplot_args))
        $boxplot_core_vertical
    end sub
    """
GLE_DRAW_SUB["bp_horiz"] = """
    sub bp_horiz $(prod(boxplot_args))
        $boxplot_core_horizontal
    end sub
    """

###############################################################
####
#### Heatmap Subroutine
####
###############################################################

"""
    add_sub_marker!(f, m)

Internal function to add a heatmap subroutine to the GLE script.
"""
function add_sub_heatmap!(f::Figure, hm::Heatmap, hashid::UInt)
    #=
    if zij <= 1 then
        col$ = firstcolor
    else if zij <=2 then
        col$ = secondcolor
    ...
    else
        col$ = missingcolor
    end
    =#
    ifpart = """
        if zij <= 1 then
        \t\t    cij\$ = \"$(col2str(hm.cmap[1]))\""""
    for k ∈ 2:length(hm.cmap)-1
        ifpart *= """
            \n\t\telse if zij <= $k then
            cij\$ = \"$(col2str(hm.cmap[k]))\""""
    end
    ifpart *= """
        \n\t\telse
            cij\$ = \"$(col2str(hm.cmiss))\"
        end if"""

    boxpart = ifelse(hm.transpose, """
        amove xg((i-1)*bh) yg(1-j*bw)
        \t\tbox xg(bh)-xg(0) yg(bw)-yg(0) nobox fill cij\$""", """
        amove xg((j-1)*bw) yg(1-i*bh)
        \t\tbox xg(bw)-xg(0) yg(bh)-yg(0) nobox fill cij\$""")

    f.subroutines["hm_$hashid"] = """
        sub hm_$hashid j ds\$ bw bh
            local zij = 0
            local cij = \"\"
            for i = 1 to ndata(ds\$)
                zij = datayvalue(ds\$,i)
                $ifpart
                $boxpart
            next i
        end sub
        """
    return nothing
end

###############################################################
####
#### Palette a vector of colors
####
###############################################################

function add_sub_palette!(f::Figure, vc::Vector{<:Color})
    nc   = length(vc)-1
    incr = 1/(length(vc)-1)
    bot = vc[1]
    top = vc[2]
    core = """
        local r = 0
        \tlocal g = 0
        \tlocal b = 0
        \tif (z <= $incr) then
        \t    r = $(round3d(bot.r))*(1-z*$nc)+$(round3d(top.r))*z*$nc
        \t    g = $(round3d(bot.g))*(1-z*$nc)+$(round3d(top.g))*z*$nc
        \t    b = $(round3d(bot.b))*(1-z*$nc)+$(round3d(top.b))*z*$nc
        """
    for i ∈ 2:length(vc)-1
        bot = vc[i] # bottom color
        top = vc[i+1]   # top color
        core *= """
            \telse if ($(incr*(i-1)) < z) and (z <= $(incr*i)) then
            \t    r = $(round3d(bot.r))*(1-(z-$(incr*(i-1)))*$nc)+$(round3d(top.r))*(z-$(incr*(i-1)))*$nc
            \t    g = $(round3d(bot.g))*(1-(z-$(incr*(i-1)))*$nc)+$(round3d(top.g))*(z-$(incr*(i-1)))*$nc
            \t    b = $(round3d(bot.b))*(1-(z-$(incr*(i-1)))*$nc)+$(round3d(top.b))*(z-$(incr*(i-1)))*$nc
            """
    end
    core *= """\tend if"""

    pname = "cmap_$(hash(vc))"
    pname ∈ keys(f.subroutines) && return nothing
    f.subroutines[pname] = """
        sub $pname z
            $core
            return rgb(r,g,b)
        end sub
        """
    return nothing
end

###############################################################
####
#### Plot a 3D line
####
###############################################################

function add_sub_plot3!(f::Figure)
    "plot3" ∈ keys(f.subroutines) && return nothing
    xs = "(xo-xmin)/xspan"
    ys = "(yo-ymin)/yspan"
    f.subroutines["plot3"] = """
        sub plot3 data\$ xmin xspan ymin yspan showline showmarker marker\$ mscale
            xo = 0
            yo = 0
            zo = 0
            io = 0
            fopen data\$ file read
            until feof(file)
                fread file x y z
                xo = x
                yo = y
                zo = z
                io = io+1
                if (io<=1) then
                    amove xg3d($xs,$ys,zo) yg3d($xs,$ys,zo)
                    if (showmarker > 0) then
                        marker marker\$ mscale
                    end if
                else
                    if (showline >= 1) then
                        aline xg3d($xs,$ys,zo) yg3d($xs,$ys,zo)
                    else
                        amove xg3d($xs,$ys,zo) yg3d($xs,$ys,zo)
                    end if
                    if (showmarker > 0) then
                        marker marker\$ mscale
                    end if
                end if
            next
            fclose file
        end sub
        """
    return nothing
end


text:
module TaskHeritableStorage # end

"""
    fork_task_local_storage(key, value) -> (key,value) or nothing

Whenever a new `Task` is created, all `(key,value)` pairs in the parent Task's
`task_local_storage` are passed to this function, which gives a chance to _fork_ those pairs
into the new child `Task`'s task local storage.

Whatever new `(key, value)` pair is returned will be inserted into the new Task's
`task_local_storage`. If `nothing` is returned, the child Task's storage will not be
modified. This is the default behavior.

This function is called from inside the newly created `Task` before anything else happens,
so within this function `current_task()` will return the new Task.
"""
function fork_task_local_storage(key, value)
    nothing
end

# ----------------------
# Replace the definition of Core._Task() to clone task_heritable_storage on construction.

function _fork_task_local_storage_from_parent(parent_tls)
    tls = task_local_storage()
    # Fork all (k,v) pairs from the parent's TLS into this task's TLS
    for (k, v) in parent_tls
        out = fork_task_local_storage(k,v)
        if out !== nothing
            (newkey, newval) = out
            tls[newkey] = newval
        end
    end
end

# DANGEROUS: Override the Core Task() constructor to copy the task heritable storage
@warn "About to replace Core._Task() definition, to enable TaskHeritableStorage. The following warning is expected:"
#  EEP: type-piracy!  (In the future, this would be implemented in julia itself) ☠
function __Task(@nospecialize(f), reserved_stack::Int, completion_future)
    return ccall(:jl_new_task, Ref{Task}, (Any, Any, Int), f, completion_future, reserved_stack)
end
function Core._Task(@nospecialize(f), reserved_stack::Int, completion_future)
    if current_task().storage !== nothing
        let parent_tls = task_local_storage()
            wrapped = () -> begin
                _fork_task_local_storage_from_parent(parent_tls);
                f();
            end
            return __Task(wrapped, reserved_stack, completion_future)
        end
    else
        return __Task(f, reserved_stack, completion_future)
    end
end

end