text:
eye(T::Type, n::Integer) = Matrix{T}(LinearAlgebra.I(n))
eye(n::Integer) = eye(Float64, n)

"""
weighted sampling: takes (unnormalized) log probs and returns index
"""
function wsample_logprob(logProbs::AbstractVector{<:Real})
  log_p_max = maximum(logProbs)
  p = exp.(logProbs .- log_p_max)
  return Distributions.wsample(p)
end

"""
log 1 minus. For example, `log1m(.3) == log(1 - .3) == log1p(-.3) == log(.7)`.
"""
log1m(x::T) where {T <: Real} = log1p(-x)

"""
Convert a vector of indices to a one-hot matrix, given also the number of
categories.
"""
function to_onehot(indvec::Vector{Int}, categories::Int)
  N = length(indvec)
  onehot_matrix = zeros(Int, N, categories)
  for n in 1:N
    onehot_matrix[n, indvec[n]] = 1
  end
  return onehot_matrix
end

"""
log pdf of gaussian mixture model, for a vector of x of size N and (m: mixture
locations, s: mixture scales, w: mixture weights) each vectors of size K
"""
function gmmlogpdf(m, s, w, x; dims)
  return logsumexp(normlogpdf.(m, s, x) .+ log.(w), dims=dims)
end

"""
log pdf of gaussian mixture model, for a scalar x N and (m: mixture locations,
s: mixture scales, w: mixture weights) each vectors of size K
"""
function gmmlogpdf(m, s, w, x::Real)
  return logsumexp(normlogpdf.(m, s, x) .+ log.(w))
end

"""
Partitions a list `xs` into `(good, bad)` according
to a condition.
"""
function partition(condition::Function, xs::Vector{T}) where T
  good = T[]
  bad = T[]

  for x in xs
    if condition(x)
      append!(good, [x])
    else
      append!(bad, [x])
    end
  end

  return good, bad
end

function quantiles(X, q; dims, drop=false)
  Q = mapslices(x -> quantile(x, q), X, dims=dims)
  out = drop ? dropdims(Q, dims=dims) : Q
  return out
end

function update_mean!(current_mean, x, iter)
  current_mean .+= (x - current_mean) / iter
end

function update_cov!(current_cov, current_mean, x, iter)
  d = x - current_mean
  current_cov .= current_cov * (iter - 1)/iter + (d*d') * (iter - 1)/iter^2
  current_cov .= Matrix(LinearAlgebra.Symmetric(current_cov))
end

"""
x: real vector of dim K - 1
return: simplex of dim K
"""
function simplex_invtransform(x::AbstractVector{<:Real})
  K = length(x) + 1
  z = logistic.(x - log.(K .- collect(1:K-1)))
  p = zeros(K)
  for k in 1:K-1
    p[k] = (1 - sum(p[1:k-1])) * z[k]
  end
  p[K] = 1 - sum(p[1:K-1])
  return p  # simplex
end

"""
Given simplex p (length K), return real vector of length K-1.
"""
function simplex_transform(p::AbstractVector{<:Real})
  K = length(p)
  ks = collect(1:K-1)
  z = [p[k] / (1 - sum(p[1:k-1])) for k in ks]
  return logit.(z) + log.(K .- ks)
end

"""
x: real vector of dim K - 1
return: log abs value of determinant of x
"""
function simplex_logabsdet(x::AbstractVector{<:Real})
  K = length(x) + 1
  ks = collect(1:K-1)
  z = logistic.(x - log.(K .- ks))
  p = simplex_invtransform(x)
  return sum(log1p.(-z) + log.(p[ks]))
end

# DIC.
deviance(loglikelihood::Real) = -2 * loglikelihood
function dic(loglikelihood::AbstractVector{<:Real})
  D = deviance.(loglikelihood)
  return mean(D) + 0.5 * var(D)
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
import Pkg; Pkg.activate(".")
using SimpleSDMLayers
using Plots
using StatsPlots
using Statistics
using StatsBase
using GBIF

## World-scale example

# Get world temperature
temperature = worldclim(1)
# Verify size
size(temperature)
# Map temperature
temperature_plot = plot(temperature, xguide = "Longitude", yguide = "Latitude", colorbar_title = "Average Temperature (°C)")
# Save result
savefig(temperature_plot, "fig/temperature.png")

## Common manipulations

# Clip to Europe
temperature_europe = temperature[left=-11.2, right=30.6, bottom=29.2, top=71.0];
# Verify size again
size(temperature_europe)
# Map result
temperature_europlot = plot(temperature_europe, frame = :box, ticks = false, colorbar = false)
# Save
savefig(temperature_europlot, "fig/temperature_europe.png")

# Coarsen example
temperature_europe_coarse = coarsen(temperature_europe, Statistics.mean, (4, 4))
plot(temperature_europe_coarse, frame = :box, ticks = false, colorbar = false)
# Save
savefig("fig/coarsen.png")

# Sliding window example
averaged = slidingwindow(temperature_europe, Statistics.mean, 100.0)
plot(averaged, frame = :box, ticks = false, colorbar = false)
# Save
savefig("fig/slidingwindow.png")

## GBIF integration

# Get taxon & occurrences
kingfisher = GBIF.taxon("Megaceryle alcyon", strict=true)
kf_occurrences = occurrences(kingfisher)

# Get some more occurrences
for i in 1:9
  occurrences!(kf_occurrences)
end

# Filter for correct coordinates
# filter!(GBIF.have_ok_coordinates, kf_occurrences)
# @info kf_occurrences

# Clip layer to occurrences
temperature_clip = clip(temperature, kf_occurrences)

# Plot occurrences
contour(temperature_clip, fill=true, colorbar_title = "Average temperature (°C)",
                  xguide = "Longitude", yguide = "Latitude")
scatter!(longitudes(kf_occurrences), latitudes(kf_occurrences),
         label = "Kingfisher occurrences", legend = :bottomleft,
         c=:white, msc=:orange)
# Save
savefig("fig/occurrences.png")

## BIOCLIM example

include("bioclim.jl")

# Get all worldclim variables
vars = worldclim(1:19)

# Get prediction for each variable
vars_predictions = [bioclim(v, kf_occurrences) for v in vars]

# Get minimum prediction per site
sdm_predictions = reduce(min, vars_predictions)

# Set value to nothing if prediction is zero
replace!(x -> isnothing(x) || iszero(x) ? nothing : x, sdm_predictions.grid)

# Filter predictions with threshold
# threshold = quantile(filter(!isnothing, sdm_predictions.grid), 0.05)
# replace!(x -> isnothing(x) || x <= threshold ? nothing : x, sdm_predictions.grid)

# Map background
plot(temperature_clip, c = :lightgrey,
     xguide = "Longitude", yguide = "Latitude")
# Map predictions
plot!(clip(sdm_predictions, kf_occurrences), c = :viridis,
      clim = (minimum(sdm_predictions), maximum(sdm_predictions)),
      colorbar_title = "Predicted suitability score")
# Map occurrences
# scatter!(longitudes(kf_occurrences), latitudes(kf_occurrences),
#          label = "Kingfisher occurrences", legend = :bottomleft,
#          c=:white, msc=:orange)
# Save result
savefig("fig/bioclim.png")

ext: jl


text:
function norm2_values(x)
    if x isa Dict
        vals = collect(values(x))
        return vals
    elseif x isa AbstractArray || x isa Tuple
        return flatten(map(norm2_values, x))
    elseif x isa Function
        return []
    elseif Yota.isstruct(x)
        sub_iters = []
        for fld in fieldnames(typeof(x))
            push!(sub_iters, norm2_values(getfield(x, fld)))
        end
        return flatten(sub_iters)
    else
        return x
    end
end

function norm2(x)
    Statistics.norm(collect(norm2_values(x)), 2)
end

"""
Walk fields of struct `s` applying callback function `cb`.
If `cb` returns nothing, recursively call `walkstruct()` on each field.
"""
function walkstruct(cb::Function, s)
    rec = cb(s)
    if rec == nothing
        for p_name in propertynames(s)
            prop = getproperty(s, p_name)
            walkstruct(cb, prop)
        end
    end
end


text:
"""
Arguments:

    * `f`: a row kernel
    * `tbl`: An `AbstractTable`
    * `arg_flds`: A mapping from numeric indices to names of argument columns

Returns:

    * `T::DataType`: the inferred type of applying `f` to an argument
    with signature the inner `eltype`s of the argument columns
    * `row_itr`: An iterator over the rows of the argument columns

"""
@noinline function _preprocess(f, tbl, arg_flds)
    # Extract columns
    row_itr = eachrow(tbl, arg_flds...)

    # Determine the type of the iterator (possibly after unwrapping Nullable
    # types)
    inner_eltypes = map(eltype, eltypes(tbl, arg_flds...))

    # See if inference knows the return type for rowwise kernel
    T = Core.Inference.return_type(f, (Tuple{inner_eltypes...}, ))

    # Branch here on non-concrete types on a semantically correct path.
    if !isleaftype(T)
        # TODO: Make this a real code path based on the type-inference
        # independent map implementation in Base.
        @printf("WARNING: Failed to type-infer expression: found %s", T)
        T = Any
    end
    return T, row_itr
end


text:
"""
   forwardrecurrence!(v, A, B, C, x)

evaluates the orthogonal polynomials at points `x`,
where `A`, `B`, and `C` are `AbstractVector`s containing the recurrence coefficients
as defined in DLMF,
overwriting `v` with the results.
"""
function forwardrecurrence!(v::AbstractVector{T}, A::AbstractVector, B::AbstractVector, C::AbstractVector, x, p0=one(T)) where T
    N = length(v)
    N == 0 && return v
    length(A)+1 ≥ N && length(B)+1 ≥ N && length(C)+1 ≥ N || throw(ArgumentError("A, B, C must contain at least $(N-1) entries"))
    p1 = convert(T, N == 1 ? p0 : muladd(A[1],x,B[1])*p0) # avoid accessing A[1]/B[1] if empty
    _forwardrecurrence!(v, A, B, C, x, convert(T, p0), p1)
end

Base.@propagate_inbounds _forwardrecurrence_next(n, A, B, C, x, p0, p1) = muladd(muladd(A[n],x,B[n]), p1, -C[n]*p0)
# special case for B[n] == 0
Base.@propagate_inbounds _forwardrecurrence_next(n, A, ::Zeros, C, x, p0, p1) = muladd(A[n]*x, p1, -C[n]*p0)
# special case for Chebyshev U
Base.@propagate_inbounds _forwardrecurrence_next(n, A::AbstractFill, ::Zeros, C::Ones, x, p0, p1) = muladd(getindex_value(A)*x, p1, -p0)

# this supports adaptivity: we can populate `v` for large `n`
function _forwardrecurrence!(v::AbstractVector, A::AbstractVector, B::AbstractVector, C::AbstractVector, x, p0, p1)
    N = length(v)
    N == 0 && return v
    v[1] = p0
    N == 1 && return v
    v[2] = p1
    @inbounds for n = 2:N-1
        p1,p0 = _forwardrecurrence_next(n, A, B, C, x, p0, p1),p1
        v[n+1] = p1
    end
    v
end

forwardrecurrence(N::Integer, A::AbstractVector, B::AbstractVector, C::AbstractVector, x) =
    forwardrecurrence!(Vector{promote_type(eltype(A),eltype(B),eltype(C),typeof(x))}(undef, N), A, B, C, x)

"""
clenshaw!(c, A, B, C, x)

evaluates the orthogonal polynomial expansion with coefficients `c` at points `x`,
where `A`, `B`, and `C` are `AbstractVector`s containing the recurrence coefficients
as defined in DLMF,
overwriting `x` with the results.
"""
clenshaw!(c::AbstractVector, A::AbstractVector, B::AbstractVector, C::AbstractVector, x::AbstractVector) =
    clenshaw!(c, A, B, C, x, Ones{eltype(x)}(length(x)), x)

"""
clenshaw!(c, A, B, C, x, ϕ₀, f)

evaluates the orthogonal polynomial expansion with coefficients `c` at points `x`,
where `A`, `B`, and `C` are `AbstractVector`s containing the recurrence coefficients
as defined in DLMF and ϕ₀ is the zeroth polynomial,
overwriting `f` with the results.
"""
function clenshaw!(c::AbstractVector, A::AbstractVector, B::AbstractVector, C::AbstractVector, x::AbstractVector, ϕ₀::AbstractVector, f::AbstractVector)
    f .= ϕ₀ .* clenshaw.(Ref(c), Ref(A), Ref(B), Ref(C), x)
end

Base.@propagate_inbounds _clenshaw_next(n, A, B, C, x, c, bn1, bn2) = muladd(muladd(A[n],x,B[n]), bn1, muladd(-C[n+1],bn2,c[n]))
Base.@propagate_inbounds _clenshaw_next(n, A, ::Zeros, C, x, c, bn1, bn2) = muladd(A[n]*x, bn1, muladd(-C[n+1],bn2,c[n]))
# Chebyshev U
Base.@propagate_inbounds _clenshaw_next(n, A::AbstractFill, ::Zeros, C::Ones, x, c, bn1, bn2) = muladd(getindex_value(A)*x, bn1, -bn2+c[n])

# allow special casing first arg, for ChebyshevT in OrthogonalPolynomialsQuasi
Base.@propagate_inbounds _clenshaw_first(A, B, C, x, c, bn1, bn2) = muladd(muladd(A[1],x,B[1]), bn1, muladd(-C[2],bn2,c[1]))

"""
    clenshaw(c, A, B, C, x)

evaluates the orthogonal polynomial expansion with coefficients `c` at points `x`,
where `A`, `B`, and `C` are `AbstractVector`s containing the recurrence coefficients
as defined in DLMF.
`x` may also be a single `Number`.
"""

function clenshaw(c::AbstractVector, A::AbstractVector, B::AbstractVector, C::AbstractVector, x::Number)
    N = length(c)
    T = promote_type(eltype(c),eltype(A),eltype(B),eltype(C),typeof(x))
    @boundscheck check_clenshaw_recurrences(N, A, B, C)
    N == 0 && return zero(T)
    @inbounds begin
        bn2 = zero(T)
        bn1 = convert(T,c[N])
        N == 1 && return bn1
        for n = N-1:-1:2
            bn1,bn2 = _clenshaw_next(n, A, B, C, x, c, bn1, bn2),bn1
        end
        bn1 = _clenshaw_first(A, B, C, x, c, bn1, bn2)
    end
    bn1
end

clenshaw(c::AbstractVector, A::AbstractVector, B::AbstractVector, C::AbstractVector, x::AbstractVector) =
    clenshaw!(c, A, B, C, copy(x))

###
# Chebyshev T special cases
###

"""
   clenshaw!(c, x)

evaluates the first-kind Chebyshev (T) expansion with coefficients `c` at points `x`,
overwriting `x` with the results.
"""
clenshaw!(c::AbstractVector, x::AbstractVector) = clenshaw!(c, x, x)

"""
   clenshaw!(c, x, f)

evaluates the first-kind Chebyshev (T) expansion with coefficients `c` at points `x`,
overwriting `f` with the results.
"""
clenshaw!(c::AbstractVector, x::AbstractVector, f::AbstractVector) = _clenshaw!(MemoryLayout(c), MemoryLayout(x), MemoryLayout(f), c, x, f)

function _clenshaw!(_, _, _, c::AbstractVector, x::AbstractVector, f::AbstractVector)
    f .= clenshaw.(Ref(c), x)
end

"""
    clenshaw(c, x)

evaluates the first-kind Chebyshev (T) expansion with coefficients `c` at  the points `x`.
`x` may also be a single `Number`.
"""
function clenshaw(c::AbstractVector, x::Number)
    N,T = length(c),promote_type(eltype(c),typeof(x))
    if N == 0
        return zero(T)
    elseif N == 1 # avoid issues with NaN x
        return first(c)*one(x)
    end

    y = 2x
    bk1,bk2 = zero(T),zero(T)
    @inbounds begin
        for k = N:-1:2
            bk1,bk2 = muladd(y,bk1,c[k]-bk2),bk1
        end
        muladd(x,bk1,c[1]-bk2)
    end
end

clenshaw(c::AbstractVector, x::AbstractVector) = clenshaw!(c, copy(x))


text:
"""
    GridData <: AbstractArray

Simulation data specific to a single grid.
"""
abstract type GridData{Y,X,R,T,N,I} <: AbstractArray{T,N} end

function (::Type{G})(d::GridData{Y,X,R,T,N}) where {G<:GridData,Y,X,R,T,N}
    args = init(d), mask(d), proc(d), opt(d), boundary(d), padval(d),
        source(d), dest(d), sourcestatus(d), deststatus(d)
    G{Y,X,R,T,N,map(typeof, args)...}(args...)
end
function ConstructionBase.constructorof(::Type{T}) where T<:GridData{Y,X,R} where {Y,X,R}
    T.name.wrapper{Y,X,R}
end

GridDataOrReps = Union{GridData, Vector{<:GridData}}

# Array interface
Base.size(d::GridData) = size(init(d))
Base.axes(d::GridData) = axes(init(d))
Base.eltype(d::GridData) = eltype(init(d))
Base.firstindex(d::GridData) = firstindex(init(d))
Base.lastindex(d::GridData) = lastindex(init(d))

# Getters
init(d::GridData) = d.init
mask(d::GridData) = d.mask
radius(d::GridData{<:Any,<:Any,R}) where R = R
proc(d::GridData) = d.proc
opt(d::GridData) = d.opt
boundary(d::GridData) = d.boundary
padval(d::GridData) = d.padval
source(d::GridData) = d.source
dest(d::GridData) = d.dest
sourcestatus(d::GridData) = d.sourcestatus
deststatus(d::GridData) = d.deststatus

gridsize(d::GridData) = size(init(d))
gridsize(A::AbstractArray) = size(A)
gridsize(nt::NamedTuple) = gridsize(first(nt))
gridsize(nt::NamedTuple{(),Tuple{}}) = 0, 0
gridaxes(d::GridData) = map(Base.OneTo, gridsize(d))
gridview(d::GridData) = sourceview(d)
sourceview(d::GridData) = view(parent(source(d)), map(a -> a .+ radius(d), gridaxes(d))...)
destview(d::GridData) = view(parent(dest(d)), map(a -> a .+ radius(d), gridaxes(d))...)
source_array_or_view(d::GridData) = source(d) isa OffsetArray ? sourceview(d) : source(d)
dest_array_or_view(d::GridData) = dest(d) isa OffsetArray ? destview(d) : dest(d)

function Base.copy!(grid::GridData{<:Any,<:Any,R}, A::AbstractArray) where R
    pad_axes = map(ax -> ax .+ R, axes(A))
    copyto!(parent(source(grid)), CartesianIndices(pad_axes), A, CartesianIndices(A))
    _updatestatus!(grid)
    return grid
end
function Base.copy!(A::AbstractArray, grid::GridData{<:Any,<:Any,R}) where R
    pad_axes = map(ax -> ax .+ R, axes(A))
    copyto!(A, CartesianIndices(A), parent(source(grid)), CartesianIndices(pad_axes))
    return grid
end
function Base.copy!(dst::GridData{<:Any,<:Any,RD}, src::GridData{<:Any,<:Any,RS}) where {RD,RS}
    dst_axes = map(s -> RD:s + RD, gridsize(dst))
    src_axes = map(s -> RS:s + RS, gridsize(src))
    copyto!(parent(source(dst)), CartesianIndices(dst_axes),
            parent(source(src)), CartesianIndices(src_axes)
    )
    return grid
end

@propagate_inbounds Base.getindex(d::GridData, I...) = getindex(source(d), I...)

# Swap source and dest arrays. Allways returns regular SimData.
_swapsource(d::Tuple) = map(_swapsource, d)
function _swapsource(grid::GridData)
    src = grid.source
    dst = grid.dest
    @set! grid.dest = src
    @set! grid.source = dst
    srcstatus = grid.sourcestatus
    dststatus = grid.deststatus
    @set! grid.deststatus = srcstatus
    return @set grid.sourcestatus = dststatus
end

#=
Initialise the block status array.
This tracks whether anything has to be done in an area of the main array.
=#
_updatestatus!(grid::GridData) = _updatestatus!(opt(grid), grid)
function _updatestatus!(opt::SparseOpt, grid)
    blocksize = 2 * radius(grid)
    src = parent(source(grid))
    for i in CartesianIndices(src)
        # Mark the status block if there is a non-zero value
        if !can_skip(opt, src[i])
            bi = _indtoblock.(Tuple(i), blocksize)
            @inbounds sourcestatus(grid)[bi...] = true
            @inbounds deststatus(grid)[bi...] = true
        end
    end
    return nothing
end
_updatestatus!(opt, grid) = nothing

#=
Find the maximum radius required by all rules
Add padding around the original init array, offset into the negative
So that the first real cell is still 1, 1
=#
function _addpadding(init::AbstractArray{T,N}, r, padval) where {T,N}
    h, w = size(init)
    paddedsize = h + 4r, w + 2r
    paddedindices = -r + 1:h + 3r, -r + 1:w + r
    sourceparent = fill(convert(eltype(init), padval), paddedsize)
    source = OffsetArray(sourceparent, paddedindices...)
    # Copy the init array to the middle section of the source array
    for j in 1:size(init, 2), i in 1:size(init, 1)
        @inbounds source[i, j] = init[i, j]
    end
    return source
end

function _build_status(opt::SparseOpt, source, r)
    hoodsize = 2r + 1
    blocksize = 2r
    nblocs = _indtoblock.(size(source), blocksize) .+ 1
    sourcestatus = zeros(Bool, nblocs)
    deststatus = zeros(Bool, nblocs)
    sourcestatus, deststatus
end
_build_status(opt::PerformanceOpt, init, r) = nothing, nothing

"""
    ReadableGridData <: GridData

    ReadableGridData(grid::GridData)
    ReadableGridData{Y,X,R}(init::AbstractArray, mask, opt, boundary, padval)

Simulation data and storage passed to rules for each timestep.

# Type parameters

- `Y`: number of rows
- `X`: number of columns
- `R`: grid padding radius
"""
struct ReadableGridData{
    Y,X,R,T,N,I<:AbstractArray{T,N},M,P<:Processor,Op<:PerformanceOpt,Ov,PV,S,D,SSt,DSt
} <: GridData{Y,X,R,T,N,I}
    init::I
    mask::M
    proc::P
    opt::Op
    boundary::Ov
    padval::PV
    source::S
    dest::D
    sourcestatus::SSt
    deststatus::DSt
end
function ReadableGridData{Y,X,R}(
    init::I, mask::M, proc::P, opt::Op, boundary::Ov, padval::PV, source::S,
    dest::D, sourcestatus::SSt, deststatus::DSt
) where {Y,X,R,I<:AbstractArray{T,N},M,P,Op,Ov,PV,S,D,SSt,DSt} where {T,N}
    ReadableGridData{Y,X,R,T,N,I,M,P,Op,Ov,PV,S,D,SSt,DSt}(
        init, mask, proc, opt, boundary, padval, source, dest, sourcestatus, deststatus
    )
end
# Generate simulation data to match a ruleset and init array.
@inline function ReadableGridData{X,Y,R}(
    init::AbstractArray, mask, proc, opt, boundary, padval
) where {Y,X,R}
    # We add one extra row and column of status blocks so
    # we dont have to worry about special casing the last block
    if R > 0
        source = _addpadding(init, R, padval)
        dest = _addpadding(init, R, padval)
    else
        if opt isa SparseOpt
            opt = NoOpt()
        end
        source = deepcopy(init)
        dest = deepcopy(init)
    end
    sourcestatus, deststatus = _build_status(opt, source, R)

    grid = ReadableGridData{X,Y,R}(
        init, mask, proc, opt, boundary, padval, source, dest, sourcestatus, deststatus
    )
    _updatestatus!(grid)
    return grid
end

Base.parent(d::ReadableGridData) = parent(source(d))

"""
    WritableGridData <: GridData

    WritableGridData(grid::GridData)

Passed to rules as write grids, and can be written to directly as an array,
or preferably using `add!` etc. All writes handle updates to SparseOpt()
and writing to the correct source/dest array.

Reads are _always from the source array_, as rules must not be sequential between
cells. This means using e.g. `+=` is not supported, instead use `add!`.
"""
struct WritableGridData{
    Y,X,R,T,N,I<:AbstractArray{T,N},M,P<:Processor,Op<:PerformanceOpt,Ov,PV,S,D,SSt,DSt
} <: GridData{Y,X,R,T,N,I}
    init::I
    mask::M
    proc::P
    opt::Op
    boundary::Ov
    padval::PV
    source::S
    dest::D
    sourcestatus::SSt
    deststatus::DSt
end
function WritableGridData{Y,X,R}(
    init::I, mask::M, proc::P, opt::Op, boundary::Ov, padval::PV,
    source::S, dest::D, sourcestatus::SSt, deststatus::DSt
) where {Y,X,R,I<:AbstractArray{T,N},M,P,Op,Ov,PV,S,D,SSt,DSt} where {T,N}
    WritableGridData{Y,X,R,T,N,I,M,P,Op,Ov,PV,S,D,SSt,DSt}(
        init, mask, proc, opt, boundary, padval, source, dest, sourcestatus, deststatus
    )
end

Base.parent(d::WritableGridData) = parent(dest(d))
@propagate_inbounds function Base.setindex!(d::WritableGridData, x, I...)
    _setdeststatus!(d, x, I...)
    dest(d)[I...] = x
end

_setdeststatus!(d::WritableGridData, x, I...) = _setdeststatus!(d::WritableGridData, opt(d), x, I...)
function _setdeststatus!(d::WritableGridData{Y,X,R}, opt::SparseOpt, x, I...) where {Y,X,R}
    blockindex = _indtoblock.(I .+ R, 2R)
    @inbounds deststatus(d)[blockindex...] = !(opt.f(x))
    return nothing
end
_setdeststatus!(d::WritableGridData, opt::PerformanceOpt, x, I...) = nothing


text:
using Printf
using Grids

struct TrainingSample
    # @todo make training sample 3 fields: game, position, steps
    image::AbstractImage
    to_play::Integer # This information is encoded in the image, and is only needed for logging purposes.
    actions::Vector{<:Integer}
    targets::Vector{Target}
    image_label::String
    # @todo change image_state to initial_state
    image_state::AbstractEnvState
    states::Vector{AbstractEnvState}
    target_states::Vector{AbstractEnvState}
    function TrainingSample(image, to_play, actions, targets, states, target_states; image_label, image_state)
        @assert length(targets) == length(actions) + 1 "Num targets ($(length(targets))) must be 1 more than num actions ($(length(actions)))"
        @assert size(image)[end] == 3 "Image/Observation must have 3 channels (has $(size(image)[end])"
        new(image, to_play, actions, targets, image_label, image_state, states, target_states)
    end
end

struct TrainingSampleSVG
    sample::TrainingSample
end

svg(x::TrainingSample) = TrainingSampleSVG(x)

actions(x::TrainingSample) = x.actions
targets(x::TrainingSample) = x.targets
image(x::TrainingSample) = x.image
image_channels(x::TrainingSample) = x |> image |> channels

function TrainingSample(game::MuZeroGame, position::Integer, num_unroll_steps::Integer, td_steps::Integer)
    image = make_image(game, position)
    image_state = get_state(game, position)
    image_label = make_image_label(game, position)
    to = min(length(history(game)), position + num_unroll_steps - 1)
    actions = map(x->x.idx, history(game)[position:to])
    targets = make_target(game, position, num_unroll_steps, td_steps)
    num_targets = length(actions)+1
    player = to_play(position)
    states = map(idx->get_state(game, idx), position:(to+1))
    game_copy = copy_at(game.game, position)
    target_states = map(target->begin
        greedy_action = argmax(policy(target))
        step!(game_copy, greedy_action)
        get_state(game_copy)
    end, targets)
    TrainingSample(image, player, actions, targets[1:num_targets], states, target_states; image_label=image_label, image_state=image_state)
end

function TrainingSample(config::MuZeroConfig, game::MuZeroGame)
    TrainingSample(game, sample_position(game), config.num_unroll_steps, config.td_steps)
end

tuple(x::TrainingSample) = x.image, x.actions, x.targets

get_states(x::TrainingSample) = x.states
get_target_states(x::TrainingSample) = x.target_states

function getdf(x::TrainingSample)
    inftype_col = "Inf Type"
    action_col = "Action"
    reward_col = "Target Reward"
    value_col = "Target Value"
    policy_col = "Target Policy"
    d = Dict(
        inftype_col=>String[],
        action_col=>Union{Int,Nothing}[],
        reward_col=>Float32[],
        value_col=>Float32[],
        policy_col=>String[],
    )
    # @todo change this to x.action_space_size
    num_policy_actions = x.targets |> first |> policy |> length
    for i in 1:num_policy_actions
        d["π[$i]"] = Float32[]
    end

    states = get_states(x)
    actions = [nothing, x.actions...]
    targets = x.targets
    for (action_idx, zipped) in enumerate(zip(states, actions, targets))
        (s,a,target) = zipped
        push!(d[inftype_col], action_idx == 1 ? "Initial" : "Recurrent")
        push!(d[action_col], a)
        push!(d[reward_col], reward(target))
        push!(d[value_col], value(target))

        pred = policy(target)
        limit = 3
        from = 1
        to = min(limit, length(pred))
        actions = reverse(sortperm(pred))[from:to]
        prs = pred[actions]
        action_prs = map(x->begin
            action, pr = x
            pr_s = @sprintf "%.2f" pr
            "$action ($(pr_s[1:end]))"
        end, zip(actions, prs))
        result_policy = join(action_prs, ", ")
        push!(d[policy_col], result_policy)
        # π = policy(target)
        # for i in 1:num_policy_actions
        #     push!(d["π[$i]"], round(π[i], sigdigits=2))
        # end
    end

    DataFrame(
        inftype_col=>d[inftype_col],
        action_col=>d[action_col],
        reward_col=>d[reward_col],
        value_col=>d[value_col],
        policy_col=>d[policy_col],
    )
end

function getvrgrid(target)
    value = @sprintf "%.1f" target.value
    reward = @sprintf "%.1f" target.reward
    vr = Grid(
        [
            LabelRect("value: $value", width=200, halign=:left),
            LabelRect("reward: $reward", width=200, halign=:left),
        ],
        :,
        1,
        margin_left=0,
        line_thickness=0
    )
    Grid(
        [vr],
        :,
        1,
        header="TARGET V/R",
        margin_left=10,
        line_thickness=0
    )
end

function getinputgrid(s, header)
    player_name = s.to_play.name
    state_grid = getgrid(s; footer="to play: $player_name", margin_left=0)
    Grid(
        [state_grid],
        :,
        1,
        header=header,
        margin_left=0,
        line_thickness=0
    )
end

function getpolicygrid(s, target)
    pred = target.policy
    pred_states = Rect[]
    limit = 3
    from = 1
    to = min(limit, length(pred))
    pred_actions = reverse(sortperm(pred))[from:to]
    prs = pred[pred_actions]
    for (a,pr) in zip(pred_actions, prs)
        pr_str = @sprintf "%.2f" pr
        push!(
            pred_states,
            getgrid(
                s,
                a,
                0f0,
                getnextstate(s, a);
                margin_left=0,
                footer="action: $a",
                subfooter=(text=pr_str, color="#666"),
            )
        )
    end
    Grid(
        pred_states,
        :,
        1;
        header="TARGET π",
        margin_left=0,
        line_thickness=0,
    )
end

function getinferencegrid(s, target, header, input_header)
    input_grid = getinputgrid(s, input_header)
    vr_grid = getvrgrid(target)
    policy_grid = getpolicygrid(s, target)
    policy_col = Grid([Grid(3,3), Grid(3,3)], :, 1, header=(text="TARGET π", offset=(10,0)))
    Grid(
        [input_grid, vr_grid, policy_grid],
        header=header,
        margin_left=0,
        line_thickness=0,
    )
end

function getrecurrentinfgrid(sample)
    state_targets = zip(sample.states[2:end],sample.targets[2:end])
    map(st->begin
        s, target = st
        getinferencegrid(
            s,
            target,
            "RECURRENT INFERENCE",
            "HIDDEN STATE"
        )
    end, state_targets)
end

function Base.show(io::IO, m::MIME"image/svg+xml", x::TrainingSampleSVG)
    sample = x.sample
    initial_grid = getinferencegrid(
        sample.image_state,
        sample.targets[1],
        "INITIAL INFERENCE",
        "INPUT"
    )
    recurrent_grids = getrecurrentinfgrid(sample)
    g = Grid(
        [initial_grid, recurrent_grids...],
        :,
        1,
        header=(text="TRAINING SAMPLE"),
        line_thickness=0,
        margin=15,
    )
    show(io, m, g)
end

function Base.show(io::IO, x::TrainingSample)
    action_hl = Highlighter((data, i, j)->(j==2) && data[i,j] == nothing , crayon"#666666");
    pretty_table(
        io,
        getdf(x);
        show_row_number=true,
        nosubheader=true,
        title="Training Sample",
        highlighters=(action_hl),
    )
end

function gridsvg(x::TrainingSample)
    first_input_state, remaining_states... = get_states(x)
    first_target_state, remaining_target_states... = x.target_states
    first_target, remaining_targets... = x.targets
    state_targets = zip(get_states(x), x.targets)
    state_type = typeof(x.image_state)
    initial_grid = GridSVG([
        gridsvg(first_input_state),
        targetsvg(first_input_state, first_target_state, first_target)
    ], dim=2)
    grids = [
        initial_grid,
    ]
    for (input_state, action, target, target_state) in zip(remaining_states, x.actions, remaining_targets, remaining_target_states)
        push!(grids, GridSVG([
            gridsvg(input_state, title="hiddenstate + action")
            targetsvg(input_state, target_state, target)
        ], dim=2))
    end
    GridSVG(grids, main_title="Training Sample")
end

function get_player_label(x::TrainingSample, position)
    if x.to_play == 1
        return isodd(position) ? "A" : "B"
    end
    isodd(position) ? "B" : "A"
end


text:
module EMLpfr
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	include("pfr/pfr.jl")
end


text: function f(x,x) x end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# This file was generated, do not modify it. # hide
mysim = simulate(m(3), 1000)

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    sample(object, nsamples, [weights], replace=false)

Generate `nsamples` samples from spatial `object`
uniformly or using `weights`, with or without
replacement depending on `replace` option.
"""
function sample(object::Union{AbstractDomain,AbstractData}, nsamples::Int,
                weights::AbstractVector=[]; replace=false)
  if isempty(weights)
    sample(object, UniformSampler(nsamples, replace))
  else
    sample(object, WeightedSampler(nsamples, weights, replace))
  end
end


text:
using LightXML
import HTTP

#************************ LOCALS TO CONFIGURE!!!! **************************
const email= ""
#Enviroment variable that need to be setup
const umls_user = get(ENV, "UMLS_USER", "")
const umls_pswd = get(ENV, "UMLS_PSSWD", "")
const search_term="(obesity[MeSH Major Topic]) AND (\"2010\"[Date - Publication])"
const max_articles = 2
const overwrite_db=true
const verbose = false

#************************ MYSQL **********************************************
const host="127.0.0.1" #If want to hide - use enviroment variables instead
const mysql_usr="root"
const mysql_pswd=""
const dbname="pubmed_processes_test"
const dbname_pmid ="pmid_processes_test"
const medline_file = 1
const medline_year = 2018
#*****************************************************************************

const conn = DBUtils.init_mysql_database(host, mysql_usr, mysql_pswd, dbname)
PubMed.create_tables!(conn)

@testset "Search and Save" begin
    println("-----------------------------------------")
    println("       Testing Search and Save")

    pubmed_search_and_save!(email, search_term, max_articles,
    conn, verbose)

    #query the article table and make sure the count is correct
    all_pmids = PubMed.all_pmids(conn)
    @test length(all_pmids) == max_articles

    conn_pmid = DBUtils.init_mysql_database(host, mysql_usr, mysql_pswd, dbname_pmid)
    PubMed.create_pmid_table!(conn_pmid)
    pubmed_pmid_search_and_save!(email, search_term, max_articles, conn_pmid, verbose)

    #query the article table and make sure the count is correct
    all_pmids = BioMedQuery.PubMed.all_pmids(conn_pmid)
    @test length(all_pmids) == max_articles

    MySQL.disconnect(conn_pmid)

end

credentials_set = get(ENV, "TRAVIS_SECURE_ENV_VARS", "true")=="true"

@testset "UMLS" begin

    if credentials_set
        println("-----------------------------------------")
        println("       Testing MESH2UMLS")
        append = false

        map_mesh_to_umls_async!(conn, umls_user, umls_pswd; append_results=append, timeout=1)
        all_pairs_query = db_query(conn, "SELECT mesh FROM mesh2umls;")
        all_pairs = all_pairs_query[1]
        @test length(all_pairs) > 0

        println("-----------------------------------------")
        println("       Testing Occurrences")
        umls_concept = "Disease or Syndrome"

        labels2ind, occur = umls_semantic_occurrences(conn, umls_concept)

        @test length(keys(labels2ind)) > 0
        @test length(findall(x->x=="Obesity", collect(keys(labels2ind)))) ==1

    end
end

@testset "Medline Load" begin
    println("-----------------------------------------")
    println("       Testing Medline Loader")

    PubMed.create_tables!(conn) #drop and re-create pubmed article tables

    load_medline!(conn, @__DIR__, start_file=medline_file, end_file=medline_file, year=medline_year, test=true)

    path = joinpath(@__DIR__,"medline","raw_files",Processes.get_file_name(medline_file, medline_year, true))
    doc = parse_file(path)

    raw_articles = root(doc)

    all_pmids = PubMed.all_pmids(conn)
    @test length(all_pmids) == length(collect(child_elements(raw_articles)))
    res = MySQL.query(conn, "SELECT DISTINCT orcid FROM author_ref;", DataFrame)
    @test size(res)[1] > 2

    rm(joinpath(dirname(@__FILE__),"medline"), recursive=true)

end

MySQL.disconnect(conn)
mysql_conn = MySQL.connect(host, mysql_usr, mysql_pswd)
MySQL.execute!(mysql_conn, "DROP DATABASE IF EXISTS $dbname;")
MySQL.execute!(mysql_conn, "DROP DATABASE IF EXISTS $dbname_pmid;")
MySQL.disconnect(mysql_conn)

println("------------End Test Processes MySQL-----------")
println("-----------------------------------------------")


text:
# function predCurrFactorBeliefs(fgl::G,
#                                fc::Graphs.ExVertex  ) where G <: AbstractDFG
#   #
#   # TODO update to use ls and lsv functions
#   prjcurvals = Dict{String, Array{BallTreeDensity,1}}()
#   for v in getNeighbors(fgl, fc)
#     pred = kde!(evalFactor(fgl, fc, v.index))
#     curr = kde!(getVal(v))
#     prjcurvals[v.attributes["label"]] = [curr; pred]
#   end
#   return prjcurvals, collect(keys(prjcurvals))
# end

# function plotFactorBeliefs(fgl::G,
#                            flbl::Symbol ) where G <: AbstractDFG
#   #
#   if !haskey(fgl.fIDs, flbl)
#     println("no key $(flbl)")
#     return nothing
#   end
#   # for fc in fgl.f
#     # if fc[2].attributes["label"] == flbl

#     fc = fgl.g.vertices[fgl.fIDs[flbl]]  # fc = fgl.f[fgl.fIDs[flbl]]
#       prjcurvals, lbls = predCurrFactorBeliefs(fgl, fc)
#       if length(lbls) == 3
#         return vstack(
#         plotKDE(prjcurvals[lbls[1]]),
#         plotKDE(prjcurvals[lbls[2]]),
#         plotKDE(prjcurvals[lbls[3]]),
#         )
#       elseif length(lbls) == 2
#         return vstack(
#         plotKDE(prjcurvals[lbls[1]]),
#         plotKDE(prjcurvals[lbls[2]]),
#         )
#       elseif length(lbls) == 1
#         return plotKDE(prjcurvals[lbls[1]])
#       end

#     # end
#   # end
#   nothing
# end
# plotFactorBeliefs(fgl::G, flbl::T) where {G <: AbstractDFG, T <: AbstractString} = plotFactorBeliefs(fgl, Symbol(flbl))

# function whosWith(cliq::Graphs.ExVertex)
#   println("$(cliq.attributes["label"])")
#   for pot in cliq.attributes["data"].potentials
#       println("$(pot)")
#   end
#   nothing
# end

# function whosWith(bt::BayesTree, frt::String)
#     whosWith(whichCliq(bt,frt))
# end

# function plotUpMsgAtCliq(fg::G,
#                          cliq::Graphs.ExVertex  ) where G <: AbstractDFG
#   #
#   for id in keys(cliq.attributes["data"].debug.outmsg.p)
#       print("$(getVariable(fg,id).label), ") #fg.v[id].
#   end
#   println("")
#   sleep(0.1)
#   potens = getData(cliq).debug.outmsg.p
#   vArrPotentials(potens)
# end

# function plotUpMsgAtCliq(fg::G,
#                          bt::BayesTree,
#                          lbl::Union{String,Symbol}  ) where G <: AbstractDFG
#   #
#   plotUpMsgAtCliq(fg, whichCliq(bt, Symbol(lbl)) )
# end

# function dwnMsgsAtCliq(fg::G,
#                        cliq::Graphs.ExVertex  ) where G <: AbstractDFG
#   #
#   for id in keys(cliq.attributes["data"].debugDwn.outmsg.p)
#       print("$(getVariable(fg,id).label), ") # fg.v[id].
#   end
#   println("")
#   sleep(0.1)
#   potens = cliq.attributes["data"].debugDwn.outmsg.p
#   potens
# end

# function dwnMsgsAtCliq(fg::G,
#                        bt::BayesTree,
#                        lbl::Union{String,Symbol}  ) where G <: AbstractDFG
#   #
#   dwnMsgsAtCliq(fg, whichCliq(bt, Symbol(lbl)) )
# end

## =============================================================================
## Delete end v0.3.x
## =============================================================================

export drawPoses, drawLandms, drawPosesLandms
export drawSubmaps, drawMarginalContour, drawFeatTrackers, investigatePoseKDE
export drawHorBeliefsList, drawUpMsgAtCliq, drawFactorBeliefs
export drawPose2DMC!

# """
#     $(SIGNATURES)
#
# Standardize the length colors used by RoMEPlotting.
#
# Notes
# - Duplicated in KernelDensityEstimatePlotting
# """
# function getColorsByLength(len::Int=10)::Vector{String}
#   COLORS = String["red";"green";"blue";"black";"deepskyblue";"yellow";"magenta"]
#   morecyan = ["cyan" for i in (length(COLORS)+1):len]
#   retc = [COLORS; morecyan]
#   return retc[1:len]
# end

@deprecate drawPoses(x...; kwargs...) plotSLAM2DPoses(x...; kwargs...)
@deprecate drawLandms(x...; kwargs...) plotSLAM2DLandmarks(x...; kwargs...)
@deprecate drawPosesLandms(x...; kwargs...) plotSLAM2D(x...; kwargs...)
@deprecate drawSubmaps(x...; kwargs...) plotSLAM2DSubmaps(x...; kwargs...)
@deprecate drawMarginalContour(x...; kwargs...) plotMarginalContour(x...; kwargs...)
@deprecate drawFeatTrackers(x...; kwargs...) plotFeatTrackers(x...; kwargs...)
@deprecate drawHorBeliefsList(x...; kwargs...) plotHorBeliefsList(x...; kwargs...)
@deprecate drawUpMsgAtCliq(x...; kwargs...) plotUpMsgAtCliq(x...; kwargs...)
@deprecate drawFactorBeliefs(x...; kwargs...) plotFactorBeliefs(x...; kwargs...)
@deprecate drawPose2DMC!(x...; kwargs...) plotPose2DMC!(x...; kwargs...)

# deprecated
function investigatePoseKDE(p::BallTreeDensity, p0::BallTreeDensity)
    return investigateMultidimKDE(p, p0)
end

function investigatePoseKDE(p::Array{BallTreeDensity,1})
    return investigateMultidimKDE(p)
end

function investigatePoseKDE(p::BallTreeDensity)
    return investigateMultidimKDE(p)
end

@deprecate investigatePoseKDE(x...; kwargs...) plotPose(x...; kwargs...)

# function vArrPotentials(potens::Dict{Symbol,EasyMessage})
#   vv = Array{Gadfly.Compose.Context,1}(length(potens))
#   i = 0
#   oned=false
#   for p in potens
#       i+=1
#       pb = kde!(p[2].pts, p[2].bws)
#       if size(p[2].pts,1) > 3
#         # vv[i] = plotKDE(pb)
#         error("can't handle higher dimensional plots here yet")
#       elseif size(p[2].pts,1) > 1
#         vv[i] = investigateMultidimKDE(pb)
#       else
#         vv[i] = plotKDE(pb)
#       end
#   end
#   return vv
# end

# function draw(em::EasyMessage;xlbl="X")
#   p = Union{}
#   if size(em.pts,1) == 1
#     p=plotKDE(kde!(em),xlbl=xlbl)
#   else
#     p=plotKDE(kde!(em))
#   end
#   return p
# end

export drawOneMC!
export drawMCMCDebug
export drawTreeUpwardMsgs
export drawFrontalDens
export drawUpMCMCPose2D!
export drawDwnMCMCPose2D!
export drawLbl
export drawAnalysis
export drawAllPose2DBeliefs
export drawComicStripLM
export drawComicStrip

@deprecate drawOneMC!(x...; kwargs...) plotOneMC!(x...; kwargs...)
@deprecate drawMCMCDebug(x...; kwargs...) plotMCMCDebug(x...; kwargs...)
@deprecate drawTreeUpwardMsgs(x...; kwargs...) plotTreeUpwardMsgs(x...; kwargs...)
@deprecate drawFrontalDens(x...; kwargs...) plotFrontalDens(x...; kwargs...)
@deprecate drawUpMCMCPose2D!(x...; kwargs...) plotUpMCMCPose2D!(x...; kwargs...)
@deprecate drawDwnMCMCPose2D!(x...; kwargs...) plotDwnMCMCPose2D!(x...; kwargs...)
@deprecate drawLbl(x...; kwargs...) plotLbl(x...; kwargs...)
@deprecate drawAnalysis(x...; kwargs...) plotAnalysis(x...; kwargs...)
@deprecate drawAllPose2DBeliefs(x...; kwargs...) plotAllPose2DBeliefs(x...; kwargs...)
@deprecate drawComicStripLM(x...; kwargs...) plotComicStripLM(x...; kwargs...)
@deprecate drawComicStrip(x...; kwargs...) plotComicStrip(x...; kwargs...)

# function plotKDE(fgl::FactorGraph,
#                  vsym::Vector{Symbol};
#                  axis=nothing,
#                  dims=nothing,
#                  c=getColorsByLength(length(vsym)),
#                  levels::Int=4,
#                  title::Union{Nothing, T}=nothing,
#                  overlay=nothing  ) where {T <: AbstractString}
#   #
#   @warn "plotKDE for FactorGraph is deprecated, use DistributedFactorGraphs objects instead."
#   verts = map((x)->getKDE(getVariable(fgl, x)), vsym)
#   plotKDE(verts, dims=dims, c=c, axis=axis, levels=levels, title=title, overlay=overlay )
# end
#
# function plotKDE(fgl::FactorGraph,
#                  vsym::Symbol;
#                  axis=nothing,
#                  dims=nothing,
#                  c=nothing,
#                  levels=4,
#                  title::Union{Nothing, T}=nothing) where {T <: AbstractString}
#   #
#   @warn "plotKDE for FactorGraph is deprecated, use DistributedFactorGraphs objects instead."
#   plotKDE(fgl, Symbol[vsym;], dims=dims, c=c, axis=axis, levels=levels, title=title)
# end

#
## =============================================================================
## Old Functions -- Possibly Deprecated
## =============================================================================
function plotHorBeliefsList(fgl::G,
                            lbls::Array{Symbol,1};
                            nhor::Int=-1,
                            gt=nothing,
                            N::Int=200,
                            extend=0.1  ) where G <: AbstractDFG
  #
  len = length(lbls)
  pDens = BallTreeDensity[]
  for lb in lbls
    ptkde = getKDE(getVariable(fgl,lb))
    push!(pDens, ptkde )
  end

  if nhor<1
    nhor = round(Int,sqrt(len))
  end
  vlen = ceil(Int, len/nhor)
  vv = Array{Gadfly.Compose.Context,1}(undef, vlen)
  conslb = deepcopy(lbls)
  vidx = 0
  for i in 1:nhor:len
    pH = BallTreeDensity[]
    gtvals = Dict{Int,Array{Float64,2}}()
    labels = String[]
    for j in 0:(nhor-1)
      if i+j <= len
        push!(pH, pDens[i+j])
        push!(labels, string(lbls[i+j]))
        if gt != nothing gtvals[j+1] = gt[lbls[i+j]]  end
      end
    end
    vidx+=1
    if gt !=nothing
      vv[vidx] = KernelDensityEstimatePlotting.drawHorDens(pH, N=N, gt=gtvals, lbls=labels, extend=extend)
    else
      vv[vidx] = KernelDensityEstimatePlotting.drawHorDens(pH, N=N, lbls=labels, extend=extend)
    end
  end
  vv
end


text:
"""
    Search

Provides the types for objects and scenario management. This module also includes plot recipes for the various object types.

# Types
- `SimObject`: the abstract type that defines a basic interface for all objects within the sim environment.
    - `Target`: the subtype of `SimObject` that specializes for the search targets.
    - `Agent`: the subtype of `SimObject` that specializes for the searching agent.
    - `Sensor`: the subtype of `SimObject` that specializes for search sensors.
- `Scenario`: the type that contains an `Agent`, `Target`, `Sensors`, and other parameters.
"""
module Search

# using Distributions
# using StatsPlots
using Interpolations
using RecipesBase
using Plots

export Target, Agent, Sensor, Scenario, run!

"""
    SimObject

Abstract SimObject, which requires the following interface:
# Fields
- `f::Function`: the policy update function.
- `x::Array`: the position history (dimension N x T where N = length(position) and T = length(t)
- `v::Array`: the velocity history (dimension N x T where N = length(position) and T = length(t)
- `t::Array`: the time history (dimension T x 1)
"""
abstract type SimObject end

"""
    Sensor(f::Function, x::Array, v::Array, t::Array, r::Number)

Subtype of SimObject. It adds the field `r`.
# Fields
- `f::Function`: the policy update function.
- `x::Array`: the position history (dimension N x T where N = length(position) and T = length(t)
- `v::Array`: the velocity history (dimension N x T where N = length(position) and T = length(t)
- `t::Array`: the time history (dimension T x 1)
- `r::Number`: the detection radius around the sensor.
"""
mutable struct Sensor <: SimObject
    f::Function
    x::Array
    v::Array
    t::Array
    r::Real
end

"""
    Agent(f::Function, x::Array, v::Array, t::Array, slims::Tuple, s::Array)

Subtype of SimObject. It adds the field `slims` and `s`.
# Fields
- `f::Function`: the policy update function.
- `x::Array`: the position history (dimension N x T where N = length(position) and T = length(t)
- `v::Array`: the velocity history (dimension N x T where N = length(position) and T = length(t)
- `t::Array`: the time history (dimension T x 1)
- `slims::Tuple`: the upper and lower limits on the Agent speed.
- `s::Array`: the array of sensors available to the Agent.
"""
mutable struct Agent <: SimObject
    f::Function
    x::Array
    v::Array
    t::Array
    slims::Tuple
    s::Array
end

"""
    Target(f::Function, x::Array, v::Array, t::Array, slims::Tuple)

Subtype of SimObject. It adds the field `slims`.
# Fields
- `f::Function`: the policy update function.
- `x::Array`: the position history (dimension N x T where N = length(position) and T = length(t).
- `v::Array`: the velocity history (dimension N x T where N = length(position) and T = length(t).
- `t::Array`: the time history (dimension T x 1).
- `slims::Tuple`: the upper and lower speed limit of the target.
"""
mutable struct Target <: SimObject
    f::Function
    x::Array
    v::Array
    t::Array
    slims::Tuple
end

"""
    Scenario(a,t,s)

Type to hold the scenario parameters and objects.
# Fields
- `agent::Agent`
- `target::Target`
- `sensors::Array`
- `t::Tuple`
- `Δt::Real`
"""
mutable struct Scenario
    agent::Agent
    target::Target
    sensors::Array
    t::Tuple
    Δt::Real
end

"`Sensor(x,r)`: Initialize a Sensor with the default policy `policyDoNothing!` at time zero."
Sensor(x,r) = Sensor(policyDoNothing!, x, r)

"`Sensor(x,t,r)`: Initialize a Sensor with the default policy `policyDoNothing!` at time `t`."
Sensor(x,t,r) = Sensor(policyDoNothing!, x, [t], r)

"`Sensor(f::Function,x,r)`: Initialize a Sensor with zero initial velocity at time zero."
Sensor(f::Function,x,r) = Sensor(f, x, x.*0, [0.0], r)

"`Sensor(f::Function,x,t,r)`: Initialize a Sensor with zero initial velocity at time `t`."
Sensor(f::Function,x,t,r) = Sensor(f, x, x.*0, [t], r)

"Initialize an Agent with zero initial velocity and default speed limits."
Agent(f::Function,x, sensors; slims = (0,120)) = Agent(f, x, x.*0, [0.0], slims, sensors)

"`Target(f::Function,x,v,slims)`: Initialize a Target at time zero."
Target(f::Function,x,v,slims) = length(x)!=length(v) ? error("Position and Speed have unequal dimension.") : Target(f,x,v,[0.0],slims)
Target(x,v,slims) = Target(policyDoNothing!, x, v, slims)
Target(x; slims = (2,6)) = Target(x[1:length(x)/2], x[length(x)/2+1:end], slims)
"""
    Target(;slims, θlims)

Randomly initalize a Target at time zero with speed and direction drawn from uniform distribution over `slims` and `θlims`
"""
function Target(;slims=(2,6), θlims = (0,2π))
    x = [0.0, 0.0]
    θ = round(rand(range(θlims[1], stop = θlims[end], length = 100)), digits = 2)
    speed = round(rand(range(slims[1],stop = slims[end], length = 100)), digits = 2)
    v = speed .* [cos(θ),sin(θ)]
    Target(x,v,slims)
end

"""
    Scenario(::Agent, ::Target, ::Array, t, Δt)

Contain the objects within a scenario.
"""
function Scenario(t,Δt)
    a = Agent
end

"""
    step!(s::SimObject, Δt)

Propogate forward the position of a `SimObject` based on the velocity and time.
```
    x = x + v⋅Δt
    v = v
    t = t + Δt
```
"""
function step!(s::SimObject, Δt)
    s.t = vcat(s.t, s.t[end]+Δt)
    s.x = hcat(s.x, s.x[:,end] + s.v[:,end] .* Δt)
    s.v = hcat(s.v, s.v[:,end])
    return s
end

"`step!(s::Sensor,Δt)`: A specialized step! for `Sensor`, as they do not move."
function step!(s::Sensor, Δt)
    s.t = vcat(s.t, s.t[end]+Δt)
    s.x = hcat(s.x, s.x[:,end])
    s.v = hcat(s.v, s.v[:,end])
end

"""
    policyDoNothing!(s::SimObject)

The default policy to make no changes from the initial conditions.
"""
function policyDoNothing!(s::SimObject)
    return
end

function policySpiral!(a::Agent, Δt)
    a.v[:,end] = sum(a.slims)/4 * [cos(a.t[end]*2π), sin(a.t[end]*2π)]
end

"""
    policyTurnPeriodically!(t::Target)

Turn at regular intervals. Parameters are defined within the function body.
"""
function policyTurnPeriodically!(t::SimObject, Δt)
    if mod(t.t[end],.5)==0
        s = (2,6)
        Δθ = .5
        θ = atan(t.v[1,end],t.v[2,end])
        θ += 2*(rand()-.5)*Δθ
        speed = round(rand(range(s[1],stop = s[end], length = 100)), digits = 2)
        t.v = speed .* [cos(θ),sin(θ)]
    end
end

"`checkContact(t::Target,s::Sensor)::Bool`: Return `true` if the distance between `t` and `s` is less than or equal to detection range `s.r`, otherwise `false`."
checkContact(t::Target,s::Sensor)::Bool = distance(t,s) <= s.r ? true : false

"""
    distance(s1::SimObject, s2::SimObject)::Real

Return the distance between two `SimObject`'s.
# Input
- `s1::SimObject`
- `s2::SimObject`

# Output
The 2-norm distance between s1 and s2.
"""
function distance(s1::SimObject, s2::SimObject)::Real
    sqrt(sum((s1.x - s2.x).^2))
end

"""
    plot(t::Agent; interval = .25)

Plot the position history of an Agent.
"""
@recipe function plot(t::Agent; interval = .25)
    tick = t.t[1]:interval:t.t[end]
    tockx = LinearInterpolation(t.t,t.x[1,:])
    tocky = LinearInterpolation(t.t,t.x[2,:])
    @series begin
        linewidth := 0
        markershape --> :+
        tockx.(tick), tocky.(tick)
    end
    xguide --> "x (nm)"
    yguide --> "y (nm)"
    # title --> $(t.t[end])
    aspect_ratio --> :equal
    linecolor --> cgrad(:winter,rev=true)
    # colorbar --> true
    linewidth --> 2
    line_z --> t.t
    legend := false
    background_color_inside --> :lightblue
    n = length(t.t)
    seriesalpha --> (range(0, 1, length = n)).^2
    markershape --> :square
    markeralpha --> ifelse.(t.t .≈ 0, 1, 0)
    markercolor --> :red
    @series t.x[1,:], t.x[2,:]

end

"""
    plot(t::Target; interval = .25)

Plot the position history of a Target.
"""
@recipe function plot(t::Target;interval=.25)
    tick = t.t[1]:interval:t.t[end]
    tockx = LinearInterpolation(t.t,t.x[1,:])
    tocky = LinearInterpolation(t.t,t.x[2,:])
    @series begin
        linewidth := 0
        markershape --> :+
        tockx.(tick), tocky.(tick)
    end
    xguide --> "x (nm)"
    yguide --> "y (nm)"
    # title --> $(t.t[end])
    aspect_ratio --> :equal
    linecolor --> cgrad(:CMRmap,rev=true)
    # colorbar --> true
    linewidth --> 2
    line_z --> t.t
    legend := false
    background_color_inside --> :lightblue
    # n = length(t.t)
    # seriesalpha --> (range(0, 1, length = n))
    markershape --> :diamond
    markeralpha --> ifelse.(t.t .≈ 0, 1, 0)
    markercolor --> :yellow
    t.x[1,:], t.x[2,:]
end

"""
    run!(s::Scenario)

Loop through the Scenario time steps, apply policy and then step for each object
"""
function run!(s::Scenario)
    for t in s.t[1]:s.Δt:s.t[end]
        for obj in vcat(s.sensors, s.target, s.agent)
            policy!(obj)
        end
        for obj in vcat(s.sensors, s.target, s.agent)
            step!(obj)
        end
    end
end

# mutable struct Grid
#     x::AbstractVector
#     y::AbstractVector
#     v::AbstractMatrix
# end

# function Grid(x::AbstractVector, y::AbstractVector, d::Sampleable{Multivariate, Continuous})
#     v = [pdf(d,[i,j]) for i in x, j in y]
#     Grid(x,y,v)
# end

# mutable struct SearchScenario
#     initialDistribution::Sampleable{Multivariate, Continuous}
#     searchRadius::Number
#     searchSchedule::Array
#     AOU::Grid
# end

# function clip!(g::Grid,ϵ=1e-9)
#     g.v[abs.(g.v) ≈ 0] .= 0
# end

# function normalize!(g::Grid)
#     g.v = g.v ./ sum(g.v)
# end

# function clipNormalize!(g::Grid)
#     clip!(g)
#     normalize!(g)
# end

# import Base.+
# function +(a::Grid,b::Grid)
#     if a.x != b.x
#         error("x-axis does not match")
#     end
#     if a.y != b.x
#         error("y-axis does not match")
#     end
#     Grid(a.x, a.y, a.v .+ b.v)
# end

# import Base.*
# function *(a::Grid,b::Grid)
#     if a.x != b.x
#         error("x-axis does not match")
#     end
#     if a.y != b.x
#         error("y-axis does not match")
#     end
#     Grid(a.x, a.y, a.v .* b.v)
# end

# function initializeScenario(x::AbstractVector, y::AbstractVector, initDist::Sampleable{Multivariate, Continuous},searchRadius,searchSched::Array)::SearchScenario
#     SearchScenario(initDist,searchRadius,searchSched,Grid(x,y,initDist));
# end

# function initializeScenario(initDist::Sampleable{Multivariate, Continuous},searchRadius,searchSched::Array,model::Function;N = 1000)::SearchScenario
#     targets = Vector{Target}(undef,N)
#     for i in 1:N
#         targets[i] = Target(rand(initDist));
#     end
#     scenario = SearchScenario(initDist,searchDist,searchSched,model,targets);
# end

# function step!(scenario::SearchScenario)
#     map(scenario.targetModel!,scenario.targets)
#     dots = map(t -> t.x,targets)
#     p = scatter(dots[:,1],dots[:,2])
#     display(p)
# end

# function isFound(target::Target,search::MvNormal)
#     norm(target.x - search.μ) <= search.Σ[1]*1.177 ? true : false
# end

end # module


text:
"""
    ArcEager()

Arc-Eager transition system for dependency parsing.

# Transitions

| Transition  | Definition                                    |
|:----------- |:--------------------------------------------- |
| LeftArc(l)  | (σ\\|s, b\\|β, A) → (σ, b\\|β, A ∪ (b, l, s)) |
| RightArc(l) | (σ\\|s, b\\|β, A) → (σ, b\\|β, A ∪ (b, l, s)) |
| Reduce      | (σ\\|s, β,  A) → (σ, β,   A)                  |
| Shift	      | (σ,  b\\|β, A) → (σ\\|b, β, A)                |

# Preconditions

| Transition  | Condition                        |
|:----------- |:-------------------------------- |
| LeftArc(l)  | ¬[s = 0], ¬∃k∃l'[(k, l', i) ϵ A] |
| RightArc(l) | ¬∃k∃l'[(k, l', j) ϵ A]           |
| Reduce      | ∃k∃l[(k, l, i) ϵ A]              |

# References

[Nivre 2003](http://stp.lingfil.uu.se/~nivre/docs/iwpt03.pdf), [Nivre 2008](https://www.aclweb.org/anthology/J08-4003.pdf).
"""
struct ArcEager <: AbstractTransitionSystem end

initconfig(::ArcEager, graph::DependencyTree) =
    ArcEagerConfig(graph)
initconfig(::ArcEager, deptype, words) =
    ArcEagerConfig{deptype}(words)

transition_space(::ArcEager, labels=[]) =
    isempty(labels) ? [LeftArc(), RightArc(), Reduce(), Shift()] :
    [LeftArc.(labels)..., RightArc.(labels)..., Reduce(), Shift()]

projective_only(::ArcEager) = true

struct ArcEagerConfig{T} <: AbstractParserConfiguration{T}
    c::StackBufferConfiguration{T}
end

@stackbufconfig ArcEagerConfig

leftarc(cfg::ArcEagerConfig, args...; kwargs...) =
    ArcEagerConfig(leftarc_reduce(cfg.c, args...; kwargs...))

rightarc(cfg::ArcEagerConfig, args...; kwargs...) =
    ArcEagerConfig(rightarc_shift(cfg.c, args...; kwargs...))

reduce(cfg::ArcEagerConfig) = ArcEagerConfig(reduce(cfg.c))

shift(cfg::ArcEagerConfig) = ArcEagerConfig(shift(cfg.c))

isfinal(cfg::ArcEagerConfig) = all(a -> head(a) >= 0, cfg.c.A)
hashead(cfg::ArcEagerConfig, k) = head(token(cfg, k)) != -1

"""
    static_oracle(cfg, tree, transition)

Default static oracle function for arc-eager dependency parsing.

See [Goldberg & Nivre 2012](https://www.aclweb.org/anthology/C12-1059.pdf).
(Also called Arc-Eager-Reduce in [Qi & Manning 2017](https://nlp.stanford.edu/pubs/qi2017arcswift.pdf)).
"""
function static_oracle(cfg::ArcEagerConfig, gold_tree, arc)
    l = i -> arc(token(gold_tree, i))
    gold_arc = (a, b) -> has_arc(gold_tree, a, b)
    if stacklength(cfg) >= 1
        s = last(stack(cfg))
        if bufferlength(cfg) >= 1
            b = first(buffer(cfg))
            if gold_arc(b, s)
                return LeftArc(l(s)...)
            elseif gold_arc(s, b)
                return RightArc(l(b)...)
            end
        end
        if all(k -> k > 0 && hashead(cfg, k), [s ; dependents(gold_tree, s)])
            return Reduce()
        end
    end
    return Shift()
end

"""
    static_oracle_prefer_shift(cfg, gold_tree, arc=untyped)

Static oracle for arc-eager dependency parsing. Similar to the
"regular" static oracle, but always Shift when ambiguity is present.

See [Qi & Manning 2017](https://nlp.stanford.edu/pubs/qi2017arcswift.pdf).
"""
function static_oracle_prefer_shift(cfg::ArcEagerConfig, gold_tree, arc=untyped)
    l = i -> arc(token(gold_tree, i))
    gold_arc = (a, b) -> has_arc(gold_tree, a, b)
    (σ, s), (b, β) = popstack(cfg), shiftbuffer(cfg)
    if gold_arc(b, s)
        return LeftArc(l(s)...)
    elseif gold_arc(s, b)
        return RightArc(l(b)...)
    end
    must_reduce = false
    for k in stack(cfg)
        if gold_arc(k, b) || gold_arc(b, k)
            must_reduce = true
            break
        elseif head(token(cfg, k)) < 0
            break
        end
    end
    has_right_children = any(k -> s in rightdeps(gold_tree, k), buffer(cfg))
    if !must_reduce || s > 0 && !hashead(cfg, s) || has_right_children
        return Shift()
    else
        return Reduce()
    end
end

"""
    dynamic_oracle(t, cfg::ArgEagerConfig, tree)

Dynamic oracle function for arc-eager parsing.

For details, see [Goldberg & Nivre 2012](https://aclweb.org/anthology/C12-1059).
"""
dynamic_oracle(t, cfg::ArcEagerConfig, tree) = cost(t, cfg, tree) == 0

# see figure 2 in goldberg & nivre 2012 "a dynamic oracle..."
possible_transitions(cfg::ArcEagerConfig, graph::DependencyTree, arc=untyped) =
    possible_transitions(cfg, arc)

function possible_transitions(cfg::ArcEagerConfig, arc=untyped)
    s, b = last(stack(cfg)), first(buffer(cfg))
    l = i -> arc(token(cfg, i))
    transitions = [LeftArc(l(s)...), RightArc(l(b)...), Reduce(), Shift()]
    return filter(t -> is_possible(t, cfg), transitions)
end

function cost(t::LeftArc, cfg::ArcEagerConfig, gold)
    # left arc cost: num of arcs (k,l',s), (s,l',k) s.t. k ϵ β
    σ, s = popstack(cfg)
    b, β = shiftbuffer(cfg)
    if has_arc(gold, b, s)
        0
    else
        count(k -> has_arc(gold, k, s) || has_arc(gold, s, k), β)
    end
end

function cost(t::RightArc, cfg::ArcEagerConfig, gold)
    # right arc cost: num of gold arcs (k,l',b), s.t. k ϵ σ or k ϵ β,
    #                 plus num of gold arcs (b,l',k) s.t. k ϵ σ
    σ, s = popstack(cfg)
    b, β = shiftbuffer(cfg)
    if has_dependency(gold, s, b)
        0
    else
        count(k -> has_arc(gold, k, b), [σ ; β]) + count(k -> has_arc(gold, b, k), σ)
    end
end

function cost(t::Reduce, cfg::ArcEagerConfig, gold)
    # num of gold arcs (s,l',k) s.t. k ϵ b|β
    σ, s = popstack(cfg)
    count(k -> has_arc(gold, s, k), buffer(cfg))
end

function cost(t::Shift, cfg::ArcEagerConfig, gold)
    # num of gold arcs (k,l',b), (b,l',k) s.t. k ϵ s|σ
    b, β = shiftbuffer(cfg)
    count(k -> has_arc(gold, k, b) || has_arc(gold, b, k), stack(cfg))
end

function is_possible(::LeftArc, cfg::ArcEagerConfig)
    s = last(stack(cfg))
    return s != 0 && !hashead(token(cfg, s))
end

is_possible(::RightArc, cfg::ArcEagerConfig) =
    !hashead(token(cfg, first(buffer(cfg))))

is_possible(::Reduce, cfg::ArcEagerConfig) =
    hashead(token(cfg, last(stack(cfg))))

is_possible(::Shift, cfg::ArcEagerConfig) = true

==(cfg1::ArcEagerConfig, cfg2::ArcEagerConfig) = cfg1.c == cfg2.c

Base.getindex(cfg::ArcEagerConfig, i) = arc(cfg, i)


text:
module AccessControl

using HttpServer
using MbedTLS
using JSON
using ConnectionPools

import Redis.get, LoggedDicts.set!

export
    acpaths,                                      # Dict(path => handler) for authentication functionality
    session_create!, session_delete!,             # Common session functions
    session_read, session_write!,                 # Client-side sessions
    read_sessionid, session_get, session_set!,    # Server-side sessions
    notfound!, badrequest!, redirect!,            # Generic handlers
    login_form, pwdreset_form, logout_link, logout_pwdreset_links,    # AuthN HTML
    login!, logout!, process_pwdreset!,           # AuthN handlers
    is_logged_in, is_not_logged_in,               # AuthN utils
    has_role                                      # AuthZ utils

### Config - to be updated by the app's call to AccessControl.configure()
config            = Dict{Symbol, Any}()
config[:acdata]   = nothing                        # No default data store for access control data

config[:session]  = Dict(:datastore => :cookie,    # One of: :cookie::Symbol, ld::LoggedDict, cp::ConnectionPool
                         :cookiename => "id",      # Store the session (client-side) or session ID (server-side) in the "id" cookie
                         :id_length => 32,         # Length of the session ID in bytes
                         :max_n_sessions => 1,     # Max number of simultaneous sessions for a given user
                         :timeout => 600)          # Max number of seconds between requests in the same session

config[:login]    = Dict(:max_attempts => 5,           # Max number of allowed login attempts
                         :lockout_duration => 1800,    # Duration (sec) of lockout after max_attempts failed login attempts
			 :success_redirect => "/",     # Redirect location on successful login. Can be a function of the user.
			 :fail_msg => "Username and/or password incorrect.")    # Alert message on failed login attempt

config[:logout]   = Dict(:redirect => "/")         # Redirect location on logout

config[:pwdreset] = Dict(:max_attempts => 5,           # Max number of allowed attempts at password reset
                         :lockout_duration => 1800,    # Duration (sec) of lockout after max_attempts failed attempts
			 :success_redirect => "/",     # Redirect location on successful password reset. Can be a function of the user.
			 :fail_msg => "Password incorrect.")    # Alert message on failed password reset attempt

config[:securecookie]   = Dict{Symbol, Any}(:cookie_max_age => 10 * 60 * 1000,    # Duration of a session's validity in milliseconds
                                            :key_length     => 32,                # Key length for AES 256-bit cipher in CBC mode
                                            :block_size     => 16,                # IV  length for AES 256-bit cipher in CBC mode
                                            :cookie_attr    => Dict("Max-Age" => "600", "Secure" => "", "HttpOnly" => ""))

### Includes
# Common
include("utils.jl")
include("configure.jl")
include("generic_handlers.jl")
# Sessions
include("sessions/secure_cookies.jl")
include("sessions/clientside_sessions.jl")
include("sessions/serverside_common.jl")
include("sessions/serverside_loggeddict.jl")
include("sessions/serverside_redis.jl")
# Authentication
include("backends/backends_common.jl")
include("backends/backend_loggeddict.jl")
include("authentication/authn_html.jl")
include("authentication/authn_handlers.jl")
include("authentication/authn_utils.jl")
include("authentication/passwordhash/password_hash.jl")
include("authentication/passwordhash/pbkdf2.jl")
# Authorization
include("authorization/authz_utils.jl")

### Start-up scripts
update_securecookie_config!(:cookie_max_age)
update_securecookie_config!(:key_length)
update_securecookie_config!(:block_size)

# Access control paths
acpaths = Dict{AbstractString, Function}()
acpaths["/login"]            = login!
acpaths["/logout"]           = logout!              # Redirect user according to cfg
acpaths["/reset_password"]   = reset_password!      # Display password reset form
acpaths["/process_pwdreset"] = process_pwdreset!    # Handle credentials supplied by user as part of password reset

end # module


text:
using BrainFlow

# specify logging library to use
BrainFlow.enable_dev_logger(BrainFlow.BOARD_CONTROLLER)

params = BrainFlowInputParams()
board_shim = BrainFlow.BoardShim(BrainFlow.SYNTHETIC_BOARD, params)

BrainFlow.prepare_session(board_shim)
BrainFlow.start_stream(board_shim)
sleep(5)
BrainFlow.stop_stream(board_shim)
data = BrainFlow.get_current_board_data(32, board_shim)
BrainFlow.release_session(board_shim)

eeg_channels = BrainFlow.get_eeg_channels(BrainFlow.SYNTHETIC_BOARD)
sampling_rate = BrainFlow.get_sampling_rate(BrainFlow.SYNTHETIC_BOARD)

data_first_channel = data[eeg_channels[1], :]
println("Original Data First Channel")
println(data_first_channel)
BrainFlow.perform_rolling_filter(data_first_channel, 3, BrainFlow.MEAN)
println("After Rolling Filter")
println(data_first_channel)

data_second_channel = data[eeg_channels[2], :]
println("Original Data Second Channel")
println(data_second_channel)
BrainFlow.perform_wavelet_denoising(data_second_channel, "db4", 2)
println("After Wavelet Denoising")
println(data_second_channel)


text:
println("Playing Penney's Game Against the computer.")

if randbool()
    mach = autobet()
    println(@sprintf "The computer bet first, chosing %s." pgdecode(mach))
    println("Now you can bet.")
    human = pgencode(humanbet())
else
    println("You bet first.")
    human = pgencode(humanbet())
    mach = autobet(human)
end
print(@sprintf "You bet %s " pgdecode(human))
println(@sprintf "and the computer bet %s." pgdecode(mach))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Optim

export CFADistributionEM, CFADistributionEMRidge

"Confirmatory factor analysis distribution object with a Sigma formulation."
type CFADistributionEM
    Sigma_X::SparseMatrixCSC
    A::SparseMatrixCSC
    Sigma_L::AbstractMatrix
end
Base.length(d::CFADistributionEM) = Int(size(d.A)[1] + nnz(d.A) + (size(d.A)[2]-1)*size(d.A)[2]/2)

function Base.rand(d::CFADistributionEM, N::Int)
    P,K = size(d.A)
    distL = MvNormal(zeros(K), d.Sigma_L)
    L = rand(distL, N)
    X = d.A*L
    for i in 1:P
        X[i,:] .+= transpose(rand(Normal(0, d.Sigma_X.nzval[i]), N))
    end
    X
end
Base.rand(d::CFADistributionEM) = rand(d::CFADistributionEM, 1)

"Likelihood of model given the data represented by the covariance S from N samples"
function Distributions.loglikelihood(d::CFADistributionEM, S::AbstractMatrix, N::Int64)
    #@assert minimum(d.Sigma_X.nzval) > 0.0 "Invalid CFADistributionEM! (Theta_X has elements <= 0)"
    #@assert minimum(eig(d.Sigma_L)[1]) > 0.0 "Invalid CFADistributionEM! (Theta_L has eigen value <= 0)"

    Theta_X = deepcopy(d.Sigma_X)
    Theta_X.nzval[:] = 1 ./ d.Sigma_X.nzval
    Theta_L = inv(d.Sigma_L)

    # Woodbury transformation of original form allows us to only compute a KxK inversion
    Theta = Theta_X - Theta_X*d.A*inv(Theta_L + d.A'*Theta_X*d.A)*d.A'*Theta_X

    @assert minimum(eig(Theta)[1]) > 0.0 "Invalid CFADistributionEM! (Theta has eigen value <= 0)"

    logdet(Theta) - trace(S*Theta)
end

function normalize_Sigma_L!(d::CFADistributionEM)
    scaling = sqrt(diag(d.Sigma_L))
    Base.cov2cor!(d.Sigma_L, scaling)
    vals = nonzeros(d.A)
    K = size(d.A)[2]
    for col = 1:K
        for j in nzrange(d.A, col)
            vals[j] *= scaling[col]
        end
    end
end

" Run the EM algorithm to find the MLE fit."
function Distributions.fit_mle(::Type{CFADistributionEM}, maskSigma_L::AbstractMatrix, A::SparseMatrixCSC, S::AbstractMatrix, N::Int64; iterations=1000, show_trace=true, ftol=1e-10)
    P,K = size(A)
    @assert size(maskSigma_L) == (K,K) "Sigma_L mask must be of size K x K (where K = size(A)[2])"
    #@assert size(maskL) == (K,K) "Sigma_L mask must be of size K x K (where K = size(A)[2])"

    # just to make sure it is a nice float data type inside our loop
    maskL = zeros(Float64, K, K)
    copy!(maskL, maskSigma_L)

    # init our parameters
    d = CFADistributionEM(speye(P), deepcopy(A), eye(K))

    # find the latent vars each observed is attached to
    latentInds = Array{Int64}[find(d.A[i,:]) for i in 1:P]

    Theta_X = deepcopy(d.Sigma_X)
    lastScore = -Inf
    timeFindT = 0.0
    timeFindW = 0.0
    timeFindQ = 0.0
    timeUpdateSigma_X = 0.0
    timeUpdateA = 0.0
    timeUpdateSigma_L = 0.0
    stopCount = 0
    for count in 1:iterations

        ## E-step

        t = time()
        for i in 1:P
            Theta_X.nzval[i] = 1/d.Sigma_X.nzval[i]
        end
        Theta_L = inv(d.Sigma_L)
        Theta = Theta_X - Theta_X*d.A*inv(Theta_L + d.A'*Theta_X*d.A)*d.A'*Theta_X
        #@assert maximum(abs(Theta .- inv(d.Sigma_X + d.A*d.Sigma_L*d.A'))) < 1e-6
        T = Theta*d.A*d.Sigma_L
        timeFindT += time() - t

        # find W
        t = time()
        W = d.Sigma_L - d.Sigma_L*d.A'*Theta*d.A*d.Sigma_L
        timeFindW += time() - t

        # find Q (the expected covariance of the hidden variables)
        t = time()
        Q = T'*S*T + W
        #Base.cov2cor!(Q, sqrt(diag(Q)))
        timeFindQ += time() - t

        ## M-step

        # update A and Sigma_X
        t = time()
        tmpTS = T'*S
        # for i in 1:P
        #     gi = grouping[i]
        #     d.A[i,gi] = tmpTS[gi,i]/Q[gi,gi]
        #     d.Sigma_X[i,i] = S[i,i] - (tmpTS[gi,i]^2)/Q[gi,gi]
        # end
        #At = d.A'
        #rows = rowvals(At)
        #vals = nonzeros(At)
        for row = 1:P

            # update Sigma_X
            inds = latentInds[row]
            tmp = tmpTS[inds,row]
            #println(inds, " ", collect(nzrange(At, col)))
            newLoadings = inv(Q[inds,inds])*tmp
            #println(S[col,col] - tmpTS[inds,col]'*newLoadings)
            d.Sigma_X.nzval[row] = (S[row,row] - tmp'*newLoadings)[1]

            # Update A
            for (i,ind) in enumerate(inds)
                d.A[row,ind] = newLoadings[i]
            end
        end
        #d.A = At'
        timeUpdateA += time() - t

        # update Sigma_L
        t = time()
        d.Sigma_L[:,:] = Q .* maskL # note we enforce the Sigma_L mask here

        timeUpdateSigma_L += time() - t

        # print our progress
        if count % 1000 == 0 && show_trace
            s = loglikelihood(d, S, N)
            println("likelihood $count = ", s, " ", s - lastScore)

            println("timeFindT = $timeFindT")
            println("timeFindW = $timeFindW")
            println("timeFindQ = $timeFindQ")
            println("timeUpdateSigma_X = $timeUpdateSigma_X")
            println("timeUpdateA = $timeUpdateA")
            println("timeUpdateSigma_L = $timeUpdateSigma_L")
            println(d.Sigma_L[1,1], " ", d.Sigma_L[1,2])
            println()

            if s - lastScore < ftol
                stopCount += 1
            end
            if stopCount > 5
                break
            end
            lastScore = s
        end
    end
    normalize_Sigma_L!(d)
    d
end

type CFADistributionEMRidge
    rho::Float64
end

" Run the EM algorithm to find the MAP fit."
function fit_map(prior::CFADistributionEMRidge, ::Type{CFADistributionEM}, A::SparseMatrixCSC, S::AbstractMatrix, N::Int64; iterations=1000, show_trace=true, ftol=1e-10)
    P,K = size(A)

    # init our parameters
    d = CFADistributionEM(speye(P), deepcopy(A), eye(K))

    Theta_X = deepcopy(d.Sigma_X)
    lastScore = -Inf
    timeFindT = 0.0
    timeFindW = 0.0
    timeFindQ = 0.0
    timeUpdateSigma_X = 0.0
    timeUpdateA = 0.0
    timeUpdateSigma_L = 0.0
    stopCount = 0
    for count in 1:iterations

        ## E-step

        t = time()
        for i in 1:P
            Theta_X.nzval[i] = 1/d.Sigma_X.nzval[i]
        end
        Theta_L = inv(d.Sigma_L)
        Theta = Theta_X - Theta_X*d.A*inv(Theta_L + d.A'*Theta_X*d.A)*d.A'*Theta_X
        #Theta = inv(Sigma_X + A*Sigma_L*A')
        T = Theta*d.A*d.Sigma_L
        timeFindT += time() - t

        # find W
        t = time()
        W = d.Sigma_L - d.Sigma_L*d.A'*Theta*d.A*d.Sigma_L
        timeFindW += time() - t

        # find Q (the expected covariance of the hidden variables)
        t = time()
        Q = T'*S*T + W
        timeFindQ += time() - t

        ## M-step

        # update A and Sigma_X
        t = time()
        tmpTS = T'*S
        rows = rowvals(d.A)
        vals = nonzeros(d.A)
        for col = 1:K
            for j in nzrange(d.A, col)
                row = rows[j]
                d.Sigma_X.nzval[row] = S[row,row] - (tmpTS[col,row]^2)/Q[col,col]
                vals[j] = tmpTS[col,row]/Q[col,col]
            end
        end
        timeUpdateA += time() - t

        # update Sigma_L
        t = time()

        Base.cov2cor!(Q, sqrt(diag(Q)))
        d.Sigma_L[:,:] = Q + prior.rho*eye(K)
        Base.cov2cor!(d.Sigma_L, sqrt(diag(d.Sigma_L)))
        #Base.cov2cor!(d.Sigma_L, sqrt(prior.rho ./ (diag(Q) + prior.rho)))
        timeUpdateSigma_L += time() - t

        # print our progress
        if count % 1000 == 0 && show_trace
            s = loglikelihood(d, S, N)
            println("likelihood $count = ", s, " ", s - lastScore)

            println("timeFindT = $timeFindT")
            println("timeFindW = $timeFindW")
            println("timeFindQ = $timeFindQ")
            println("timeUpdateSigma_X = $timeUpdateSigma_X")
            println("timeUpdateA = $timeUpdateA")
            println("timeUpdateSigma_L = $timeUpdateSigma_L")
            println(d.Sigma_L[1,1], " ", d.Sigma_L[1,2])
            println()

            if abs(s - lastScore) < ftol
                stopCount += 1
            end
            if stopCount > 5
                break
            end
            lastScore = s
        end
    end
    normalize_Sigma_L!(d)
    d
end


text:
module NCTSSOS

using DynamicPolynomials
using MultivariatePolynomials
using JuMP
using MosekTools
using LightGraphs
using ChordalGraph
using MetaGraphs
using LinearAlgebra
using SparseArrays

export nctssos_first, nctssos_higher!, cs_nctssos_first, cs_nctssos_higher!

include("clique_merge.jl")
include("ncupop.jl")
include("nccpop.jl")
include("mixncpop.jl")

end


text:
using Memoization
using Test

##

@testset "Memoization" begin

    local n

    @memoize f(x::T, y=nothing, args... ; z::Union{Int,Nothing}, w=nothing, kwargs...) where {T<:Int} = (n+=1; ((x,y,args...), (;z=z,w=w,kwargs...)))

    # args
    @test (n=0; f(1,    z=1) == ((1,nothing),(z=1,w=nothing)) && n==1)
    @test (n=0; f(1,    z=1) == ((1,nothing),(z=1,w=nothing)) && n==0)

    @test (n=0; f(1,1,  z=1) == ((1,1),      (z=1,w=nothing)) && n==1)
    @test (n=0; f(1,1,  z=1) == ((1,1),      (z=1,w=nothing)) && n==0)

    @test (n=0; f(1,1,1,z=1) == ((1,1,1),    (z=1,w=nothing)) && n==1)
    @test (n=0; f(1,1,1,z=1) == ((1,1,1),    (z=1,w=nothing)) && n==0)

    # kwargs
    @test (n=0; f(1,z=1,w=1)     == ((1,nothing),(z=1,w=1))     && n==1)
    @test (n=0; f(1,z=1,w=1)     == ((1,nothing),(z=1,w=1))     && n==0)

    @test (n=0; f(1,z=1,w=1,p=1) == ((1,nothing),(z=1,w=1,p=1)) && n==1)
    @test (n=0; f(1,z=1,w=1,p=1) == ((1,nothing),(z=1,w=1,p=1)) && n==0)

    # cache clear
    Memoization.empty_cache!(f)
    @test (n=0; f(1,z=1) == ((1,nothing),(z=1,w=nothing)) && n==1)
    Memoization.empty_cache!(f)
    @test (n=0; f(1,z=1) == ((1,nothing),(z=1,w=nothing)) && n==1)

    # Dict vs. IdDict cache vs. custom cache
    @memoize IdDict g(x) = (n+=1; x)
    @test (n=0; g([1,2])==[1,2] && n==1)
    @test (n=0; g([1,2])==[1,2] && n==1)
    @memoize   Dict g(x) = (n+=1; x)
    @test (n=0; g([1,2])==[1,2] && n==1)
    @test (n=0; g([1,2])==[1,2] && n==0)
    @memoize Dict() g(x) = (n+=1; x)
    @test (n=0; g([1,2])==[1,2] && n==1)
    @test (n=0; g([1,2])==[1,2] && n==0)

    # redefinition
    # should only clear cache in toplevel, hence @eval
    @eval @memoize h(x) = x
    @test @eval h(2)==2 && h(2)==2
    @eval @memoize h(x) = 2x
    @test @eval h(2)==4 && h(2)==4

    # inference
    @test @inferred((@memoize foo(x) = x)(2)) == 2
    # this is broken because @inferred is failing despite @code_warntype giving that its inferred:
    @test_broken @inferred((@memoize foo(;x) = x)(x=2)) == 2

    # closures
    function make_func(x)
        @memoize func(y) = (n+=1; (x,y))
    end
    k = make_func(1)
    @test (n=0; k(2)==(1,2) && n==1)
    @test (n=0; k(2)==(1,2) && n==0)
    k′ = make_func(2)
    @test (n=0; k′(2)==(2,2) && n==1)
    @test (n=0; k′(2)==(2,2) && n==0)
    @test (n=0; k(2)==(1,2) && n==0)

    # callables
    @eval struct Baz{T} end
    @memoize (::Baz{T})(x::X) where {X, T<:Int} = (n+=1; Int)
    @memoize (::Baz{T})(x::X) where {X, T<:String} = (n+=1; String)
    bazint = Baz{Int}()
    bazstring = Baz{String}()
    @test (n=0; bazint(2)==Int && n==1)
    @test (n=0; bazint(2)==Int && n==0)
    @test (n=0; bazstring(2)==String && n==1)
    @test (n=0; bazstring(2)==String && n==0)
    Memoization.empty_cache!(Baz{Int})
    @test (n=0; bazint(2)==Int && n==1)
    Memoization.empty_cache!(Baz{String})
    @test (n=0; bazstring(2)==String && n==1)
    Memoization.empty_cache!(Baz)
    @test (n=0; bazint(2)==Int && n==1)
    @test (n=0; bazstring(2)==String && n==1)

    # unnamed args
    @memoize uarg(::Type{T}) where {T} = (n+=1; T)
    @test (n=0; uarg(Int)==Int && n==1)
    @test (n=0; uarg(Int)==Int && n==0)
    @test (n=0; uarg(Float64)==Float64 && n==1)
    @test (n=0; uarg(Float64)==Float64 && n==0)

end


text:
__precompile__()
module CosmoCentral

using QuadGK
using Base: @kwdef
using Parameters
using NumericalIntegration
using Interpolations
using PyCall
using Dierckx
using HDF5
using LoopVectorization
using Statistics
using JSON
using JSON3
using FFTW
using SpecialFunctions
using LatinHypercubeSampling
using Random
using DelimitedFiles
using LinearAlgebra
numpy = pyimport("numpy")

const classy = PyNULL()

function __init__()
    copy!(classy, pyimport("classy"))
end

include("CosmologicalStructures.jl")
include("Background.jl")
include("Density.jl")
include("Bias.jl")
include("MathUtils.jl")
include("WeightFunctions.jl")
include("BoltzmannSolver.jl")
include("PowerSpectrum.jl")
include("AngularCoefficients.jl")
include("FileIO.jl")
include("Forecaster.jl")
include("Derivator.jl")
include("LHSampler.jl")
include("FFTLog.jl")
include("SourceFunctions.jl")
include("TransferFunctions.jl")
include("IntrinsicAlignment.jl")
include("Covariance.jl")
include("Fisher.jl")
include("MCMCUtils.jl")

end # module


text:
module ODEOperatorsTests

using Gridap.ODEs.ODETools
using Test

import Gridap.ODEs.ODETools: test_ode_operator

include("ODEOperatorMocks.jl")

op = ODEOperatorMock{Float64,Constant}(1.0,2.0,3.0,1)

u = ones(2)
u_t = ones(2)*2.0

@assert(length(u) == 2)
@assert(length(u_t) == 2)

cache = allocate_cache(op)
update_cache!(cache,op,0.0)

r = allocate_residual(op,u,cache)
@test r == zeros(2)

J = allocate_jacobian(op,u,cache)
@test J == zeros(2,2)

t = 0.0
residual!(r,op,t,(u,u_t),cache)
_r = zeros(2)
_r[1] = u_t[1] - op.a * u[1]
_r[2] = u_t[2] - op.b * u[1] - op.c * u[2]
@test all(r .== _r)

J .= 0
jacobian!(J,op,t,(u,u_t),1,1.0,cache)
_J = zeros(2,2)
_J[1,1] = -op.a
_J[2,1] = -op.b
_J[2,2] = -op.c
@test all(J .== _J)

jacobian!(J,op,t,(u,u_t),2,1.0,cache)
_J[1,1] += 1.0
_J[2,2] += 1.0
@test all(J .== _J)
_J
J

@test test_ode_operator(op,t,u,u_t)

end #module


text:
function mark_kross_2!(r)
    for i in [Nord, West, Sud, Ost]
        go_to_border_and_draw(i)
    end
end

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs


text:
server = LanguageServerInstance(IOBuffer(), IOBuffer(), dirname(Pkg.Types.Context().env.project_file), first(Base.DEPOT_PATH))
server.runlinter = true
server.jr_endpoint = nothing

LanguageServer.initialize_request(init_request, server, nothing)

testtext = """
module testmodule
struct testtype
    a
    b::Float64
    c::Vector{Float64}
end
function testfunction(a, b::Float64, c::testtype)
    return c
end
end
testmodule
f(a,b,c) = 1
f()
f(1,)
f(1,2,)
"""
LanguageServer.textDocument_didOpen_notification(LanguageServer.DidOpenTextDocumentParams(LanguageServer.TextDocumentItem("testdoc", "julia", 0, testtext)), server, nothing)

doc = LanguageServer.getdocument(server, LanguageServer.URI2("testdoc"))
LanguageServer.parse_all(doc, server)

res = LanguageServer.textDocument_hover_request(LanguageServer.TextDocumentPositionParams(LanguageServer.TextDocumentIdentifier("testdoc"), LanguageServer.Position(3, 11)), server, nothing)
@test res.contents.value == string(LanguageServer.sanitize_docstring(StaticLint.CoreTypes.Float64.doc), "\n```julia\nCore.Float64 <: Core.AbstractFloat\n```")

res = LanguageServer.textDocument_hover_request(LanguageServer.TextDocumentPositionParams(LanguageServer.TextDocumentIdentifier("testdoc"), LanguageServer.Position(7, 12)), server, nothing)
@test occursin(r"c::testtype", res.contents.value)

res = LanguageServer.textDocument_hover_request(LanguageServer.TextDocumentPositionParams(LanguageServer.TextDocumentIdentifier("testdoc"), LanguageServer.Position(9, 1)), server, nothing)
@test res.contents.value == "Closes module definition for `testmodule`\n"

res = LanguageServer.textDocument_signatureHelp_request(LanguageServer.TextDocumentPositionParams(LanguageServer.TextDocumentIdentifier("testdoc"), LanguageServer.Position(12, 2)), server, nothing)
@test res.activeParameter == 0

res = LanguageServer.textDocument_signatureHelp_request(LanguageServer.TextDocumentPositionParams(LanguageServer.TextDocumentIdentifier("testdoc"), LanguageServer.Position(13, 4)), server, nothing)
@test res.activeParameter == 1

res = LanguageServer.textDocument_signatureHelp_request(LanguageServer.TextDocumentPositionParams(LanguageServer.TextDocumentIdentifier("testdoc"), LanguageServer.Position(14, 6)), server, nothing)
@test res.activeParameter == 2


text:
using Inpaintings, Test

my_tests = [
    "test_peaks.jl"
]

println("Running tests:")

@testset "Inpaintings.jl" begin
    @testset "$my_test" for my_test in my_tests
        println("  * $(my_test) *")
        include(my_test)
        println("\n\n")
    end
end

println("Tests finished")


text:
using Optim
using ParametricInversionTestUtils
using fdunc


text:
using Sockets
include("sender.jl")
include("os_updater.jl")
include("../strings.jl")

function use_cryptomat(p::Player)
    print_dict(p, "cryptomat_0")
    global note_max_length = 5
    while true
        user_input = readline(p.socket)
        if (user_input == "1") || (user_input == "3")
            sendSecret(p, parse(Int, user_input), p.msg)
            print_dict(p, "cryptomat_0")
        elseif user_input == "u"
            print_dict(p, "cryptomat_1")
            p.msg = readline(p.socket)
            print_dict(p, "cryptomat_0")
        elseif user_input == "c"
            p.msg = ""
        elseif user_input == "p"
            print_dict(p, "cryptomat_0")
            continue
        elseif user_input == "o"
            updateOS(p)
            print_dict(p, "cryptomat_0")
        elseif user_input == "g"
            print_dict(p, "crpytomat_generating_token")
            p.token = sendSecret(p, 1, p.msg, true)
            print_dict(p, "cryptomat_get_token")
        elseif user_input == "l"
            break
        elseif user_input == "◈"
            if !(isfile("data/.note"))
                print_dict(p, "cryptomat_2")
            else
                print_dict(p, "cryptomat_3")
                f = open_file_try("data/.note", "r")
                write(p.socket, "$(read(f, String))\n")
                close(f)
            end
        elseif user_input == "🕐"
            update = readline(p.socket)
            try
                p.dimension = parse(Int, update)
            catch
                continue
            end
        else
            print_dict(p, "cryptomat_4")
        end
    end
    print_dict(p, "cryptomat_5")
end


text:
function Range(arg0::GroundStation, arg1::jboolean, arg2::AbsoluteDate, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::ObservableSatellite)
    return Range((GroundStation, jboolean, AbsoluteDate, jdouble, jdouble, jdouble, ObservableSatellite), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function get_station(obj::Range)
    return jcall(obj, "getStation", GroundStation, ())
end

function is_two_way(obj::Range)
    return jcall(obj, "isTwoWay", jboolean, ())
end


text:
"""
Group alignment of functions using Bayesian method

    group_warping_bayes(f; iter=20000, times=5, powera=1)
    :param f: array (M,N) of N functions
    :param iter: number of MCMC iterations
    :param times: time slicing
    :param powera: MCMC parameter

    Returns Dict containing
    :return f_q: registered srvfs
    :return gam_q: warping function
    :return f_a: registered functions
    :return gam: warping functions
"""
function group_warping_bayes(f; iter=20000, times=5, powera=1)

    tau = ceil(times*.4);
    gp = collect(1:size(f,2));
    # Default setting shall work for many situations.
    # If convergence issues arise then adjust proposal variance tau.
    if (times == 2)
        warn("Small times may lead to convergence issues")
    end

    burnin = NaN;
    kappa = 1000;
    alpha = 1;
    beta = 0.001;
    var_const = 10000;
    thin = 1;
    scale_data = true;
    cut = 5*times;

    m = size(f, 1)-1;
    if mod(m, times) != 0
        error(@sprintf("Number of points on q function = %d is not a mulitple of times = %d", m, times))
    end
    timet = collect(linspace(0, 1, m+1));
    n = size(f, 2)
    qt_matrix = zeros(m, n);
    qt_fitted_matrix = zeros(m, n);

    for j = 1:n
        qt_matrix[:, j] = sign(diff(f[:,j])).*sqrt(abs(diff(f[:,j]))./diff(timet));
        if scale_data
            rescale = sqrt(m/sum(qt_matrix[:, j].^2));
            qt_matrix[:, j] = rescale.*qt_matrix[:, j];
        end
    end

    row = collect(1:times:m);
    L = length(row);

    res_dp = DP_mean(f, times);
    mu_5 = copy(res_dp["estimator2"]);
    mu_5 = vec(mu_5);
    match_matrix = copy(res_dp["match_matrix"]);

    MAP = copy(mu_5);
    best_match_matrix = copy(match_matrix);
    dist_vec = 100.*ones(n);
    best_vec = copy(dist_vec);
    sumdist = zeros(iter);
    kappa_collect = zeros(iter);
    log_collect = zeros(iter);
    logmax = 0;
    mu_prior  = ones(m);
    cov_prior = diagm(var_const.*ones(m));
    mu_q = zeros(round(Integer, iter/thin), m);
    mu_q_standard = zeros(mu_q);

    burnin = round(Integer, 0.5*iter/thin);
    AVG = length(collect(burnin:(round(Integer, iter/thin))));
    res = itermatch(iter, n, m, mu_5, match_matrix, qt_matrix,
                    qt_fitted_matrix, L, tau, times, kappa, alpha, beta,
                    powera, best_vec, dist_vec, best_match_matrix, mu_prior,
                    var_const, sumdist, thin, mu_q, mu_q_standard, logmax,
                    burnin, AVG);

    mu_q_standard = res["mu_q_standard"];
    mu_q = res["mu_q"];
    MAP = res["MAP"];
    best_match_matrix = res["best_match_matrix"];
    kappafamily = res["kappafamily"];
    sumdist = res["sumdist"];
    log_posterior = res["log_posterior"];
    dist = res["dist"]

    mu_est = mean(mu_q_standard[burnin:round(Integer, iter/thin),:], 1);
    rescale = sqrt(m/sum(mu_est.^2));
    mu_est *= rescale;
    mu_est2 = mean(mu_q[burnin:round(Integer, iter/thin),:], 1);
    rescale = sqrt(m/sum(mu_est2.^2));
    mu_est2 *= rescale;

    bayes_warps = res["bayes_warps"];
    gam_q = zeros(m+1, n);
    gam_a = zeros(m+1, n);
    f_a = zeros(m+1, n);
    f_q = zeros(m+1, n);

    for t = 1:n
         f_i = interpolate(([row;m+1],), best_match_matrix[:,t], Gridded(Linear()))
         gam_q[:, t] = f_i[1:(m+1)];
         f_o = interpolate((collect(1:size(f,1)),), f[:,t], Gridded(Linear()))
         tmp = f_o[linspace(1,m,times*(m+1)-1)];
         f_q[:, t] = tmp[round(Integer, (gam_q[:, t]-1)*times+1)];
         f_i = interpolate(([row;m+1],), bayes_warps[:,t], Gridded(Linear()))
         gam_a[:, t] = f_i[1:(m+1)];
         f_a[:, t] = tmp[round(Integer, (gam_a[:, t]-1)*times+1)];
    end

    out = Dict("f_q" => f_q, "gam_q" => gam_q, "gam_a" => gam_a, "f_a" => f_a);

end


text:
using DashBootstrapComponents, DashHtmlComponents

row = dbc_row(
    [
        dbc_col(html_div("One of three columns")),
        dbc_col(html_div("One of three columns")),
        dbc_col(html_div("One of three columns")),
    ],
    className = "g-0",
);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
"""
    manifest(c::Pkg.Types.Context)
Retrieves the manifest of a Context.
"""
manifest(c::Pkg.Types.Context) = c.env.manifest

"""
    project(c::Pkg.Types.Context)
Retrieves the project of a Context.
"""
project(c::Pkg.Types.Context) = c.env.project

"""
    isinproject(context, package::Union{String,UUID})
Checks whether a package is in the dependencies of a given context, e.g. is directly loadable.
"""
function isinproject end

"""
    isinmanifest(context, package::Union{String,UUID})
Checks whether a package is in the manifest of a given context, e.g. is either directly loadable or is a dependency of an loadable package.
"""
function isinmanifest end

"""
    find_parent(c, package::Union{String,UUID})
Finds all loadable packages for which `package` is a dependency.
"""
find_parent(c::Pkg.Types.Context, name::String, out = Set{UUID}()) = find_parent(c, packageuuid(c, name), out)
function find_parent(c::Pkg.Types.Context, uuid::UUID, out = Set{UUID}())
    for pkg in manifest(c)
        if uuid in values(deps(packageuuid(pkg), c))
            if isinproject(c, packagename(pkg))
                push!(out, packageuuid(pkg))
            else
                find_parent(c, packageuuid(pkg), out)
            end
        end
    end
    return out
end

@static if VERSION < v"1.1"
    isinmanifest(context::Pkg.Types.Context, module_name::String) = module_name in keys(manifest(context))
    isinmanifest(context::Pkg.Types.Context, uuid::UUID) = any(get(p[1], "uuid", "") == string(uuid) for (u,p) in manifest(context))

    isinproject(context::Pkg.Types.Context, package_name::String) = haskey(deps(project(context)), package_name)
    isinproject(context::Pkg.Types.Context, package_uuid::UUID) = any(u == package_uuid for (n,u) in deps(project(context)))

    function packageuuid(c::Pkg.Types.Context, name::String)
        for pkg in manifest(c)
            if first(pkg) == name
                return UUID(last(pkg)[1]["uuid"])
            end
        end
    end
    packageuuid(pkg::Pair{Any,Any}) = last(pkg) isa String ? UUID(last(pkg)) : UUID(first(last(pkg))["uuid"])
    packageuuid(pkg::Pair{String,Any}) = last(pkg) isa String ? UUID(last(pkg)) : UUID(first(last(pkg))["uuid"])

    function packagename(c::Pkg.Types.Context, uuid)
        for (n,p) in c.env.manifest
            if get(first(p), "uuid", "") == string(uuid)
                return n
            end
        end
        return nothing
    end

    function deps(uuid::UUID, c::Pkg.Types.Context)
        if any(p[1]["uuid"] == string(uuid) for (n,p) in manifest(c))
            return manifest(c)[string(uuid)][1].deps
        else
            return Dict{Any,Any}()
        end
    end
    deps(d::Dict{String,Any}) = get(d, "deps", Dict{String,Any}())
    deps(pe::PackageEntry) = get(pe[1], "deps", Dict{String,Any}())
    path(pe::PackageEntry) = get(pe[1], "path", nothing)
    version(pe::PackageEntry) = get(pe[1], "version", nothing)

    function frommanifest(c::Pkg.Types.Context, uuid)
        for (n,p) in c.env.manifest
            if get(first(p), "uuid", "") == string(uuid)
                return p
            end
        end
        return nothing
    end
else
    isinmanifest(context::Pkg.Types.Context, module_name::String) = any(p.name == module_name for (u,p) in manifest(context))
    isinmanifest(context::Pkg.Types.Context, uuid::UUID) = haskey(manifest(context), uuid)

    isinproject(context::Pkg.Types.Context, package_name::String) = haskey(deps(project(context)), package_name)
    isinproject(context::Pkg.Types.Context, package_uuid::UUID) = any(u == package_uuid for (n,u) in deps(project(context)))

    function packageuuid(c::Pkg.Types.Context, name::String)
        for pkg in manifest(c)
            if last(pkg).name == name
                return first(pkg)
            end
        end
    end
    packageuuid(pkg::Pair{String,UUID}) = last(pkg)
    packageuuid(pkg::Pair{UUID,PackageEntry}) = first(pkg)
    packagename(c::Pkg.Types.Context, uuid::UUID) = manifest(c)[uuid].name

    function deps(uuid::UUID, c::Pkg.Types.Context)
        if haskey(manifest(c), uuid)
            return deps(manifest(c)[uuid])
        else
            return Dict{String,Base.UUID}()
        end
    end
    deps(pe::PackageEntry) = pe.deps
    deps(proj::Pkg.Types.Project) = proj.deps
    deps(pkg::Pair{String,UUID}, c::Pkg.Types.Context) = deps(packageuuid(pkg), c)
    path(pe::PackageEntry) = pe.path
    version(pe::PackageEntry) = pe.version
    frommanifest(c::Pkg.Types.Context, uuid) = manifest(c)[uuid]
end

"""
    tryaccess(root::Module, target::Symbol, visited = Set{Module}())
Traverses all submodules of `root` to search for the module `target`.
"""
function tryaccess(root::Module, target::Symbol, visited = Set{Module}())
    root in visited && return nothing
    push!(visited, root)
    try
        return getfield(root, target)
    catch err
    end
    for n in names(root, all = true, imported = true)
        !isdefined(root, n) && continue
        x = getfield(root, n)
        if x isa Module
            y = tryaccess(x, target, visited)
            if y isa Module && nameof(y) == target
                return y
            end
        end
    end
    return nothing
end

function can_access(m::Module, s::Symbol)
    try
        return Base.eval(m, :($m.$s))
    catch
        return nothing
    end
end

function change_env(c, pe)
    if path(pe) isa String
        env_path = path(pe)
        Pkg.API.activate(env_path)
    elseif !is_stdlib(c, packageuuid(pe)) && ((Pkg.API.dir(packagename(pe)) isa String) &&!isempty(Pkg.API.dir(packagename(pe))))
        env_path = Pkg.API.dir(packagename(pe))
        Pkg.API.activate(env_path)
    end
end

function sha2_256_dir(path, sha = sha = zeros(UInt8, 32))
    (uperm(path) & 0x04) != 0x04 && return
    startswith(path, ".") && return
    if isfile(path) && endswith(path, ".jl")
        s1 = open(path) do f
            sha2_256(f)
        end
        sha .+= s1
    elseif isdir(path)
        for f in readdir(path)
            sha = sha2_256_dir(joinpath(path, f), sha)
        end
    end
    return sha
end

function sha_pkg(pe::PackageEntry)
    path(pe) isa String && isdir(path(pe)) && isdir(joinpath(path(pe), "src")) ? sha2_256_dir(joinpath(path(pe), "src")) : nothing
end

function _lookup(tr::PackageRef{N}, depot::Dict{String,ModuleStore}) where N
    if haskey(depot, tr.name[1])
        if N == 1
            return depot[tr.name[1]]
        else
            return _lookup(tr, depot[tr.name[1]], 2)
        end
    end
end

function _lookup(tr::PackageRef{N}, m::ModuleStore, i) where N
    if i < N && haskey(m.vals, tr.name[i])
        _lookup(tr, m.vals[tr.name[i]], i + 1)
    elseif i == N && haskey(m.vals, tr.name[i])
        return m.vals[tr.name[i]]
    end
end


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

JLLWrappers.@generate_main_file_header("flex")
JLLWrappers.@generate_main_file("flex", UUID("48a596b8-cc7a-5e48-b182-65f75e8595d0"))
end  # module flex_jll


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

JLLWrappers.@generate_main_file_header("libidn2")
JLLWrappers.@generate_main_file("libidn2", UUID("e3d30ef0-22f3-5ccc-b703-2e063d8d1f30"))
end  # module libidn2_jll


text:
@testset "delta" begin
    A = rand(10,5)
    B = rand(20,5)
    d = KernelFunctions.Delta()
    @test pairwise(d,A,dims=1) == Matrix(I,size(A,1),size(A,1))
    @test pairwise(d,A,B,dims=1) == zeros(size(A,1),size(B,1))
    @test d(1,2) == 0
    @test d(1,1) == 1
    @test_throws DimensionMismatch d(rand(3),rand(4))
end


text:
# using Base.Test
# using MIToS.Utils

print("""

Test eachline for ASCIIString
=============================
""")

let example = ""
  @test_throws ErrorException eachline(example)
end

let example = "Hola"
  @test collect(eachline(example)) == ["Hola"]
end

let example = "Hola\n"
  @test collect(eachline(example)) == ["Hola\n"]
end

let example = "\n"
  @test collect(eachline(example)) == ["\n"]
end

let example = "Hola\nMundo"
  @test collect(eachline(example)) == ["Hola\n", "Mundo"]
end

let example = "Hola\nMundo\n"
  @test collect(eachline(example)) == ["Hola\n", "Mundo\n"]
end

let example = "Hola\nMundo\n\n"
  @test collect(eachline(example)) == ["Hola\n", "Mundo\n", "\n"]
end