text:
"""
    flts(A::AbstractArray{<:Real}, y::Vector{<:Real})

Fast least trimmed squares by an algorithm from Rousseeuw et al. - Computing LTS Regression for Large Data Sets (1999)
The basic idea is to find an outlier-free subset `H` of regressors that minimizes the sum of squared residuals `Q`, resulting in a robust version of least squares up to 50% outliers

# Args

* A::AbstractArray{<:Real}: Vector or Matrix ∈ R(n x p) containing regressors
* y::Vector{<:Real}: Output Vector ∈ R(n)

# Keyword arguments

* h::Int: The size of the subset H, has to be between 0.5(n + p + 1)) <= h <= n and defaults to 0.5(n + p + 1)
* outliers::Real: Amount of outliers, between 0 <= outliers <= 0.5. A valid definition of h overwrites an outlier definition!
* N::Int: Number of initial p-subsets, defaults to 500 (minimum is 10). This is the main parameter regarding speed, to estimate a reasonable N, the probability of having at least one outlier free p-subset is given by: 1-(1-(1-outliers)^p)^N > 0. If verbose = true information about this gets printed.
* maxiter::Int: Maximum number of C-steps applied in the final optimization, defaults to 100
* ΔQmin::Real: Minimum change of the objective function (squared residuals) that defines convergence
* return_set::Bool: If true, the subset H and the regarding value Q are also returned in the form of (H, θ, Q)
* verbose::Bool: Print some more information

# Return
The parameters θ that minimize Q and optionally a tuple (H, θ, Q), which also contains the selected subset H and the minimal Q

# Example
Taken from the original paper
```jldoctest
julia> x = rand(Normal(0,10), 1000)
julia> y = x .+ 1 + randn(1000)
julia> y[801:end] = rand(Normal(0,5), 200)
julia> x[801:end] = rand(Normal(50,5), 200)
julia> A = [x ones(1000)]
julia> flts(xb,y, return_set = true)

([656, 336, 293, 763, 281, 248, 269, 231, 678, 372  …  734, 220, 434, 658, 174, 465, 745, 371, 759, 732], [1.0003968687209952, 1.0058265326498548], 750.0607838759548)
```
"""
function flts(A::AbstractArray{<:Real}, y::AbstractVector{<:Real};
    h::Integer = 0, outliers::Real = -1, N::Integer = 500, maxiter::Integer = 100, ΔQmin::Real = 0.0001, return_set::Bool = false, verbose::Bool = false)
    # Input checking
    n = length(y)
    size(A, 1) == n || throw(DimensionMismatch("Both inputs A and y should have the same number of rows"))
    N >= 10 || throw(DomainError("N needs to be >= 10"))
    if ndims(A) == 2
        p = size(A, 2)
    else
        # make sure matrix multiplication works
        p = 1
        A = reshape(A, :, 1)
    end

    # Get the subset length h
    if (round(0.5(n + p + 1)) <= h <= n)
        verbose && @info "h was set to: h = $h. Breakdown point is at ~$(round(100 - h / n * 100))% outliers"
    elseif (0.0 <= outliers <= 0.5)
        h = Int(round((1-outliers)*n))
        verbose && @info "h was set to: h = $h. Breakdown point is at ~$(round(100 - h / n * 100))% outliers"
    else
        h = Int(round(0.5(n + p + 1)))
        verbose && @info "h was set to default: h = $h. Breakdown point is at ~$(round(100 - h / n * 100))% outliers"
    end
    # Information about possible success
    verbose && @info "Chance to find an outlier-free p subset is at least ~$((1-(1-(h/n)^p)^N) * 100) %"

    # Initial subset H1 using p-subsets (Option b from paper)
    hs = fill(h, N)
    f(h) = get_initial_H(A, y, p, n, h)
    # Array of Tuples (H, θ, Q)
    initials = map(f, hs)

    # Carry out two C-steps, take top 10 candidates for further optimization
    g1(x) = optimize_H(A, y, h, x, 2, 0)
    opts = map(g1, initials)
    sort!(opts, by = x -> x[3])
    candidates = opts[1:10]

    # Optimze until convergence and return the best
    g2(x) = optimize_H(A, y, h, x, maxiter, ΔQmin)
    results = map(g2, candidates)
    sort!(results, by = x -> x[3])

    # Return the winner
    winner = results[1]
    if return_set
        return winner
    else
        return winner[2]
    end
end

# Helper function that applies a certain number of C-steps to a subset H and checks for convergence
function optimize_H(A, y, h::Int, initial, maxiter::Int, ΔQmin)
    Q_old = initial[3]
    # Create object in outer scope
    opt = nothing
    # Apply C-steps until convergence / maxiter
    for i in 1:maxiter
        opt = C_step(A, y, initial[2], h)
        ΔQ = Q_old - opt[3]
        if ΔQ < ΔQmin
            break
        end
        Q_old = opt[3]
    end
    return opt
end

# Helper function that gets an initial set H from a random p subset
function get_initial_H(A, y, p, n, h)
    inds = 1:n
    # sample a random p subset
    J = sample(inds, p, replace = false)
    # if subset does not specify a unique hyperplane add random observation
    i = 1
    while ((p+i+1) < n) && (rank(A[J, :]) < p)
        J = sample(inds, p+i, replace = false)
        i += 1
    end

    # Estimate first parameters
    θ_J = \(A[J, :], y[J])
    # Use a C-step for an initial H, θ, Q
    H_initial, θ_initial, Q_initial = C_step(A, y, θ_J, h)
    return H_initial, θ_initial, Q_initial
end

# The C-step, workinghorse that guarentees converging parameters
function C_step(A, y, θ_old, h)
    # Use the old parameters to get the residuals
    residuals = y .- A * θ_old
    # Sort the residuals and return a new subset with the smallest residuals
    # H_new = sortperm(abs.(residuals))[1:h] # should be more efficient, but apperantly is not
    H_new = sortperm(abs.(residuals))[1:h]
    # Get new parameters using oridinary least squares
    θ_new = \(A[H_new, :], y[H_new])
    Q_new = get_Q(A, y, H_new, θ_new)
    return H_new, θ_new, Q_new
end

# get the objective function Q
function get_Q(A, y, H, θ)
    residuals = y .- A * θ
    return sum(abs2, residuals[H])
end


text:
"""

    bmm(A, B)

Perform a batch matrix-matrix product of matrices stored in `A` and `B`. size(A,2) ==
size(B,1) and size(A)[3:end] and size(B)[3:end] must match.

If A is a (m,n,b...) tensor, B is a (n,k,b...) tensor, and the output is a (m,k,b...)
tensor.

"""
function bmm(a, b)
    sa,sb = size(a),size(b)
    @assert sa[2] == sb[1] && sa[3:end] == sb[3:end]
    a3,b3 = reshape(a,sa[1],sa[2],:), reshape(b,sb[1],sb[2],:)
    c3 = similar(a,sa[1],sb[2],size(a3,3))
    bmm!(a3, b3, c3)
    reshape(c3, sa[1], sb[2:end]...)
end

function bmm(a::KnetArray{T}, b::KnetArray{T}) where T
    sa,sb = size(a),size(b)
    @assert sa[2] == sb[1] && sa[3:end] == sb[3:end]
    a3,b3 = reshape(a,sa[1],sa[2],:), reshape(b,sb[1],sb[2],:)
    c3 = similar(a,sa[1],sb[2],size(a3,3))
    bmm!('N','N',one(T),a3,b3,zero(T),c3)
    reshape(c3, sa[1], sb[2:end]...)
end

function bmm!(transA::AbstractChar, transB::AbstractChar, alpha::Number, A::KnetArray{T}, B::KnetArray{T}, beta::Number, C::KnetArray{T}) where {T}
    cublasop(c::Char)=(if c=='N'; 0; elseif c=='T'; 1; elseif c=='C'; 2; else error("Unknown cublas op $c"); end)
    if ndims(A) != 3 || ndims(B) != 3
        throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    end
    ma,ka,bsa = size(A)
    kb,nb,bsb = size(B)
    if bsa != bsb
        throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    end
    bs = bsa
    if transA == 'N'
        m=ma; k=ka;
    else
        m=ka; k=ma;
    end
    if transB == 'N'
        k == kb || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
        n=nb; k==kb;
    else
        k == nb || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
        n=kb; k==nb;
    end
    (m == size(C,1) && n == size(C,2) && bs == size(C,3)) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    lda,ldb,ldc=ma,ka,ma
    transa = cublasop(transA); transb = cublasop(transB)
    alpha = T[alpha]; beta = T[beta]
    strideA, strideB, strideC = m*k, k*n, m*n
    if T<:Float64
        @cublas(cublasDgemmStridedBatched, (Cptr, UInt32, UInt32, Cint, Cint, Cint, Ptr{T}, Ptr{T}, Cint, Clonglong, Ptr{T}, Cint, Clonglong, Ptr{T}, Ptr{T}, Cint, Clonglong, Cint), cublashandle(), transa, transb, m, n, k, alpha, A, lda, strideA, B, ldb, strideB, beta, C, ldc, strideC, bs)
    elseif T<:Float32
        @cublas(cublasSgemmStridedBatched, (Cptr, UInt32, UInt32, Cint, Cint, Cint, Ptr{T}, Ptr{T}, Cint, Clonglong, Ptr{T}, Cint, Clonglong, Ptr{T}, Ptr{T}, Cint, Clonglong, Cint), cublashandle(), transa, transb, m, n, k, alpha, A, lda, strideA, B, ldb, strideB, beta, C, ldc, strideC, bs)
    else
        error("CUBLAS does not support $T")
    end
    return C
end

function bmm!(A, B, C)
    if ndims(A) != 3 || ndims(B) != 3
        throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    end
    ma,ka,bsa = size(A)
    kb,nb,bsb = size(B)

    (bsa == bsb && ka == kb && size(C,1) == ma && size(C,2) == nb && size(C,3) == bsa) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))

    for i=1:bsa
        C[:, :, i] = view(A, :, :, i) * view(B, :, :, i)
    end
    return C
end

@primitive bmm(x1,x2),dy,y  bmm(dy, permutedims(x2, (2,1,(3:ndims(x2))...)))  bmm(permutedims(x1, (2,1,(3:ndims(x1))...)), dy)
@zerograd bmm!(transA::AbstractChar, transB::AbstractChar, alpha::Number, A::KnetArray, B::KnetArray, beta::Number, C::KnetArray)
@zerograd bmm!(A, B, C)

ext: jl


text:
# older but much faster BVH
#

struct Triangles{IntT,FloatT}
  indices::Vector{IntT}
  P::Vector{SVector{3,FloatT}}
end

export
insert, query, makeroot, isleaf, visit, visitbf,
updateBVH, createBVH, @createBVH_nplex,
selfintersect, selfintersects, @selfintersect, @selfintersects,
DAABBNode, EmptyTree, triangle_intersect,
remove, getparent, getsibling, getleft, getright, isleft, isright, verify, getAABB

const Nullable{T} = Union{Nothing,T}
const isnull = isnothing

using Base.Iterators
#using SIMD
import Base.<

#========== macros ==========#

#== this will make a block of
quote
  baseindex_ = 3 * (j - 1)
  j1 = A[baseindex_ + 1]
  j2 = A[baseindex_ + 2]
  j3 = A[baseindex_ + 3]
end
==#
macro simplex_indices(A,dim,j)
  Expr(:block,
       Expr(:(=),
            Expr(:escape, Symbol("baseindex_")),
            Expr(:call, :*, eval(dim), Expr(:call, :-, Expr(:escape, j), 1))),
       map( i->Expr(:(=),
                    Expr(:escape, Symbol(j,i)),
                    Expr(:ref,
                         Expr(:escape, A),
                         Expr(:call, :+, Expr(:escape,Symbol("baseindex_")), i))
                   ),1:eval(dim))...,
       Expr(:(=),
            Expr(:escape, Symbol(j,"vertindices")),
            Expr(:tuple, map( i->Expr(:escape, Symbol(j,i)), 1:eval(dim))...)))
end

# use to reference symbols made by other macros in same calling context
macro mref(prefix,vname)
  Expr(:escape, Symbol(prefix,vname))
end

#==
transforms
simplex_P(P,3,i)
to
quote
  pi1 = P[i1]
  pi2 = P[i2]
  pi3 = P[i3]
end
==#
macro simplex_P(P,dim,j)
  Expr(:block,
       map( i->Expr(:(=),
                    Expr(:escape, Symbol("p",j,i)),
                    Expr(:ref,
                         Expr(:escape, P),
                         Expr(:escape,Symbol(j,i)))
                   ),1:eval(dim))...)
end

macro stup(prefix,dim)
  Expr(:tuple, map(i->Expr(:escape, Symbol(prefix,i)), 1:eval(dim))...)
end

macro sargs(prefix,dim)
  Expr(:..., Expr(:tuple, map(i->Expr(:escape, Symbol(prefix,i)), 1:eval(dim))...))
end

macro sargs(prefix,iter,dim)
  Expr(:..., Expr(:tuple, map(i->Expr(:escape, Symbol(prefix,iter,i)), 1:eval(dim))...))
end

macro ifenable( isenabled, jlti, body )
  if eval(isenabled)
    :(if $(esc(jlti))
        $(esc(body))
      end)
  else
    esc(body)
  end
end

#============  DynamicAABBTree =============#

mutable struct DAABBNode{T,F}
  aabb::AABB{F}
  data::T
  parent::Nullable{DAABBNode{T,F}}
  left::Nullable{DAABBNode{T,F}}
  right::Nullable{DAABBNode{T,F}}
end

"""
this special type indicates a tree with nothing in it.
insert will dispatch on this type or a DAABBNode
"""
struct EmptyTree
end

const ChildNode=Nullable{DAABBNode}

DAABBNode(aabb::AABB{F}, data::T, parent::DAABBNode{T,F}, childleft::DAABBNode{T,F}, childright::DAABBNode{T,F}) where {T,F} =
DAABBNode(aabb,data,parent,childleft,childright)

DAABBNode(aabb::AABB{F}, data::T, parent::DAABBNode{T,F}) where {T,F} =
DAABBNode(aabb,data,parent,nothing,nothing)

DAABBNode(aabb::AABB{F}, data::T) where {T,F} =
DAABBNode(aabb,data,nothing,nothing,nothing)

function getAABB(notree::EmptyTree)
  AABB(0.,0.,0.,0.,0.,0.)
end

function getAABB(tree::DAABBNode{T,F}) where {T,F}
  tree.aabb
end

function getparent(kid::EmptyTree)
  EmptyTree()
end

function getparent(kid::DAABBNode{T,F}) where {T,F}
  if isnull(kid.parent)
    EmptyTree()
  else
    kid.parent
  end
end

function getsibling(parent::DAABBNode{T,F}, kid::DAABBNode{T,F}) where {T,F}
  if !isnull(parent.right) && (parent.right == kid)
    parent.left
  else
    parent.right
  end
end

getsibling(parent::EmptyTree, kid::DAABBNode{T,F}) where {T,F} = EmptyTree()

getsibling(kid::DAABBNode{T,F}) where {T,F} = getsibling(getparent(kid),kid)

getleft(parent::EmptyTree) = EmptyTree()
getright(parent::EmptyTree) = EmptyTree()

getleft(parent::DAABBNode{T,F}) where {T,F} = isnull(parent.left) ? EmptyTree() : parent.left
getright(parent::DAABBNode{T,F}) where {T,F} = isnull(parent.right) ? EmptyTree() : parent.right

isleft(kid::DAABBNode{T,F}) where {T,F} = (getleft(getparent(kid)) == getsibling(getsibling(kid)))
isright(kid::DAABBNode{T,F}) where {T,F} = (getright(getparent(kid)) == getsibling(getsibling(kid)))

"""
return an empty tree
insert will dispatch on this type or a DAABBNode
"""
function makeroot()
  EmptyTree()
end

"""
insert a new leaf with data.  returns the new root
"""
function insert( ancestor::DAABBNode{T,F}, aabb::AABB{F}, fataabb::AABB{F}, data::T ) where {T,F}

  new_ancestor_aabb = union(fataabb,ancestor.aabb)

  if isleaf(ancestor)
    new_ancestor = DAABBNode(new_ancestor_aabb, zero(data))
    new_ancestor.parent = ancestor.parent
    new_ancestor.left = ancestor
    new_child = DAABBNode(aabb, data, new_ancestor)
    new_ancestor.right = new_child
    ancestor.parent = new_ancestor
    new_ancestor
  else
    ancestor.aabb = new_ancestor_aabb
    vol = volume(fataabb)

    #assert( !isnull(ancestor.left) ) # we keep tree in this state. no null children.
    #assert( !isnull(ancestor.right) )
    childL = ancestor.left
    childR = ancestor.right

    volleft = volume(union(fataabb, childL.aabb))
    volright = volume(union(fataabb, childR.aabb))

    if volleft < volright
      childL = insert(childL, aabb, fataabb, data)
      ancestor.left = childL
    else
      childR = insert(childR, aabb, fataabb, data)
      ancestor.right = childR
    end
    ancestor
  end

end

function insert( ancestor::DAABBNode{T,F}, aabb::AABB{F}, data::T, slop::F=0.1 ) where {T,F}
  fataabb = inflate(aabb,slop)
  insert(ancestor, aabb, fataabb, data)
end

function insert( root::EmptyTree, aabb::AABB{F}, data::T, slop::F=0.1 ) where {T,F}
  DAABBNode(aabb, data)
end

function refitup( root::EmptyTree, shapec::ST, slop::S ) where {ST,S}
  root
end

function refitup( ancestor::DAABBNode{T,F}, shapec::ST, slop::F ) where {T,F,ST}
  if isleaf(ancestor)
    ancestor.aabb = AABB( shapec, ancestor.data )
    refitup( getparent(ancestor), shapec, slop )
  else
    #assert(!isnull(ancestor.left))
    #assert(!isnull(ancestor.right))
    childL = ancestor.left
    childR = ancestor.right
    laabb = childL.aabb
    raabb = childR.aabb

    if isleaf(childL)
      laabb = AABB( shapec, childL.data )
      childL.aabb = laabb
      laabb = inflate( laabb, slop )
    end
    if isleaf(childR)
      raabb = AABB( shapec, childR.data )
      childR.aabb = raabb
      raabb = inflate( raabb, slop )
    end

    ancestor.aabb = union(laabb, raabb)
    refitup( getparent(ancestor), shapec, slop )
  end
end

function remove( root::EmptyTree, parent::EmptyTree, kid::DAABBNode{T,F} ) where {T,F}
  EmptyTree()
end

function remove( root::EmptyTree, parent::DAABBNode{T,F}, kid::DAABBNode{T,F} ) where {T,F}
  # pull sibling up to parents position
  brosis = getsibling(parent,kid)
  brosis.parent = nothing
  brosis
end

function remove( grandmapa::DAABBNode{T,F}, modad::DAABBNode{T,F}, kid::DAABBNode{T,F} ) where {T,F}
  # pull sibling up to parents position, grandparent will now raise the little bastard
  brosis = getsibling(modad,kid)
  brosis.parent = grandmapa

  if grandmapa.left == modad
    grandmapa.left = brosis
  else
    grandmapa.right = brosis
  end
  grandmapa
end

function remove( root::EmptyTree, kid::DAABBNode{T,F} ) where {T,F}
  EmptyTree()
end

function remove( parent::DAABBNode{T,F}, kid::DAABBNode{T,F} ) where {T,F}
  grandparent = getparent(parent)
  remove( grandparent, parent, kid )
end

function remove( kid::DAABBNode{T,F} ) where {T,F}
  parent = getparent(kid)
  grandparent = getparent(parent)
  remove( grandparent, parent, kid  )
end

function isleaf( child::DAABBNode{T,F} ) where {T,F}
  isnull(child.left) && isnull(child.right)
end

#==== caching flat tree ===#

"""
flatten the tree into a vector with leaf nodes first and branch nodes after
can be used to accelerate queries
"""
function flattentree(root::DAABBNode{T,F}, nleaf::T) where {T,F}
  leafcache = Vector{DAABBNode{T,F}}(undef, nleaf)
  visit(root, node->leafcache[node.data] = node)
  branchcache = Vector{DAABBNode{T,F}}()
  nodemap = Dict{DAABBNode{T,F},T}()
  for i in 1:nleaf
    nodemap[leafcache[i]] = i
  end

  id = nleaf + 1
  cachebranch(node) = begin
    node.data = id
    push!(branchcache,node)
    nodemap[node] = id
    id += 1
  end

  visit(root, identity, cachebranch)
  vcat(leafcache,branchcache)
end

create_childcache(flattree::Vector{DAABBNode{T,F}}) where {T,F} = map(node->( isnull(node.left) ? 0 : node.left.data, isnull(node.right) ? 0 : node.right.data), flattree);

create_aabbcache(flattree::Vector{DAABBNode{T,F}}) where {T,F} = map(node->node.aabb,flattree)

struct TreeCache{IntT,FloatT}
  aabbcache::Vector{AABB{FloatT}}
  nleaf::IntT
  childcache::Vector{Tuple{IntT,IntT}}
end

"""
create a flat memory cached verion of the tree for quick queries
"""
function cachetree( root::DAABBNode{T,F}, nleaf::T ) where {T,F}
  flattree = flattentree(root, nleaf)
  TreeCache( create_aabbcache(flattree),
             nleaf,
             create_childcache(flattree))
end

"""
return all leaves where aabb overlaps with the queried one. uses tree cache
known bug: a single node tree will miss that hit
"""
function query( rootid::T, aabb::AABB{F}, hits::Vector{T},
                      treecache::TreeCache{T,F} ) where {T,F}

  nleaf = treecache.nleaf

  nschunk = 32
  nscapacity = nschunk
  #!me should make a new type and use that... AutoVector
  nodestack = Vector{T}(undef,nscapacity)
  inodestack = 1
  nodestack[inodestack] = rootid

  while(inodestack > 0)
    ancestor = nodestack[inodestack]
    inodestack -= 1

    if overlaps(aabb,treecache.aabbcache[ancestor])

      if ancestor <= nleaf
        push!(hits,ancestor)
      else

        childL = treecache.childcache[ancestor][1]
        childR = treecache.childcache[ancestor][2]
        if childL > 0
          inodestack += 1
          nodestack[inodestack] = childL
        end

        if childR > 0
          if inodestack == nscapacity
            nscapacity *= 2
            nodestackprev = nodestack
            nodestack = Vector{T}(undef,nscapacity)
            copyto!(nodestack, 1, nodestackprev, 1, inodestack)
          end
          inodestack += 1
          nodestack[inodestack] = childR
        end
      end
    end
  end
end

"""
return all leaves where aabb overlaps with the queried one
known bug: a single node tree will miss that hit
"""
function query( root::DAABBNode{T,F}, aabb::AABB{F}, hits::Vector{T} ) where {T,F}
  nodestack = Vector{DAABBNode{T,F}}([root])

  while(!isempty(nodestack))
    ancestor = pop!(nodestack)

    if overlaps(aabb,ancestor.aabb)

      if isleaf( ancestor )
        push!(hits,ancestor.data)
      else

        childL = ancestor.left
        childR = ancestor.right
        if !isnull(childL)
          push!(nodestack, childL)
        end

        if !isnull(childR)
          push!(nodestack, childR)
        end
      end
    end
  end
end

"""
return all aabb's in the tree the overlap the query
"""
function query( ancestor::DAABBNode{T,F}, aabb::AABB{F} ) where {T,F}
  hits = Vector{T}()
  query( ancestor, aabb, hits )
  hits
end

"""
visit nodes in depth first order
"""
function visit( ancestor::DAABBNode{T,F}, leaffcn = show, branchfcn = identity ) where {T,F}
  if isleaf(ancestor)
    leaffcn(ancestor)
  else
    branchfcn(ancestor)
    if !isnull(ancestor.left)
      visit( ancestor.left, leaffcn, branchfcn )
    end
    if !isnull(ancestor.right)
      visit( ancestor.right, leaffcn, branchfcn )
    end
  end
end

"""
visit nodes in depth first order
general tree traversal with pre, post and inorder callback methods
"""
function visitdf( ancestor::DAABBNode{T,F}, pre, post = identity, inorder = identity ) where {T,F}
  pre(ancestor)
  if !isnull(ancestor.left)
    visitdf( ancestor.left, pre, post, inorder )
  end
  inorder(ancestor)
  if !isnull(ancestor.right)
    visitdf( ancestor.right, pre, post, inorder )
  end
  post(ancestor)
end

function visitreduce( ancestor::DAABBNode{T,F}, leaffcn, leaf0, branchfcn = (a,x)->a, branch0 = 0 ) where {T,F}
  if isleaf(ancestor)
    leaf0 = leaffcn(leaf0,ancestor)
  else
    branch0 = branchfcn(branch0,ancestor)
    if !isnull(ancestor.left)
      leaf0, branch0 = visitreduce( ancestor.left, leaffcn, leaf0, branchfcn, branch0 )
    end
    if !isnull(ancestor.right)
      leaf0, branch0 = visitreduce( ancestor.right, leaffcn, leaf0, branchfcn, branch0 )
    end
  end
  (leaf0,branch0)
end

function visitbfq( q::Vector{DAABBNode{T,F}},
                        leaffcn::Function, branchfcn::Function ) where {T,F}
  if !isempty(q)
    ancestor::DAABBNode{T,F} = shift!(q)
    if isleaf(ancestor)
      leaffcn(ancestor)
    else
      branchfcn(ancestor)
      if !isnull(ancestor.left)
        push!(q,ancestor.left)
      end
      if !isnull(ancestor.right)
        push!(q,ancestor.right)
      end
    end
    visitbfq( q, leaffcn, branchfcn )
  end
end

"""
visit each node in breadfirst order
"""
function visitbf( ancestor::DAABBNode{T,F}, leaffcn = show, branchfcn = identity ) where {T,F}
  visitbfq( [ancestor], leaffcn, branchfcn )
end

"""
return the depth of the tree
"""
function depth( node::DAABBNode{T,F}, curd = 1, minmax = max ) where {T,F}
  if isleaf(node)
    curd
  else
    minmax(depth(node.left, curd+1, minmax), depth( node.right, curd+1, minmax))
  end
end

#==
function printtree( q::Vector{DAABBNode{T,F}}, depth = 1, x = 1 ) where {T,F}
            if isempty(q)
               print("\n")
            else
              node = shift!(q)
              if !isleaf(node)
                  dl = depth( node.left )
                  dr = depth( node.right )
                  for i in 1:dl; print("       "); end
                  print("[  ]")
                  for i in 1:dr; print("       "); end
                  push!(q,node.left)
                  push!(q,node.right)
              else
                  if
                    ==#

"""
verify the tree satisfies bounding volume relations.  i.e. children are contained by parents
"""
function verify( node::DAABBNode{T,F}, P::Vector{Point{F}}, tris::Vector{T} ) where {T,F}
  containment = true
  if isleaf(node)
    containment = contains(node.aabb, AABB(node.data, P, tris))
  else
    if !isnull(node.left)
      containment = containment && contains(node.aabb, node.left.aabb)
      if containment
        containment = containment && verify( node.left, P, tris )
      end
      containment = containment && node.left.parent == node
    end
    if !isnull(node.right)
      containment = containment && contains(node.aabb, node.right.aabb)
      if containment
        containment = containment && verify( node.right, P, tris )
      end
      containment = containment && node.right.parent == node
    end
  end
  containment
end

#===== mesh stuff =====#
function createBVH( P::Vector{Point{F}}, tris::Vector{Tuple{T,T,T}}, sloppercent::F=F(0.05) ) where {F,T}
  slop = reduce(max,reduce(max,P)-reduce(min,P))*sloppercent
  root = makeroot()
  for ((i1,i2,i3), i) in zip(tris, 1:length(tris))
    root = insert(root, AABB( P[i1], P[i2], P[i3] ), i, slop)
  end
  root
end

function createBVH( P::Vector{Point{F}}, tris::Vector{T}, sloppercent::F=F(0.05) ) where {F,T}
  slop = reduce(max,reduce(max,P)-reduce(min,P))*sloppercent
  root = makeroot()
  for ((i1,i2,i3), i) in zip(partition(tris,3), 1:div(length(tris),3))
    root = insert(root, AABB( P[i1], P[i2], P[i3] ), i, slop)
  end
  root
end

function createBVH( trimesh::Triangles{IntT,FloatT}, sloppercent::FloatT=FloatT(0.05) ) where {IntT,FloatT}
  createBVH( trimesh.P, trimesh.indices, sloppercent )
end

macro createBVH_nplex(P, simplices, dim, sloppercent)
  quote
    slop = reduce(max,reduce(max,$(esc(P)))-reduce(min,$(esc(P))))*$(esc(sloppercent))
    root = makeroot()
    for (@stup(i,$dim), i) in zip(partition($(esc(simplices)),$dim), 1:div(length($(esc(simplices))),$dim))
      @simplex_P($(esc(P)),$dim,i)
      root = insert(root, AABB(@sargs("p",i,$dim)), i, slop)
    end
    root
  end
end

function points(tris::Triangles{IntT,FloatT},i::IntT) where {IntT,FloatT}
  ii = (i-1)*3;
  @inbounds r = (tris.P[tris.indices[ii+1]],tris.P[tris.indices[ii+2]],tris.P[tris.indices[ii+3]]); r
end

AABB( shapec::ST, ishape::IntT ) where {IntT<:Integer,ST} = AABB( points(shapec, ishape)... )

function updateBVH(node::DAABBNode{T,F},
                           shapec::ST, slop::F,
                           reinsert::Vector{DAABBNode{T,F}}) where {T,F,ST}
  #!me should compare with how box2d does it.  they do a very quick leaf only check.
  if isleaf(node)
    parent = getparent(node)
    parentAABB = getAABB(parent)
    currAABB = AABB(shapec, node.data)
    node.aabb = currAABB
    if contains(parentAABB,currAABB)
    else
      grandmapa = getparent(parent)
      remove( grandmapa, parent, node )
      refitup( grandmapa, shapec, slop )
      #insert( root, currAABB, node.data, slop )
      push!(reinsert,node)
    end
  else
    if !isnull(node.left)
      updateBVH(  node.left , shapec, slop, reinsert )
    end
    if !isnull(node.right)
      updateBVH(  node.right , shapec, slop, reinsert )
    end
  end

end

function updateBVH( node::DAABBNode{T,F}, shapec::ST, sloppercent::F = F(0.05) ) where {T,F,ST}
  #return createBVH( shapec )

  # create a foster grandparent to deal with depth 1 child removal ( will promote new root )
  knownuniverse = AABB(-Inf,Inf,-Inf,Inf,-Inf,Inf)
  fostergparent = DAABBNode( knownuniverse, 0 )
  fostergparent.left = node;

  # need well-formed children, so create a valid, but irrelevant child that fits in AABB of main node
  fostergparent.right = DAABBNode( AABB(shapec,1), 1 )
  node.parent = fostergparent

  P = shapec.P
  slop = reduce(max,reduce(max,P)-reduce(min,P))*sloppercent
  reinsert = Vector{DAABBNode{T,F}}()
  updateBVH( node, shapec, slop, reinsert )

  root = getleft(fostergparent)
  if (typeof(root) != EmptyTree)
    root.parent = nothing
  end

  #println("verify after refit/remove: ", verify(root, shapec))
  for inode in reinsert
     root = insert( root, inode.aabb, inode.data, slop )
  end
  #println("verify after reinsert: ", verify(root, shapec))
  root
end

macro selfintersect( BVH, Pa, indsa, dima, Pquery, indsq, dimq, collidef, skipduplicates )
  quote
    nleaf = div(length($(esc(indsa))),$dima)
    tc = cachetree($(esc(BVH)),nleaf)
    rootid = $(esc(BVH)).data

    nquery = div(length($(esc(indsq))),$dimq)
    hits = map(x->Vector{Tuple{eltype($(esc(indsa))),eltype($(esc(indsa))),Any}}(),1:nquery)
    Threads.@threads for i in 1:nquery
      @simplex_indices($(esc(indsq)),$dimq,i)
      @simplex_P($(esc(Pquery)),$dimq,i)

      elhits = Vector{typeof(rootid)}()

      query(rootid, AABB(@sargs("p",i,$dimq)), elhits, tc)

      for j in elhits
        # we can only skip duplicates if we are collding same simplex types on same mesh
        @ifenable $skipduplicates j < i begin
          @simplex_indices($(esc(indsa)),$dima,j)

          #!me could use @ifenable here to select between selfintersection and non self
          if isempty(intersect(@mref(i,"vertindices"),@mref(j,"vertindices")))
            @simplex_P($(esc(Pa)),$dima,j)

            #println("BVH check ", j, ", ", i)
            #hitfound = $(esc(collidef))(@sargs("p",i,$dimq),@sargs("p",j,$dima))
            #println("check hit: ", (@sargs("p",j,$dima), @sargs("p",i,$dimq))) #!me
            hit = $(esc(collidef))(@sargs("p",j,$dima), @sargs("p",i,$dimq))
            if first(hit)
              push!(hits[i],(j,i,hit))
            end
          end
        end
      end
    end
    reduce(vcat,hits)
  end
end

function selfintersect( BVH::DAABBNode{T,F}, triangles::Triangles{T,F}, tricollide::FN ) where {T,F,FN}
  P = triangles.P
  tris = triangles.indices
  @selfintersect( BVH, P, tris, 3, P, tris, 3, tricollide, true )
end

macro selfintersects( BVH, Pa, indsa, dima, Pquery, indsq, dimq, collidef, skipduplicates )
  quote
    nleaf = div(length($(esc(indsa))),$dima)
    tc = cachetree($(esc(BVH)),nleaf)
    rootid = $(esc(BVH)).data

    nquery = div(length($(esc(indsq))),$dimq)
    #for i in 1:ntri #nothreads
    hashit = map(x->false,1:nquery) #threads
    Threads.@threads for i in 1:nquery #threads
      @simplex_indices($(esc(indsq)),$dimq,i)
      @simplex_P($(esc(Pquery)),$dimq,i)

      elhits = Vector{typeof(rootid)}()

      #query(BVH, AABB(p1, p2, p3), elhits)
      query(rootid, AABB(@sargs("p",i,$dimq)), elhits, tc) #threads

      for j in elhits
        #!me ack, this optimization only works if we are colliding same types... otherwise, super wrong.
        #!me can selectively add with another macro pass...
        #if j < i
        @ifenable $skipduplicates j < i begin
          @simplex_indices($(esc(indsa)),$dima,j)

          if isempty(intersect(@mref(i,"vertindices"),@mref(j,"vertindices")))
            @simplex_P($(esc(Pa)),$dima,j)

            #hitfound = triangle_intersect(pi1,pi2,pi3, pj1, pj2, pj3,cfunc) #nothreads
            #hitfound = $(esc(collidef))(@sargs("p",i,$dimq),@sargs("p",j,$dima)) #threads
            hitfound = $(esc(collidef))(@sargs("p",j,$dima), @sargs("p",i,$dimq)) #threads
            #println("test:", @sargs("p",i,$dimq),@sargs("p",j,$dima), hitfound)
            if first(hitfound)
              hashit[i] = true
              #return true #nothreads
            end
          end
        end
      end
    end
    #return false #nothreads
    return reduce((accum,el)->accum || el,hashit; init=false) #threads
  end
end

function selfintersects( BVH::DAABBNode{T,F}, triangles::Triangles{T,F}, tricollide::FN ) where {T,F,FN}
  P = triangles.P
  tris = triangles.indices
  @selfintersects( BVH, P, tris, 3, P, tris, 3, tricollide, false )
end

function selfintersectdumb( P::Vector{Point{F}}, tris::Vector{T}, trintersect::FN ) where {T,F,FN}
  hits = Vector{Tuple{Int64,Int64,Any}}()
  for (vertindices, i) in zip(partition(tris,3), 1:div(length(tris),3))
    (i1,i2,i3) = vertindices
    p1 = P[i1]
    p2 = P[i2]
    p3 = P[i3]
    paabb = AABB(p1, p2, p3)

    for (vertindicesj, j) in zip(partition(tris,3), 1:div(length(tris),3))
      if i < j
        (j1,j2,j3) = vertindicesj
        p1j = P[j1]
        p2j = P[j2]
        p3j = P[j3]
        paabbj = AABB(p1j, p2j, p3j)
        if isempty(intersect(vertindices,(j1,j2,j3))) && overlaps(paabb,paabbj)
          hitfound = trintersect(p1,p2,p3, p1j,p2j,p3j)
          if first(hitfound)
            push!(hits,(i,j,hitfound))
          end
        end
      end
    end
  end

  unique(hits)
end


text:
# find i such that xs[i] is nearest to x
# ASSUMES xs is sorted (from least to largest) and xs[1] <= x
function findnearestindex(xs, x)
    n = length(xs)
    # using binary search algorithm O(log n), search the last element found
    # i.e., xs[i_first] <= x < xs[min(i_first + 1, n)] and i_first >= 1
    i_first = max(1, searchsortedlast(xs, x))
    i_next = min(i_first + 1, n)
    # suffices to compare xs[i_first] and xs[i_next]
    return ((x - xs[i_first]) < (xs[i_next] - x)) ? i_first : i_next
end


text:
#===============================================================================
                    Online visualization of Markov Chains
===============================================================================#

function _arrange_θs_to_plot(v::Val{:all}, θ_hist, indices, v1::Val{1})
    d = length(θ_hist[1][1])
    nu = length(θ_hist[1])
    idx_to_plot = 1:(nu*indices[end])
    [(idx_to_plot, _arrange_θs_to_plot(v, θ_hist, indices, i, v1)) for i in 1:d]
end

_arrange_θs_to_plot(v::Val{:all}, θ_hist, indices, i, v1::Val{1}) = map(
    θ->θ[i],
    collect(Iterators.flatten(θ_hist[indices]))
)

function _arrange_θs_to_plot(v::Val{:all}, θ_hist, indices, v1)
    d = length(θ_hist[1][1])
    nu = length(θ_hist[1])
    idx_to_plot = (nu*indices[1]-nu):(nu*indices[end])
    new_idx = (indices[1]-1):indices[end]
    [(idx_to_plot, _arrange_θs_to_plot(v, θ_hist, new_idx, i, v1)) for i in 1:d]
end

_arrange_θs_to_plot(v::Val{:all}, θ_hist, indices, i, v1) = map(
    θ->θ[i],
    collect(Iterators.flatten(θ_hist[indices]))[length(θ_hist[1]):end]
)

function _arrange_θs_to_plot(v::Val{:first}, θ_hist, indices, v1::Val{1})
    d = length(θ_hist[1][1])
    [(indices, map(θ->θ[1][i], θ_hist[indices])) for i in 1:d]
end

function _arrange_θs_to_plot(v::Val{:first}, θ_hist, indices, v1)
    d = length(θ_hist[1][1])
    new_idx = (indices[1]-1):indices[end]
    [(new_idx, map(θ->θ[1][i], θ_hist[new_idx])) for i in 1:d]
end

function _arrange_ar_to_plot(rolling_ar, indices, v1::Val{1})
    nu = length(rolling_ar[1])
    [(indices, map(ra->ra[i], rolling_ar[indices])) for i in 1:nu]
end

function _arrange_ar_to_plot(rolling_ar, indices, v1)
    num_updt = length(rolling_ar[1])
    new_idx = (indices[1]-1):indices[end]
    [(new_idx, map(ra->ra[i], rolling_ar[new_idx])) for i in 1:num_updt]
end

function _arrange_ll_to_plot(ll, indices, v1::Val{1})
    nu = length(ll[1])
    idx_to_plot = 1:(nu*indices[end])
    ( idx_to_plot, collect(Iterators.flatten(ll[indices])) )
end

function _arrange_ll_to_plot(ll, indices, v1)
    nu = length(ll[1])
    idx_to_plot = (nu*indices[1]-nu):(nu*indices[end])
    new_idx = (indices[1]-1):indices[end]
    ( idx_to_plot, collect(Iterators.flatten(ll[new_idx]))[nu:end] )
end

_const(ws::eMCMC.GlobalWorkspace) = (
    θ_len = length(ws.state),
    num_updts = length(ws.state_history[1]),
    N = length(ws.state_history)
)

@recipe function f(ws::eMCMC.GlobalWorkspace, v::K; indices=1:_const(ws).N) where {K<:Union{Val{:first}, Val{:all}}}
    layout --> (_const(ws).θ_len+1, 1)
    label --> ""
    #xguide --> (K <: Val{:first} ? "mcmc iteration" : "$(_const(ws).θ_len)*mcmc iteration")
    linecolor --> :steelblue
    (
        indices === nothing ?
        [([],[]) for _=1:_const(ws).θ_len+1] :
        vcat(
            _arrange_θs_to_plot(v, ws.state_history, indices, Val(indices[1])),
            _arrange_ll_to_plot(ws.ll_history, indices, Val(indices[1]))
        )
    )
end

@recipe function f(ws::eMCMC.GlobalWorkspace, v::Val{:ar}; indices=1:_const(ws).N)
    layout --> (_const(ws).num_updts, 1)
    label --> ""
    linecolor --> :steelblue
    xguide --> "mcmc iteration"
    (
        indices === nothing ?
        [([], []) for _=1:_const(ws).num_updts] :
        _arrange_ar_to_plot(ws.stats.rolling_ar, indices, Val(indices[1]))
    )
end
#=
using Plots
gr()
#pyplot()
plot_type = Val(:all)
num_refreshes = 1
num_segm = 99
segm_len = 100
updt_len = segm_len*num_segm
st = time()
begin
    for i in 1:num_refreshes
        p1 = plot(ws, plot_type; indices=1:((i-1)*updt_len+1))
        p2 = plot(ws, Val(:ar); indices=1:((i-1)*updt_len+1))
        p = plot(p1, p2, layout=(1,2), size=(1400,1000))
        d = display(p)
        println(typeof(p1))
        println("-------------")
        println(typeof(d))
        sleep(10)
        for k in ((i-1)*updt_len+1):segm_len:(i*updt_len)
            sleep(0.001)
            plot!(p1, ws, plot_type; indices=k:min((k+segm_len)))
            plot!(p2, ws, Val(:ar); indices=k:min((k+segm_len)))
            d = display(p)
        end
        i == num_refreshes || close(d)
    end
end
elapsed = time() - st
println("elapsed time: $elapsed.")

ws.ll_history

ws.acceptance_history

struct Foo
    data::Vector{Float64}
    t::Vector{Float64}
end

Base.length(foo::Foo) = length(foo.t)

@recipe function f(ws::Foo; idx_to_add=1:length(ws))
    ws.t[idx_to_add], ws.data[idx_to_add]
end

dt = 0.01
tt = collect(0.0:dt:1.0)
xx = zeros(Float64, length(tt))
xx[2:end]  = sqrt(dt) * randn(length(tt)-1)
cumsum!(xx, xx)
foo = Foo(xx, tt)

using Plots
pyplot()
begin
    p = plot(foo; idx_to_add=1:2, color=:steelblue)
    display(p)
    for i in 2:10
        sleep(1)
        plot!(p, foo; idx_to_add=i:(i+1), color=:steelblue)
        display(p)
    end
end

group = rand(map((i->begin
                    "group $(i)"
                end),1:4),100)
plot(rand(100),layout=@layout([a b;c]),group=group,linetype=[:bar :scatter :steppre])
=#


text:
"""
aspure(krs::Union{KRatio, KRatios}, mc::Type{<:MatrixCorrection} = XPP, fc::Type{<:FluorescenceCorrection} = NullFluorescence, cc::Type{<:CoatingCorrection} = NullCoating)

Reinterprete a k-ratio or k-ratios as relative to a pure, uncoated element at the same beam energy and take-off angle as the unknown.
"""
function aspure(
    krs::KRatio, #
    mc::Type{<:MatrixCorrection} = XPP,
    fc::Type{<:FluorescenceCorrection} = NullFluorescence,
    cc::Type{<:CoatingCorrection} = NullCoating,
)
    res = krs
    if (length(keys(krs.standard)) > 1) || (krs.unkProps[:TakeOffAngle]!=krs.stdProps[:TakeOffAngle]) || (krs.unkProps[:BeamEnergy]!=krs.stdProps[:BeamEnergy])
        pureprops = Dict( :TakeOffAngle => krs.unkProps[:TakeOffAngle], :BeamEnergy => krs.unkProps[:BeamEnergy] )
        purestd = NeXLCore.pure(krs.element)
        pure = zafcorrection(mc, fc, cc, purestd, krs.xrays, pureprops[:BeamEnergy], missing)
        std = zafcorrection(mc, fc, cc, krs.standard, krs.xrays, krs.stdProps[:BeamEnergy], get(krs.stdProps, :Coating, missing))
        kpure = k(pure, std, pureprops[:TakeOffAngle], krs.stdProps[:TakeOffAngle])
        res = KRatio(krs.xrays, krs.unkProps, pureprops, purestd, krs.kratio / kpure)
    end
    return res
end
function aspure(
    krs::KRatios, #
    mc::Type{<:MatrixCorrection} = XPP,
    fc::Type{<:FluorescenceCorrection} = NullFluorescence,
    cc::Type{<:CoatingCorrection} = NullCoating,
)
    res = krs
    if (length(keys(krs.standard)) > 1) || (krs.unkProps[:TakeOffAngle]!=krs.stdProps[:TakeOffAngle]) || (krs.unkProps[:BeamEnergy]!=krs.stdProps[:BeamEnergy])
        pureprops = Dict( :TakeOffAngle => krs.unkProps[:TakeOffAngle], :BeamEnergy => krs.unkProps[:BeamEnergy] )
        purestd = NeXLCore.pure(krs.element)
        pure = zafcorrection(mc, fc, cc, purestd, krs.xrays, pureprops[:BeamEnergy], missing)
        std = zafcorrection(mc, fc, cc, krs.standard, krs.xrays, krs.stdProps[:BeamEnergy], get(krs.stdProps, :Coating, missing))
        kpure = k(pure, std, pureprops[:TakeOffAngle], krs.stdProps[:TakeOffAngle])
        res = KRatios(krs.xrays, krs.unkProps, pureprops, purestd, krs.kratios / kpure)
    end
    return res
end

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
function domain end
function differentiate end
function integrate end
function localdot end

const ∂ = differentiate
const ∫ = integrate

differentiate(a::Union{Number,AbstractArray}) = zero(a)
localdot(a::T, b::T) where {T<:Union{Number,AbstractArray}} = dot(a, b)

struct DomainMismatch <: Exception end

Base.show(io::IO, ::DomainMismatch) =
    Base.print(io, "DomainMismatch(): function space arguments have non-matching domain.")

abstract type FunctionSpace{T} end

scalartype(::Type{<:FunctionSpace{T}}) where T = scalartype(T)
Base.eltype(::Type{<:FunctionSpace{T}}) where T = T

function LinearAlgebra.isapprox(x::FunctionSpace, y::FunctionSpace;
                                atol::Real=0,
                                rtol::Real=Base.rtoldefault(eltype(x), eltype(y), atol))
    return norm(x-y) <= max(atol, rtol*max(norm(x), norm(y)))
end


text:
module TestSKL

using Random
using Test
using AutoMLPipeline.AbsTypes
using AutoMLPipeline.SKLearners
using AutoMLPipeline.SKPreprocessors
using AutoMLPipeline.Utils
using AutoMLPipeline.BaseFilters
using AutoMLPipeline.Pipelines
using AutoMLPipeline.EnsembleMethods
using AutoMLPipeline.FeatureSelectors
using Statistics
using DataFrames

const IRIS = getiris()
const X = IRIS[:,1:3] |> DataFrame
const XC = IRIS[:,1:4] |> DataFrame
const YC = IRIS[:,5] |> Vector
const Y = IRIS[:,4] |> Vector

const classifiers = [
    "LinearSVC","QDA","MLPClassifier","BernoulliNB",
    "RandomForestClassifier","LDA",
    "NearestCentroid","SVC","LinearSVC","NuSVC","MLPClassifier",
    "RidgeClassifierCV","SGDClassifier","KNeighborsClassifier",
    "GaussianProcessClassifier","DecisionTreeClassifier",
    "PassiveAggressiveClassifier","RidgeClassifier",
    "ExtraTreesClassifier","GradientBoostingClassifier",
    "BaggingClassifier","AdaBoostClassifier","GaussianNB","MultinomialNB",
    "ComplementNB","BernoulliNB"
 ]

const regressors = [
    "SVR",
    "Ridge",
    "RidgeCV",
    "Lasso",
    "ElasticNet",
    "Lars",
    "LassoLars",
    "OrthogonalMatchingPursuit",
    "BayesianRidge",
    "ARDRegression",
    "SGDRegressor",
    "PassiveAggressiveRegressor",
    "KernelRidge",
    "KNeighborsRegressor",
    "RadiusNeighborsRegressor",
    "GaussianProcessRegressor",
    "DecisionTreeRegressor",
    "RandomForestRegressor",
    "ExtraTreesRegressor",
    "GradientBoostingRegressor",
    "MLPRegressor",
    "AdaBoostRegressor"
]

function fit_test(learner::String,in::DataFrame,out::Vector)
    _learner=SKLearner(Dict(:learner=>learner))
    fit!(_learner,in,out)
    @test _learner.model != Dict()
    return _learner
end

function fit_transform_reg(model::Learner,in::DataFrame,out::Vector)
    @test sum((transform!(model,in) .- out).^2)/length(out) < 2.0
end

@testset "scikit classifiers" begin
    Random.seed!(123)
    for cl in classifiers
	fit_test(cl,XC,YC)
    end
end

@testset "scikit regressors" begin
    Random.seed!(123)
    for rg in regressors
	model=fit_test(rg,X,Y)
	fit_transform_reg(model,X,Y)
    end
end

end


text:
module LibSndFile

import SampledSignals
using SampledSignals: SampleSource, SampleSink, SampleBuf
using SampledSignals: PCM16Sample, PCM32Sample
using SampledSignals: nframes, nchannels, samplerate
using FileIO: File, Stream, filename, stream
using FileIO: add_format, add_loader, add_saver, @format_str
using Printf: @printf
using LinearAlgebra: transpose!
using libsndfile_jll: libsndfile

const supported_formats = (format"WAV", format"FLAC", format"OGG")

include("libsndfile_h.jl")
include("lengthIO.jl")
include("sourcesink.jl")
include("loadsave.jl")
include("readwrite.jl")

function __init__()
    # ogg currently not in the registry
    add_format(format"OGG", "OggS", [".ogg", ".oga"], [:LibSndFile])
    for fmt in supported_formats
        add_loader(fmt, :LibSndFile)
        add_saver(fmt, :LibSndFile)
    end
end

end # module LibSndFile


text:
@testset "AbstractValEdge" begin
    @test ValEdge <: AbstractValEdge
    @test ValDiEdge <: AbstractValEdge
end

@testset "Edge Constructors" for
    V       in TEST_VERTEX_TYPES_SMALL,
    E_VALS  in TEST_EDGEVAL_TYPES_SMALL,
    u       in V[1, 2, typemax(V)],
    v       in V[1, 2, typemax(V)]

    values = rand_sample(E_VALS)

    # We enforce that the constructor always puts the smaller vertex as src
    @testset "ValEdge($V($u), $V($v), $values)" begin

        src_should, dst_should = minmax(u, v)
        e = ValEdge(u, v, values)
        e_rev = ValEdge(v, u, values)

        @test e.src == e_rev.src == src_should
        @test e.dst == e_rev.dst == dst_should
        @test e.src <= e.dst
        @test e.values == values
        @test e isa ValEdge{V, E_VALS}
    end

    @testset "ValEdge{$V}($u, $v, $values)" begin

        src_should, dst_should = minmax(u, v)
        e = ValEdge{V}(u, v, values)
        e_rev = ValEdge{V}(v, u, values)

        @test e.src == e_rev.src == src_should
        @test e.dst == e_rev.dst == dst_should
        @test e.src <= e.dst
        @test e.values == values
        @test e isa ValEdge{V, E_VALS}
    end

    @testset "ValDiEdge($V($u), $V($v), $values)" begin

        e = ValDiEdge(u, v, values)

        @test e.src == u
        @test e.dst == v
        @test e.values == values
        @test e isa ValDiEdge{V, E_VALS}
    end

    @testset "ValDiEdge{$V}($u, $v, $values)" begin

        e = ValDiEdge{V}(u, v, values)

        @test e.src == u
        @test e.dst == v
        @test e.values == values
        @test e isa ValDiEdge{V, E_VALS}
    end
end

@testset "edge access functions" for
    V       in TEST_VERTEX_TYPES_SMALL,
    E_VALS  in TEST_EDGEVAL_TYPES_SMALL,
    u       in V[1, 2],
    v       in V[1, 2]

    values = rand_sample(E_VALS)

    e = ValEdge(u, v, values)
    e_rev = ValEdge(v, u, values)
    @testset "e == ValEdge($V($u), $V($v), $values)" begin

        @testset "src, dst" begin
            @test issetequal((src(e), dst(e)), (u, v))
            @test issetequal((src(e_rev), dst(e_rev)), (u, v))
            @test src(e) isa V
            @test dst(e) isa V
        end

        @testset "get_edgeval(e, :)" begin
            @test get_edgeval(e, :) == values
            @test get_edgeval(e, :) isa E_VALS
        end

        @testset "get_edgeval(e, $key)" for
            key in (keys(values) ∪ Base.OneTo(length(values))) # number & symbols

            @test get_edgeval(e, key) == values[key]
            key isa Integer && @test get_edgeval(e, key) isa E_VALS.types[key]
        end
    end

    e = ValDiEdge(u, v, values)
    @testset "e == ValDiEdge($V($u), $V($v), $values)" begin

        @testset "src, dst" begin
            @test src(e) == u
            @test dst(e) == v
            @test src(e) isa V
            @test dst(e) isa V
        end

        @testset "get_edgeval(e, :)" begin
            @test get_edgeval(e, :) == values
            @test get_edgeval(e, :) isa E_VALS
        end

        @testset "get_edgeval(e, $key)" for
            key in (keys(values) ∪ Base.OneTo(length(values))) # number & symbols

            @test get_edgeval(e, key) == values[key]
            @test get_edgeval(e, key) == get_edgeval(e, :)[key]
            key isa Integer && @test get_edgeval(e, key) isa E_VALS.types[key]
        end
    end

end

@testset "reverse edge" begin

    @testset "edge == $e" for
        V       in TEST_VERTEX_TYPES_SMALL,
        E_VALS  in TEST_EDGEVAL_TYPES_SMALL,
        u       in V[1, 2],
        v       in V[1, 2],
        e       in ( ValEdge(u, v, rand_sample(E_VALS)), )

        e_rev = reverse(e)
        @test src(e) == src(e_rev)
        @test dst(e) == dst(e_rev)
        @test get_edgeval(e, :) == get_edgeval(e_rev, :)
        @test typeof(e) == typeof(e_rev)
    end

    @testset "edge == $e" for
        V       in TEST_VERTEX_TYPES_SMALL,
        E_VALS  in TEST_EDGEVAL_TYPES_SMALL,
        u       in V[1, 2],
        v       in V[1, 2],
        e       in ( ValDiEdge(u, v, rand_sample(E_VALS)), )

        e_rev = reverse(e)
        @test src(e) == dst(e_rev)
        @test dst(e) == src(e_rev)
        @test get_edgeval(e, :) == get_edgeval(e_rev, :)
        @test typeof(e) == typeof(e_rev)
    end
end

@testset "hash" begin

    @test hash(ValEdge(1, 2, (a=3, b=4))) == hash(ValEdge(1, 2, (a=3, b=4)), UInt(0))
    @test hash(ValDiEdge(1, 2, (a=3, b=4))) == hash(ValDiEdge(1, 2, (a=3, b=4)), UInt(0))

    @test hash(ValEdge(1, 2, (a=3, b=4))) != hash(ValEdge(1, 2))
    @test hash(ValDiEdge(1, 2, (a=3, b=4))) != hash(ValDiEdge(1, 2))
    @test hash(ValDiEdge(1, 2, (a=3, b=4))) != hash(ValDiEdge(2, 1, (a=3, b=4)))
    @test hash(ValEdge(1, 2)) != hash(ValEdge(1, 3))
    @test hash(ValDiEdge(1, 2)) != hash(ValDiEdge(1, 3))
    @test hash(ValEdge(1, 3)) != hash(ValEdge(2, 3))
    @test hash(ValDiEdge(1, 3)) != hash(ValDiEdge(2, 3))
    @test hash(ValEdge{Int8}(1, 2, (a=3, b=4))) == hash(ValEdge{Int16}(1, 2, (a=3, b=4)))
    @test hash(ValDiEdge{Int8}(1, 3, (a=3, b=4))) != hash(ValDiEdge{Int16}(1, 2, (a=3, b=4)))

end


text:
push!(LOAD_PATH,"src/") # necessary for unregistered packages

using Documenter
using FluxDiffUtils

makedocs(
    sitename = "FluxDiffUtils.jl",
    repo = "https://github.com/jlchan/FluxDiffUtils.jl",
    modules=[FluxDiffUtils],
    pages = [
        "Home" => "index.md",
    ]
)

deploydocs(
    repo = "github.com/jlchan/FluxDiffUtils.jl",
)


text:
### A Pluto.jl notebook ###
# v0.14.2

using Markdown
using InteractiveUtils

# ╔═╡ 794a5356-6e05-44b4-ad9d-375a1fa3613b
begin
	using PlutoUI
	using DataFrames
	using Statistics
	using CategoricalArrays
	using Test
end

# ╔═╡ b3dcf3f4-9d68-11eb-01a6-317fe0e0cd7e
md"""
## Introduction to DataFrames / 3 - Missing values

ref. [Introduction to DataFrames, Bogumił Kamiński](https://github.com/bkamins/Julia-DataFrames-Tutorial/blob/master/03_missingvalues.ipynb)

$(html"
© Pascal, April 2021
")
"""

# ╔═╡ d6993a5b-8158-4ce3-accb-8e8c514229b3
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ de97ded0-3952-4188-9a43-a8f383e1ae65
md"""
### Handling missing values

A `Julia` singleton type `Missing` allows us to deal with missing values:
"""

# ╔═╡ 624c9b0a-e7e8-4ec7-8802-dcf12a71117e
missing, typeof(missing)

# ╔═╡ 0921e011-b191-4b27-ac9e-17ea3bc0105f
md"""
Arrays automatically create an appropriate union type.
"""

# ╔═╡ 1fb36dcf-0b6a-4ab0-8c0f-0b2c7f23c3c9
begin
	x = [1, 2, missing, 3];
	typeof(x)
end

# ╔═╡ c64cd43a-8a17-4d2f-8d69-0db1073b7968
md"""
`ismissing` (predicate) function checks if a passed value is missing.
"""

# ╔═╡ 58aeff5c-3f12-4bbb-82b9-84bddb2159dd
ismissing.(x)   ## using broadcasting

# ╔═╡ 60f16698-3f4a-4a15-9366-10e7f801d9f5
md"""
We can extract the type combined with `Missing` from a `Union` via `nonmissingtype`

This is useful for arrays, see:
"""

# ╔═╡ 94ccb18d-4ea6-40f7-aa90-314d2c5704cc
eltype(x), nonmissingtype(eltype(x))

# ╔═╡ dcac6464-a916-4c57-9f8a-991b98baa4db
md"""
missing comparisons produce missing.
"""

# ╔═╡ 6e2620ac-736e-4531-b580-fd4bbf6112e7
missing == missing, missing != missing , missing < missing

# ╔═╡ a1f09f3b-e49d-49e2-8b82-aed03374b896
md"""
This is also true when missings are compared with values of other types.
"""

# ╔═╡ 123870b9-a79f-4fa7-b8b8-24e81a989ad9
1 == missing, π != missing , missing < √2

# ╔═╡ 92eccb60-fbf3-464f-a207-60dc2e23ebe6
md"""
`isequal`, `isless`, and `===` produce results of type `Bool`.

**Note that `missing` is considered greater than any numeric value.**
"""

# ╔═╡ 673e843d-0ce8-4cae-bc38-37aed282e4e6
isequal(missing, missing), missing === missing, isequal(1, missing), isless(1, missing)

# ╔═╡ cf1e4d30-2c35-49b9-a917-55df01522788
md"""
In the next few examples, we see that many (*not all!*) functions handle missing.
"""

# ╔═╡ ee9e58dd-2d41-4fc8-81f1-8b9c348df8bc
map(x -> x(missing), [sin, cos, zero, sqrt])   ## part 1

# ╔═╡ 1f71ca99-9b15-43e4-8afc-15ffc5d4a1a1
map(x -> x(missing, 1), [+, - , *, /, div])    ## part 2 - arithmetic ops

# ╔═╡ d6c88a1c-80b3-49c0-84af-716c28974f42
map(x -> x([1, 2, missing]),
	[minimum, maximum, extrema, mean, float])  ## part 3

# ╔═╡ fb6deb6c-c503-4a9d-8a4c-7f1c4c24ca76
md"""
`skipmissing` returns iterator skipping missing values. We can use `collect` and `skipmissing` to create an array that excludes these missing values."
"""

# ╔═╡ 9d7df124-4942-44c2-807b-f5c39d14ec14
skipmissing([1, missing, 2, missing]) |> collect

# ╔═╡ e8ffd09a-1944-4d0f-926d-0fb1f9dda7a4
md"""
Here we use `replace` to create a new array that replaces all missing values with some value (`NaN` in this case).
"""

# ╔═╡ a108cd54-45f7-48ec-a130-441dc7cf0f5e
replace([1.0, missing, 2.0, missing], missing => NaN)

# ╔═╡ 4049c78a-2250-42f1-b548-6fd0925711bf
md"""
Another way to do this:
"""

# ╔═╡ aeacb14c-db6f-4363-948c-00588ccb4583
coalesce.([1.0, missing, 2.0, missing], NaN)

# ╔═╡ 7e385236-ff48-4f6d-af46-924c727c25d7
md"""
We can also use `recode` from `CategoricalArrays.jl` if we have a default output value.
"""

# ╔═╡ 8fa0b5c6-22a4-46ba-ac06-4c1c71a0304a
recode([1.0, missing, 2.0, missing], false, missing => true)

# ╔═╡ e9eb3ad3-9278-4d5c-8de8-4821939080c6
md"""
There are also `replace!` and `recode!` functions that work in place.

Here is an example how you can to missing inputation in a data frame.
"""

# ╔═╡ 56a979fe-dc5e-4416-8f85-47f1f9db0285
df = DataFrame(a=[1, 2, missing], b=["a", "b", missing])

# ╔═╡ d32c9670-c8f9-46a5-ac54-4f9639325a5c
md"""
we change `df.a` vector in place:
"""

# ╔═╡ 9c8574fe-c181-4017-aada-00508d8360e3
replace!(df.a, missing => 999)

# ╔═╡ cf55fe7e-c5b9-46f0-815b-785f1fc64a95
md"""
Now we overwrite `df.b` with a new vector, because the replacement type is different than what `eltype(df.b)` accepts:
"""

# ╔═╡ 9f1fb12a-39bb-479a-b115-5364c7ea999d
df.b = coalesce.(df.b, 100)

# ╔═╡ 53c20047-c44d-44b9-93a6-9b093c326c0b
df

# ╔═╡ 38f2f006-3d80-4190-bbf8-e87c8400dbb4
md"""
We can use `unique` or `levels` to get unique values with or without missings, respectively.
"""

# ╔═╡ 2284ff98-9303-4d98-9409-34b2ac28117d
unique([1, missing, 2, missing]), levels([1, missing, 2, missing])

# ╔═╡ 47687957-622a-4ca4-9b45-361b011d4af9
md"""
In this next example, we convert x to y with `allowmissing`, where y has a type that accepts missings.
"""

# ╔═╡ 82a421c5-2b80-44d5-9bfb-242877ea624d
begin
	x₁ = [1, 2, 3]
	y₁ = allowmissing(x₁)
	typeof(y₁)
end

# ╔═╡ 7080eca8-faba-4c70-a7ba-a843d3ccb459
md"""
Then, we convert back with `disallowmissing`. This would fail if `y` contained missing values!
"""

# ╔═╡ 29010c54-cd2f-41e8-81fc-5354e4ca7f6e
begin
	z₁ = disallowmissing(y₁)
	x₁, y₁, z₁
end

# ╔═╡ 7021180c-01a1-4565-be9d-53e72ffee9bc
md"""
`disallowmissing` has an `error` keyword argument that can be used to decide how it should behave when it encounters a column that actually contains a missing value.
"""

# ╔═╡ ae893d9e-9ea5-4ea9-94d8-d646b8699088
df₁ = allowmissing(DataFrame(ones(2, 3), :auto))

# ╔═╡ 8d7b3298-8f14-4d21-a569-04c61ade9bda
df₁[1, 1] = missing

# ╔═╡ 9f8c6697-33bf-4362-96b1-98ab5802dd9e
df₁

# ╔═╡ e1a6c188-2299-494e-92a7-ab4556a59c85
@test_throws MethodError disallowmissing(df₁)    ## an error is thrown

# ╔═╡ 140fe595-4d87-4282-b880-b7d73dd2ecdb
disallowmissing(df₁, error=false) # column :x1 is left untouched as it contains missing

# ╔═╡ 1dd4fba9-c80d-44c0-8a6d-ff88a6dd9545
md"""
In this next example, we show that the type of each column in x is initially `Int64`. After using `allowmissing!` to accept missing values in columns 1 and 3, the types of those columns become `Union{Int64,Missing}`.
"""

# ╔═╡ 3bcca454-4077-4714-b6c9-a2c50f0a7748
with_terminal() do
	x₂ = DataFrame(rand(Int, 2,3), :auto)
	println("Before: ", eltype.(eachcol(x₂)))

	allowmissing!(x₂, 1)       # make first column accept missings
	allowmissing!(x₂, :x3)     # make :x3 column accept missings

	println("After: ", eltype.(eachcol(x₂)))
end

# ╔═╡ 5d4237f4-8bcc-47cf-a21d-8fe3c56e5ec1
md"""
In this next example, we will use `completecases` to find all the rows of a DataFrame that have complete data.
"""

# ╔═╡ 6a3921f9-ceb7-47e2-a042-574dc01d7951
with_terminal() do
	x₂ = DataFrame(A=[1, missing, 3, 4], B=["A", "B", missing, "C"])
	println("Complete cases:\n", completecases(x₂))
end

# ╔═╡ 99efee9e-d7b5-4bdb-ad60-f6613545bbd9
md"""
We can use `dropmissing` or `dropmissing!` to remove the rows with incomplete data from a DataFrame and either create a new DataFrame or mutate the original in-place.
"""

# ╔═╡ bbffa430-3ae8-42f8-80cc-04c646c3ae89
begin
	x₂ = DataFrame(A=[1, missing, 3, 4], B=["A", "B", missing, "C"])
	y₂ = dropmissing(x₂)
	dropmissing!(x₂), describe(x₂)
end

# ╔═╡ 5fe7ddc1-b36d-4757-a213-d7a5a437ef14
md"""
Alternatively you can pass `disallowmissing` keyword argument to `dropmissing` and `dropmissing!`
"""

# ╔═╡ 5bb6d14f-af90-4a46-ab25-168cf7e57670
x₃ = DataFrame(A=[1, missing, 3, 4], B=["A", "B", missing, "C"])

# ╔═╡ 47bb179b-f865-4d48-8e4f-317af7fbf590
dropmissing!(x₃, disallowmissing=false)

# ╔═╡ ca9927e5-3ff2-4c73-8206-4d7495668312
md"""
### Making functions missing-aware

If we have a function that does not handle missing values we can wrap it using `passmissing` function so that if any of its positional arguments is missing we will get a missing value in return.

In the example below we change how string function behaves:
"""

# ╔═╡ 13e1f76a-baac-4a79-b150-3be8af1d2d54
string(missing)

# ╔═╡ a6d55bb4-8c4e-4eee-8452-df6510d3dcc4
string(missing, " ", missing)

# ╔═╡ e6d047c6-40d4-4bb4-b566-1878d59f3172
string(1, 2, 3)

# ╔═╡ 8c7187f7-ad7d-47e2-8152-26969ed66f58
lift_string = passmissing(string)

# ╔═╡ d8a357b2-a6f3-4aa4-978b-b7d5967efae9
lift_string(missing)  ## No longer string "missing" but instead the value missing

# ╔═╡ 1d4e483c-4217-49fa-bee1-6c4f729dd823
lift_string(missing, " ", missing)

# ╔═╡ adc9610d-01f6-4060-97d2-9b3a4d0970d5
lift_string(1, 2, 3)

# ╔═╡ a4615e9b-d3d9-443a-9146-e411f9bbd749
md"""
### Aggregating rows containing missing values

Create an example DataFrame containing missing values:
"""

# ╔═╡ 19982031-e9a0-4a17-aad2-18d816096f45
df₂ = DataFrame(a=[1, missing, missing], b=[1, 2, missing])

# ╔═╡ 609834b7-c4fc-4f5e-9ec3-f442e084c94c
md"""
If we just run `sum` on the rows we get two missing entries:
"""

# ╔═╡ a27323a5-d651-41d0-a588-5d5fd9395522
sum.(eachrow(df₂))

# ╔═╡ 9c886a50-b929-4ef5-bbae-a727d7b691b1
md"""
We can apply `skipmissing` on the rows to avoid this problem
"""

# ╔═╡ 1fc8067a-da59-4e6b-a8e6-7e32fcad49dd
@test_throws ArgumentError sum.(skipmissing.(eachrow(df₂)))  ## However this triggers another problem

# ╔═╡ 372b242d-81eb-4b99-84e1-7cc79f244a1c
md"""
We get an error. The problem is that the last row of `df₂` contains only missing values, and since eachrow is type unstable the `eltype` of the result of `skipmissing` is unknown (so it is marked `Any`), see:
"""

# ╔═╡ e6a28a9e-942c-407d-9b99-da9bf1e594ff
skipmissing(eachrow(df₂)[end]) |> collect

# ╔═╡ 85063cfc-9992-4530-bf15-5d101b64e170
md"""
In such cases it is useful to switch to `Tables.namedtupleiterator` which is type stable as discussed in `01_constructors.ipynb` notebook.
"""

# ╔═╡ 77330f55-8ec6-492b-9e08-00dfcbb004eb
Tables.namedtupleiterator(df₂) .|>
	skipmissing .|>
    sum

# ≡ sum.(skipmissing.(Tables.namedtupleiterator(df)))

# ╔═╡ c54cdbf1-7787-4f0f-b9ff-d4389f4a7c5e
html"""

"""

# ╔═╡ Cell order:
# ╟─b3dcf3f4-9d68-11eb-01a6-317fe0e0cd7e
# ╠═794a5356-6e05-44b4-ad9d-375a1fa3613b
# ╟─d6993a5b-8158-4ce3-accb-8e8c514229b3
# ╟─de97ded0-3952-4188-9a43-a8f383e1ae65
# ╠═624c9b0a-e7e8-4ec7-8802-dcf12a71117e
# ╟─0921e011-b191-4b27-ac9e-17ea3bc0105f
# ╠═1fb36dcf-0b6a-4ab0-8c0f-0b2c7f23c3c9
# ╟─c64cd43a-8a17-4d2f-8d69-0db1073b7968
# ╠═58aeff5c-3f12-4bbb-82b9-84bddb2159dd
# ╟─60f16698-3f4a-4a15-9366-10e7f801d9f5
# ╠═94ccb18d-4ea6-40f7-aa90-314d2c5704cc
# ╟─dcac6464-a916-4c57-9f8a-991b98baa4db
# ╠═6e2620ac-736e-4531-b580-fd4bbf6112e7
# ╟─a1f09f3b-e49d-49e2-8b82-aed03374b896
# ╠═123870b9-a79f-4fa7-b8b8-24e81a989ad9
# ╟─92eccb60-fbf3-464f-a207-60dc2e23ebe6
# ╠═673e843d-0ce8-4cae-bc38-37aed282e4e6
# ╟─cf1e4d30-2c35-49b9-a917-55df01522788
# ╠═ee9e58dd-2d41-4fc8-81f1-8b9c348df8bc
# ╠═1f71ca99-9b15-43e4-8afc-15ffc5d4a1a1
# ╠═d6c88a1c-80b3-49c0-84af-716c28974f42
# ╟─fb6deb6c-c503-4a9d-8a4c-7f1c4c24ca76
# ╠═9d7df124-4942-44c2-807b-f5c39d14ec14
# ╟─e8ffd09a-1944-4d0f-926d-0fb1f9dda7a4
# ╠═a108cd54-45f7-48ec-a130-441dc7cf0f5e
# ╟─4049c78a-2250-42f1-b548-6fd0925711bf
# ╠═aeacb14c-db6f-4363-948c-00588ccb4583
# ╟─7e385236-ff48-4f6d-af46-924c727c25d7
# ╠═8fa0b5c6-22a4-46ba-ac06-4c1c71a0304a
# ╟─e9eb3ad3-9278-4d5c-8de8-4821939080c6
# ╠═56a979fe-dc5e-4416-8f85-47f1f9db0285
# ╟─d32c9670-c8f9-46a5-ac54-4f9639325a5c
# ╠═9c8574fe-c181-4017-aada-00508d8360e3
# ╟─cf55fe7e-c5b9-46f0-815b-785f1fc64a95
# ╠═9f1fb12a-39bb-479a-b115-5364c7ea999d
# ╠═53c20047-c44d-44b9-93a6-9b093c326c0b
# ╟─38f2f006-3d80-4190-bbf8-e87c8400dbb4
# ╠═2284ff98-9303-4d98-9409-34b2ac28117d
# ╟─47687957-622a-4ca4-9b45-361b011d4af9
# ╠═82a421c5-2b80-44d5-9bfb-242877ea624d
# ╟─7080eca8-faba-4c70-a7ba-a843d3ccb459
# ╠═29010c54-cd2f-41e8-81fc-5354e4ca7f6e
# ╟─7021180c-01a1-4565-be9d-53e72ffee9bc
# ╠═ae893d9e-9ea5-4ea9-94d8-d646b8699088
# ╠═8d7b3298-8f14-4d21-a569-04c61ade9bda
# ╠═9f8c6697-33bf-4362-96b1-98ab5802dd9e
# ╠═e1a6c188-2299-494e-92a7-ab4556a59c85
# ╠═140fe595-4d87-4282-b880-b7d73dd2ecdb
# ╟─1dd4fba9-c80d-44c0-8a6d-ff88a6dd9545
# ╠═3bcca454-4077-4714-b6c9-a2c50f0a7748
# ╟─5d4237f4-8bcc-47cf-a21d-8fe3c56e5ec1
# ╠═6a3921f9-ceb7-47e2-a042-574dc01d7951
# ╟─99efee9e-d7b5-4bdb-ad60-f6613545bbd9
# ╠═bbffa430-3ae8-42f8-80cc-04c646c3ae89
# ╟─5fe7ddc1-b36d-4757-a213-d7a5a437ef14
# ╠═5bb6d14f-af90-4a46-ab25-168cf7e57670
# ╠═47bb179b-f865-4d48-8e4f-317af7fbf590
# ╟─ca9927e5-3ff2-4c73-8206-4d7495668312
# ╠═13e1f76a-baac-4a79-b150-3be8af1d2d54
# ╠═a6d55bb4-8c4e-4eee-8452-df6510d3dcc4
# ╠═e6d047c6-40d4-4bb4-b566-1878d59f3172
# ╠═8c7187f7-ad7d-47e2-8152-26969ed66f58
# ╠═d8a357b2-a6f3-4aa4-978b-b7d5967efae9
# ╠═1d4e483c-4217-49fa-bee1-6c4f729dd823
# ╠═adc9610d-01f6-4060-97d2-9b3a4d0970d5
# ╟─a4615e9b-d3d9-443a-9146-e411f9bbd749
# ╠═19982031-e9a0-4a17-aad2-18d816096f45
# ╟─609834b7-c4fc-4f5e-9ec3-f442e084c94c
# ╠═a27323a5-d651-41d0-a588-5d5fd9395522
# ╟─9c886a50-b929-4ef5-bbae-a727d7b691b1
# ╠═1fc8067a-da59-4e6b-a8e6-7e32fcad49dd
# ╟─372b242d-81eb-4b99-84e1-7cc79f244a1c
# ╠═e6a28a9e-942c-407d-9b99-da9bf1e594ff
# ╟─85063cfc-9992-4530-bf15-5d101b64e170
# ╠═77330f55-8ec6-492b-9e08-00dfcbb004eb
# ╟─c54cdbf1-7787-4f0f-b9ff-d4389f4a7c5e

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
mutable struct ILUTPreconditioner{Tv, Ti} <: AbstractPreconditioner{Tv,Ti}
    A::ExtendableSparseMatrix{Tv,Ti}
    fact::IncompleteLU.ILUFactorization{Tv,Ti}
    droptol::Float64
    function ILUTPreconditioner{Tv,Ti}(;droptol=1.0e-3) where {Tv,Ti}
        p=new()
        p.droptol=droptol
        p
    end
end

"""
```
ILUTPreconditioner(;droptol=1.0e-3,valuetype=Float64, indextype=Int64)
ILUTPreconditioner(matrix; droptol=1.0e-3)
```

Create the [`ILUTPreconditioner`](@ref) wrapping the one
from [IncompleteLU.jl](https://github.com/haampie/IncompleteLU.jl)
For using this, you need to issue `using IncompleteLU`.
"""
ILUTPreconditioner(;valuetype::Type=Float64, indextype::Type=Int64, kwargs...)=ILUTPreconditioner{valuetype,indextype}(;kwargs...)

@eval begin
    @makefrommatrix ILUTPreconditioner
end

function update!(precon::ILUTPreconditioner)
    A=precon.A
    @inbounds flush!(A)
    precon.fact=IncompleteLU.ilu(A.cscmatrix,τ=precon.droptol)
end


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 4a425262-509a-11eb-08ac-cd5457e6c813
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
			"CSV",
			"DataFrames",
			"Plots",
			"PlutoUI",
	])
	using CSV
	using DataFrames
	using Dates
	using InteractiveUtils
	using Markdown
	using Statistics
	using Plots
	using PlutoUI
end

# ╔═╡ 6c0d01d0-509a-11eb-2179-cbd8c3ec59d2
md"""
# International COVID-19 Visualization
This is a visualization of weekly reported COVID-19 cases against total reported cases. The plot is done on a log-scale causing exponential growth to be represented linearly. Since time is not explicitly represented on the graph, direct comparisons can be drawn between countries who experienced similar case trends at different times.

This project was created by Jeremiah DeGreeff with inspiration from Problem Set 1 for [MIT Course 6.S083 Spring 2020](https://www.eecs.mit.edu/academics-admissions/academic-information/subject-updates-spring-2020/6s083) and [this visualization](https://aatishb.com/covidtrends/) created by Aatish Bhatia.
"""

# ╔═╡ bf71b0fa-509a-11eb-36d7-9fb6af718625
md"""
## Data Preparation
"""

# ╔═╡ cc72d7a2-509a-11eb-3515-b70d385cd79e
begin
	data_url = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"
	md"""
	COVID-19 Cases Data Source: [JHU CSSE](https://github.com/CSSEGISandData/COVID-19)
	"""
end

# ╔═╡ 0f6c6f82-509b-11eb-21ae-bb984f52b4d4
begin
	filename = "data/covid_country_data.csv"
	download(data_url, filename)
	all_data = CSV.read(filename, DataFrame)
	rename!(all_data, 1 => "province", 2 => "country")

	# Enumerate all countries
	all_countries = collect(all_data[:, 2])
	unique_countries = unique(all_countries)

	# Extract the dates we want to plot
	format = Dates.DateFormat("m/d/Y")
	date_strings = String.(names(all_data)[5:end])
	dates = parse.(Date, date_strings, format) .+ Year(2000)

	all_data
end

# ╔═╡ c31e21ce-509b-11eb-3ed5-ab97c8adb545
begin
	date_slider = @bind date_index Slider(1:length(dates), default=1, show_value=true)
	md"""
	Days since $(dates[1]): $(date_slider)
	"""
end

# ╔═╡ d2af04be-509b-11eb-2b1b-65d6b02ba7db
begin
	double_rate_slider = @bind double_rate Slider(1:100, default=7, show_value=true)
	md"""
	Exponential growth doubling rate (in days): $(double_rate_slider)
	"""
end

# ╔═╡ da696e06-509b-11eb-2a3c-8fc8c4fb3039
begin
	legend_checkbox = @bind legend CheckBox(default=true)
	md"""
	Show legend $(legend_checkbox)
	"""
end

# ╔═╡ df9860c6-509b-11eb-307c-7bd5236d725c
md"""
## Plotting Library
"""

# ╔═╡ e64c39c4-509b-11eb-361c-771e9909c07b
# Function to apply a function over a dictionary of data sequences
apply_over_dict(func, dict, args=()) = Dict(k => isempty(args) ? func(v) : func(v, args) for (k,v) ∈ dict)

# ╔═╡ edb8991e-509b-11eb-24ea-4dca6c092374
# Function to convert aggregate data to weekly data
aggreggate_to_weekly_change(aggregate_data) = [aggregate_data[i] -  aggregate_data[max(1, i - 6)] for i ∈ 1:length(aggregate_data)]

# ╔═╡ f81bc298-509b-11eb-2a44-4f46c85811c1
# Function to prepare version of data for a log range
truncate(data, min_value = 0) = map(x -> x > min_value ? x : NaN, data)

# ╔═╡ fd74fe2e-509b-11eb-019b-493a9a1eb249
# Function to determine upper bound for a log range
max_value(data) = 10^ceil(Int64, log(10, maximum([d for (_, s) ∈ data for d ∈ s if !isnan(d)])))

# ╔═╡ 04546f7c-509c-11eb-361d-3dac17d55476
# Function to plot weekly change vs aggregate on a log scale
function plot_covid(aggregate_data, dates, date_index; min_value = 10, double_rate = missing, normalize = false)
	p = plot(xscale=:log10, yscale=:log10, leg=:bottomright)
	title!("Trajectory of COVID-19 Confirmed Cases ($(dates[date_index]))")
	xlabel!(normalize ? "Total Confirmed Cases per 100,000 Residents" : "Total Confirmed Cases")
	ylabel!(normalize ? "Recent Confirmed Cases per 100,000 Residents" : "Recent Confirmed Cases")

    # Weekly change data is the number of new cases in the past week
    weekly_change_data = apply_over_dict(aggreggate_to_weekly_change, aggregate_data)

	aggregate_data = apply_over_dict(truncate, aggregate_data, min_value)
    weekly_change_data = apply_over_dict(truncate, weekly_change_data)

	min_x = 10^floor(Int64, log(10, min_value))
	if all(s -> all(isnan, s), values(aggregate_data))
		x_range = (min_x, min_x * 10)
		y_range = (1, 10)
	else
		x_range = (min_x, max_value(aggregate_data))
		y_range = (1, max_value(weekly_change_data))
	end
	xlims!(x_range)
	ylims!(y_range)

	for (i, c) in enumerate(keys(aggregate_data))
		plot!(aggregate_data[c][1:date_index], weekly_change_data[c][1:date_index], color=i, label=false)
		if !isnan(aggregate_data[c][date_index]) && !isnan(weekly_change_data[c][date_index])
			scatter!(aggregate_data[c][date_index:date_index], weekly_change_data[c][date_index:date_index], color=i, label=c, markerstrokewidth=0)
		end
	end

	if !ismissing(double_rate)
		# Exponential data provides a reference for what pure exponential growth with a particular growth rate would look like
		exponential_aggregate = [2^(t / double_rate) for t ∈ 0:ceil(log2(x_range[2]) * double_rate)]
		exponential_weekly_change = aggreggate_to_weekly_change(exponential_aggregate)

		exponential_aggregate = truncate(exponential_aggregate)
		exponential_weekly_change = truncate(exponential_weekly_change)

		plot!(exponential_aggregate, exponential_weekly_change, color=:gray, lw=2, label="Exponential Growth (d=$(double_rate))")
	end
	p
end

# ╔═╡ 7639ff12-509c-11eb-3215-b108ab8d04af
md"""
## PlutoUI Library Overrides
"""

# ╔═╡ 782fc4be-509c-11eb-0e53-bbad102c817d
begin
	import Markdown: htmlesc, withtag

    struct MultiSelect
        options::Array{Pair{<:AbstractString,<:Any},1}
        default::Union{Missing, AbstractVector{AbstractString}}
        size::Int
    end

	MultiSelect(options::Array{<:AbstractString,1}; default=missing, size=5) = MultiSelect([o => o for o in options], default, size)
    MultiSelect(options::Array{<:Pair{<:AbstractString,<:Any},1}; default=missing, size=5) = MultiSelect(options, default, size)

	function Base.show(io::IO, ::MIME"text/html", select::MultiSelect)
        withtag(io, Symbol("select multiple"), :size => string(select.size)) do
            for o in select.options
                print(io, """<option value="$(htmlesc(o.first))"$(!ismissing(select.default) && o.first ∈ select.default ? " selected" : "")>""")
                if showable(MIME"text/html"(), o.second)
                    show(io, MIME"text/html"(), o.second)
                else
                    print(io, o.second)
                end
                print(io, "</option>")
            end
        end
    end

	Base.get(select::MultiSelect) = ismissing(select.default) ? Any[] : select.default
end

# ╔═╡ 643690d0-509c-11eb-14e0-3ff522e70c9b
begin
	select = @bind countries MultiSelect(unique_countries, size=length(unique_countries))
    md"""
    Select the countries to display.

    $(select)
    """
end

# ╔═╡ 7cca2120-509b-11eb-17a5-e36d54879f37
begin
	if !isempty(countries)
aggregate_data = Dict(c => collect(mapcols(sum, all_data[all_countries.==c, :][:, 5:end])[1, :]) for c ∈ countries)
		p = plot_covid(aggregate_data, dates, date_index; min_value=100, double_rate=double_rate)
		if !legend plot!(p, leg=false) end
		p |> as_svg
	end
end

# ╔═╡ Cell order:
# ╟─4a425262-509a-11eb-08ac-cd5457e6c813
# ╟─6c0d01d0-509a-11eb-2179-cbd8c3ec59d2
# ╟─bf71b0fa-509a-11eb-36d7-9fb6af718625
# ╟─cc72d7a2-509a-11eb-3515-b70d385cd79e
# ╟─0f6c6f82-509b-11eb-21ae-bb984f52b4d4
# ╟─643690d0-509c-11eb-14e0-3ff522e70c9b
# ╟─c31e21ce-509b-11eb-3ed5-ab97c8adb545
# ╟─d2af04be-509b-11eb-2b1b-65d6b02ba7db
# ╟─da696e06-509b-11eb-2a3c-8fc8c4fb3039
# ╟─7cca2120-509b-11eb-17a5-e36d54879f37
# ╟─df9860c6-509b-11eb-307c-7bd5236d725c
# ╟─e64c39c4-509b-11eb-361c-771e9909c07b
# ╟─edb8991e-509b-11eb-24ea-4dca6c092374
# ╟─f81bc298-509b-11eb-2a44-4f46c85811c1
# ╟─fd74fe2e-509b-11eb-019b-493a9a1eb249
# ╟─04546f7c-509c-11eb-361d-3dac17d55476
# ╟─7639ff12-509c-11eb-3215-b108ab8d04af
# ╟─782fc4be-509c-11eb-0e53-bbad102c817d

max_stars_repo_name: jrdegreeff/COVID-Visualization

max_issues_repo_name: jrdegreeff/COVID-Visualization

max_forks_repo_name: jrdegreeff/COVID-Visualization


text:
# This file was generated, do not modify it. # hide
r = report(m)

figure(figsize=(8,6))

vals_sf = r.parameter_values[:, 1]
vals_bf = r.parameter_values[:, 2]

tricontourf(vals_sf, vals_bf, r.measurements)
xticks(1:3:12, fontsize=12)
xlabel("Number of sub-features", fontsize=14)
yticks(0.4:0.2:1, fontsize=12)
ylabel("Bagging fraction", fontsize=14)

savefig("assets/literate/A-ensembles-2-heatmap.svg") # hide


text:
@testset "Subpopulation" begin

@testset "Constructor" begin
    subpopulation = Subpopulation( "Subpop" )
    @test all( [subpopulation.name == "Subpop",
        subpopulation.sourceNode == "active",
        isempty( subpopulation.timeConds ),
        isempty( subpopulation.historyConds ),
        isempty( subpopulation.attributeConds )] )
end  # @testset "Constructor"

subpopulation = Subpopulation( "Subpop" )

@testset "function setSubpopulationSourceNode!" begin
    setSubpopulationSourceNode!( subpopulation, "Node" )
    @test subpopulation.sourceNode == "Node"
    setSubpopulationSourceNode!( subpopulation, "" )
    @test subpopulation.sourceNode == "active"
end  # @testset "function setSubpopulationSourceNode!"

cond1 = MPcondition( "Attribute", !=, "value")
cond2 = MPcondition( "stARted AS", ==, "Pilot" )
cond3 = MPcondition( "TeNuRe", <=, 36 )

@testset "function addSubpopulationCondition!" begin
    addSubpopulationCondition!( subpopulation, cond1, cond2, cond3 )
    @test all( length.( [subpopulation.timeConds, subpopulation.historyConds,
        subpopulation.attributeConds] ) .== 1 )
    @test ( subpopulation.historyConds[1].attribute == "started as" ) &&
        ( subpopulation.timeConds[1].attribute == "tenure" )
end  # @testset "function addSubpopulationCondition!"

@testset "function clearSubpopulationConditions!" begin
    clearSubpopulationConditions!( subpopulation )
    @test all( isempty.( [subpopulation.timeConds, subpopulation.historyConds,
        subpopulation.attributeConds] ) )
end  # @testset "function clearSubpopulationConditions!"

@testset "function setSubpopulationConditions!" begin
    setSubpopulationConditions!( subpopulation, [cond1, cond2, cond3] )
    @test all( length.( [subpopulation.timeConds, subpopulation.historyConds,
        subpopulation.attributeConds] ) .== 1 )
end  # @testset "function setSubpopulationConditions!"

end  # @testset "Subpopulation"


text:
### MCJob

abstract MCJob

abstract GibbsJob <: MCJob

# Set defaults for possibly unspecified output options

function augment_variable_outopts!(outopts::Dict)
  destination = get!(outopts, :destination, :nstate)

  if destination != :none
    if !haskey(outopts, :monitor)
      outopts[:monitor] = [:value]
    end

    if destination == :iostream
      if !haskey(outopts, :filepath)
        outopts[:filepath] = ""
      end

      if !haskey(outopts, :filesuffix)
        outopts[:filesuffix] = "csv"
      end

      if !haskey(outopts, :flush)
        outopts[:flush] = false
      end
    end
  end
end

augment_variable_outopts!{K, V}(outopts::Vector{Dict{K, V}}) = map(augment_variable_outopts!, outopts)

function augment_parameter_outopts!(outopts::Dict)
  augment_variable_outopts!(outopts)

  if outopts[:destination] != :none
    if !haskey(outopts, :diagnostics)
      outopts[:diagnostics] = Symbol[]
    end
  end
end

augment_parameter_outopts!{K, V}(outopts::Vector{Dict{K, V}}) = map(augment_parameter_outopts!, outopts)

# initialize_output() needs to be defined for custom variable state or NState input arguments
# Thus multiple dispatch allows to extend the code base to accommodate new variable states or NStates

function initialize_output(state::BasicUnvVariableState, n::Integer, outopts::Dict)
  output::Union{VariableNState, VariableIOStream, Void}

  if outopts[:destination] == :nstate
    output = BasicUnvVariableNState(n, eltype(state))
  elseif outopts[:destination] == :iostream
    output = BasicVariableIOStream(
      (),
      n,
      filepath=outopts[:filepath],
      filesuffix=outopts[:filesuffix],
      mode="w"
    )
  elseif outopts[:destination] == :none
    output = nothing
  else
    error(":destination must be set to :nstate or :iostream or :none, got $(outopts[:destination])")
  end

  output
end

function initialize_output(state::BasicMuvVariableState, n::Integer, outopts::Dict)
  output::Union{VariableNState, VariableIOStream, Void}

  if outopts[:destination] == :nstate
    output = BasicMuvVariableNState(state.size, n, eltype(state))
  elseif outopts[:destination] == :iostream
    output = BasicVariableIOStream(
      (state.size,),
      n,
      filepath=outopts[:filepath],
      filesuffix=outopts[:filesuffix],
      mode="w"
    )
  elseif outopts[:destination] == :none
    output = nothing
  else
    error(":destination must be set to :nstate or :iostream or :none, got $(outopts[:destination])")
  end

  output
end

function initialize_output(state::BasicMavVariableState, n::Integer, outopts::Dict)
  output::Union{VariableNState, VariableIOStream, Void}

  if outopts[:destination] == :nstate
    output = BasicMavVariableNState(state.size, n, eltype(state))
  elseif outopts[:destination] == :iostream
    output = BasicVariableIOStream(
      state.size,
      n,
      filepath=outopts[:filepath],
      filesuffix=outopts[:filesuffix],
      mode="w"
    )
  elseif outopts[:destination] == :none
    output = nothing
  else
    error(":destination must be set to :nstate or :iostream or :none, got $(outopts[:destination])")
  end

  output
end

function initialize_output(state::BasicDiscUnvParameterState, n::Integer, outopts::Dict)
  output::Union{VariableNState, VariableIOStream, Void}

  if outopts[:destination] == :nstate
    output = BasicDiscUnvParameterNState(n, outopts[:monitor], outopts[:diagnostics], eltype(state)...)
  elseif outopts[:destination] == :iostream
    output = BasicDiscParamIOStream(
      (),
      n,
      outopts[:monitor],
      filepath=outopts[:filepath],
      filesuffix=outopts[:filesuffix],
      diagnostickeys=outopts[:diagnostics],
      mode="w"
    )
  elseif outopts[:destination] == :none
    output = nothing
  else
    error(":destination must be set to :nstate or :iostream or :none, got $(outopts[:destination])")
  end

  output
end

function initialize_output(state::BasicDiscMuvParameterState, n::Integer, outopts::Dict)
  output::Union{VariableNState, VariableIOStream, Void}

  if outopts[:destination] == :nstate
    output = BasicDiscMuvParameterNState(state.size, n, outopts[:monitor], outopts[:diagnostics], eltype(state)...)
  elseif outopts[:destination] == :iostream
    output = BasicDiscParamIOStream(
      (state.size,),
      n,
      outopts[:monitor],
      filepath=outopts[:filepath],
      filesuffix=outopts[:filesuffix],
      diagnostickeys=outopts[:diagnostics],
      mode="w"
    )
    mark(output)
  elseif outopts[:destination] == :none
    output = nothing
  else
    error(":destination must be set to :nstate or :iostream or :none, got $(outopts[:destination])")
  end

  output
end

function initialize_output(state::BasicContUnvParameterState, n::Integer, outopts::Dict)
  output::Union{VariableNState, VariableIOStream, Void}

  if outopts[:destination] == :nstate
    output = BasicContUnvParameterNState(n, outopts[:monitor], outopts[:diagnostics], eltype(state))
  elseif outopts[:destination] == :iostream
    output = BasicContParamIOStream(
      (),
      n,
      outopts[:monitor],
      filepath=outopts[:filepath],
      filesuffix=outopts[:filesuffix],
      diagnostickeys=outopts[:diagnostics],
      mode="w"
    )
  elseif outopts[:destination] == :none
    output = nothing
  else
    error(":destination must be set to :nstate or :iostream or :none, got $(outopts[:destination])")
  end

  output
end

function initialize_output(state::BasicContMuvParameterState, n::Integer, outopts::Dict)
  output::Union{VariableNState, VariableIOStream, Void}

  if outopts[:destination] == :nstate
    output = BasicContMuvParameterNState(state.size, n, outopts[:monitor], outopts[:diagnostics], eltype(state))
  elseif outopts[:destination] == :iostream
    output = BasicContParamIOStream(
      (state.size,),
      n,
      outopts[:monitor],
      filepath=outopts[:filepath],
      filesuffix=outopts[:filesuffix],
      diagnostickeys=outopts[:diagnostics],
      mode="w"
    )
    mark(output)
  elseif outopts[:destination] == :none
    output = nothing
  else
    error(":destination must be set to :nstate or :iostream or :none, got $(outopts[:destination])")
  end

  output
end

function initialize_task!(job::MCJob)
  # Hook inside task to allow remote resetting
  task_local_storage(:reset, job.resetplain!)

  while true
    iterate(job)
  end
end

Base.close(job::MCJob) = job.close(job)

Base.reset(job::MCJob) = job.reset!(job)
Base.reset(job::MCJob, x::Real) = job.reset!(job, x)
Base.reset{N<:Real}(job::MCJob, x::Vector{N}) = job.reset!(job, x)

save(job::MCJob, i::Integer) = job.save!(job, i)

iterate(job::MCJob) = job.iterate!(job)

Base.run(job::MCJob) = job.run!(job)
Base.run{J<:MCJob}(job::Vector{J}) = map(run, job)


text:
@deprecate get_vel_s velf(state).s
@deprecate get_vel_t velf(state).t
@deprecate get_name nothing false

@deprecate propagate(veh, state, roadway, Δt) propagate(veh, state, a, roadway, Δt)

function simulate!(
    scene::Frame{E},
    roadway::R,
    models::Dict{I,M},
    nticks::Int64,
    timestep::Float64;
    rng::AbstractRNG = Random.GLOBAL_RNG,
    callbacks = nothing
) where {E<:Entity,A,R,I,M<:DriverModel}
    Base.depwarn(
"`simulate!` without specifying a pre-allocated data structure for `scenes` is now deprecated.\n
You probably want to use the `simulate` function instead.\n
Alternatively, you can provide a pre-allocated data structure via the `scenes=` keyword",
        :simulate_no_prealloc
    )
    return simulate(scene, roadway, models, nticks, timestep, rng=rng, callbacks=callbacks)
end

"""
    # TODO: this should be removed, but where to document the function then?
    run_callback(callback::Any, scenes::Vector{F}, roadway::R, models::Dict{I,M}, tick::Int) where {F,I,R,M<:DriverModel}
    run_callback(callback::Any, rec::EntityQueueRecord{S,D,I}, roadway::R, models::Dict{I,M}, tick::Int) where {S,D,I,R,M<:DriverModel}
run callback and return whether simlation should terminate
A new method should be implemented when defining a new callback object.
"""
function run_callback(callback, scenes, actions::Union{Nothing, Vector{Frame{A}}}, roadway, models, tick) where {A<:EntityAction}
    Base.depwarn(
"Using a deprecated version of `run_callback`. Since v0.7.10, user-defined callback functions should also take an `actions` argument.
 If you have implemented `run_callback` with an actions argument, make sure the method signature is more specific than this one.\n
Check the section `Simulation > Callbacks` in the package documentation for more information
 on how to define callback functions.\n
Your function call is being forwarded to `run_callback` without actions argument.",
        :run_callback_no_actions
    )
    return run_callback(callback, scenes, roadway, models, tick)
end

# 1D state

state_1d_dep_msgs = """
    State1D, Vehicle1D, and Scene1D are deprecated, use VehicleState, Vehicle, and Scene instead.
"""

"""
    State1D

A data type to represent one dimensional states

# Fields
    - `s::Float64` position
    - `v::Float64` speed [m/s]
"""
struct State1D
    s::Float64 # position
    v::Float64 # speed [m/s]
end

function vel(s::State1D)
    Base.depwarn(state_1d_dep_msgs, :vel)
    return s.v
end

function Base.write(io::IO, ::MIME"text/plain", s::State1D)
    Base.depwarn(state_1d_dep_msgs, :write)
    @printf(io, "%.16e %.16e", s.s, s.v)
end
function Base.read(io::IO, ::MIME"text/plain", ::Type{State1D})
    Base.depwarn(state_1d_dep_msgs, :read)
    i = 0
    tokens = split(strip(readline(io)), ' ')
    s = parse(Float64, tokens[i+=1])
    v = parse(Float64, tokens[i+=1])
    return State1D(s,v)
end

"""
    Vehicle1D
A specific instance of the Entity type defined in Records.jl to represent vehicles in 1d environments.
"""
const Vehicle1D = Entity{State1D, VehicleDef, Int64}

"""
    Scene1D

A specific instance of the Frame type defined in Records.jl to represent a list of vehicles in 1d environments.

# constructors
    Scene1D(n::Int=100)
    Scene1D(arr::Vector{Vehicle1D})
"""
const Scene1D = Frame{Vehicle1D}

function Scene1D(n::Int=100)
    Base.depwarn(state_1d_dep_msgs, :Scene1D)
    Frame(Vehicle1D, n)
end

function Scene1D(arr::Vector{Vehicle1D})
    Base.depwarn(state_1d_dep_msgs, :Scene1D)
    Frame{Vehicle1D}(arr, length(arr))
end

function get_center(veh::Vehicle1D)
    Base.depwarn(state_1d_dep_msgs, :get_center)
    veh.state.s
end
function get_footpoint(veh::Vehicle1D)
    Base.depwarn(state_1d_dep_msgs, :get_footpoint)
    veh.state.s
end
function get_front(veh::Vehicle1D)
    Base.depwarn(state_1d_dep_msgs, :get_front)
    veh.state.s + length(veh.def)/2
end
function get_rear(veh::Vehicle1D)
    Base.depwarn(state_1d_dep_msgs, :get_rear)
    veh.state.s - length(veh.def)/2
end

function get_headway(veh_rear::Vehicle1D, veh_fore::Vehicle1D, roadway::StraightRoadway)
    Base.depwarn(state_1d_dep_msgs, :get_headway)
    return get_headway(get_front(veh_rear), get_rear(veh_fore), roadway)
end
function get_neighbor_fore(scene::Frame{Entity{State1D, D, I}}, vehicle_index::I, roadway::StraightRoadway) where {D, I}
    Base.depwarn(state_1d_dep_msgs, :get_neighbor_fore)
    ego = scene[vehicle_index]
    best_ind = nothing
    best_gap = Inf
    for (i,veh) in enumerate(scene)
        if i != vehicle_index
            Δs = get_headway(ego, veh, roadway)
            if Δs < best_gap
                best_gap, best_ind = Δs, i
            end
        end
    end
    return NeighborLongitudinalResult(best_ind, best_gap)
end
function get_neighbor_rear(scene::Frame{Entity{State1D, D, I}}, vehicle_index::I, roadway::StraightRoadway) where {D, I}
    Base.depwarn(state_1d_dep_msgs, :get_neighbor_rear)
    ego = scene[vehicle_index]
    best_ind = nothing
    best_gap = Inf
    for (i,veh) in enumerate(scene)
        if i != vehicle_index
            Δs = get_headway(veh, ego, roadway)
            if Δs < best_gap
                best_gap, best_ind = Δs, i
            end
        end
    end
    return NeighborLongitudinalResult(best_ind, best_gap)
end

function propagate(veh::Vehicle1D, action::LaneFollowingAccel, roadway::StraightRoadway, Δt::Float64)
    Base.depwarn(state_1d_dep_msgs, :propagate)
    a = action.a
    s, v = veh.state.s, veh.state.v

    s′ = s + v*Δt + a*Δt*Δt/2
    v′ = max(v + a*Δt, 0.)  # no negative velocities

    s′ = mod_position_to_roadway(s′, roadway)

    return State1D(s′, v′)
end

function observe!(model::LaneFollowingDriver, scene::Frame{Entity{State1D, D, I}}, roadway::StraightRoadway, egoid::I) where {D, I}
    Base.depwarn(state_1d_dep_msgs, :observe!)
    vehicle_index = findfirst(egoid, scene)

    fore_res = get_neighbor_fore(scene, vehicle_index, roadway)

    v_ego = vel(scene[vehicle_index].state)
    v_oth = vel(scene[fore_res.ind].state)
    headway = fore_res.Δs

    track_longitudinal!(model, v_ego, v_oth, headway)

    return model
end


text:
module Triplot

import GR

export tricontour, tricontourf, tripcolor

include("trimesh.jl")
include("tricontour.jl")
include("tricontourf.jl")
include("tripcolor.jl")

end


text:
using UnivariateDifferentiableOptimizationAlgorithms
using Base.Test
# write your own tests here

# include("testcompare-1D-6-variant.jl")


text:
import UnicodePlots
import UnicodePlots: BrailleCanvas, nrows
import UnicodePlots: lineplot!, scatterplot!

function render_unicode(cur_size, cur_data, cur_layout, cur_config)

  work_layout = deepcopy(cur_layout)

  cur_plot, has_legend, is_linear_x, is_linear_y, max_label_length = (
    _configure_repl_layout(cur_size, cur_data, work_layout)
  )

  if isempty(cur_data)
    return cur_plot
  end

  # repl doesn't allow legend placement
  @assert isempty(work_layout["legend"])
  delete!(work_layout, "legend")

  # no annotations in repl
  delete!(work_layout, "annotations")

  work_config = deepcopy(cur_config)
  delete!(work_config, "responsive")
  @assert isempty(work_config)

  for (cur_index, cur_datum) in enumerate(cur_data)
    @assert cur_datum["mode"] in ["lines", "markers"]

    work_datum = deepcopy(cur_datum)
    singular_mode = cur_datum["mode"][1:end-1]

    if cur_datum["mode"] == "lines"
      plotfunc = lineplot!
    else
      plotfunc = scatterplot!
    end

    cur_kwargs = Dict()
    cur_kwargs[:color] = Symbol(cur_datum[singular_mode]["color"])

    if haskey(cur_datum, "showlegend")
      @assert !cur_datum["showlegend"]
      @assert !haskey(cur_datum, "name")
      delete!(work_datum, "showlegend")
    elseif haskey(cur_datum, "name")
      cur_label = cur_datum["name"]
      if length(cur_label) > max_label_length
        cur_label = cur_label[1:(max_label_length-1)] * "…"
      end

      if has_legend
        cur_kwargs[:name] = cur_label
      end
      delete!(work_datum, "name")
    end

    cur_x = cur_datum["x"]
    cur_y = cur_datum["y"]

    is_linear_x || ( cur_x = log10.(float.(cur_x)) )
    is_linear_y || ( cur_y = log10.(float.(cur_y)) )

    plotfunc(
      cur_plot, cur_x, cur_y; cur_kwargs...
    )

    delete!(work_datum[singular_mode], "color")
    delete!(work_datum[singular_mode], "width")
    delete!(work_datum[singular_mode], "size")

    @assert isempty(work_datum[singular_mode])
    delete!(work_datum, singular_mode)

    delete!(work_datum, "mode")
    delete!(work_datum, "x")
    delete!(work_datum, "y")

    isempty(work_datum) || @show work_datum
    @assert isempty(work_datum)
  end

  @assert isempty(work_layout["shapes"])
  delete!(work_layout, "shapes")
  @assert isempty(work_layout)

  for (k,v) in work_layout
    @show (k,v)
  end

  cur_plot
end

function _configure_repl_layout(cur_size, cur_data, work_layout)
  size_x, size_y = cur_size

  has_legend = work_layout["showlegend"]
  delete!(work_layout, "showlegend")

  if isempty(cur_data)
    max_label_length = 0
  else
    max_label_length = maximum(
      length.(map(
        cur_datum -> haskey(cur_datum, "name") ? cur_datum["name"] : "", cur_data
      ))
    )

    max_label_length = min(16, max_label_length)

    if has_legend
      size_x -= max_label_length
    end
  end

  _get_lims(x) = collect(extrema(Base.Iterators.flatten(
    map(cur_datum -> cur_datum[x], cur_data)
  )))

  function _get_extrema(x)
    if haskey(work_layout["$(x)axis"], "range")
      @assert work_layout["$(x)axis"]["autorange"] == false
      cur_extrema = work_layout["$(x)axis"]["range"]
      delete!(work_layout["$(x)axis"], "range")
    else
      @assert work_layout["$(x)axis"]["autorange"] == true
      if isempty(cur_data)
        cur_extrema = (1, 10)
      else
        cur_extrema = _get_lims(x)
      end
    end

    delete!(work_layout["$(x)axis"], "autorange")
    return cur_extrema
  end

  if isempty(cur_data)
    min_x, max_x = _get_extrema("x")
    min_y, max_y = _get_extrema("y")
  else
    min_x, max_x = _get_extrema("x")
    min_y, max_y = _get_extrema("y")
  end

  function _is_linear(x)
    cur_is_linear = work_layout["$(x)axis"]["type"] == "linear"
    if !cur_is_linear
      @assert work_layout["$(x)axis"]["type"] == "log"
    end
    delete!(work_layout["$(x)axis"], "type")

    cur_is_linear
  end

  is_linear_x = _is_linear("x")
  is_linear_y = _is_linear("y")

  if is_linear_x
    origin_x = min_x
    width = max_x - min_x
  else
    origin_x = log10(min_x)
    width = log10(max_x) - log10(min_x)
  end

  if is_linear_y
    origin_y = min_y
    height = max_y - min_y
  else
    origin_y = log10(min_y)
    height = log10(max_y) - log10(min_y)
  end

  # ignore desired format for repl
  delete!(work_layout["xaxis"], "exponentformat")
  delete!(work_layout["yaxis"], "exponentformat")

  if haskey(work_layout["xaxis"], "title")
    cur_xlabel = work_layout["xaxis"]["title"]
    delete!(work_layout["xaxis"], "title")
  else
    cur_xlabel = ""
  end

  if haskey(work_layout["yaxis"], "title")
    cur_ylabel = work_layout["yaxis"]["title"]
    delete!(work_layout["yaxis"], "title")
  else
    cur_ylabel = ""
  end

  # unicode plots has no concept of ticks
  delete!(work_layout["xaxis"], "tickvals")
  delete!(work_layout["yaxis"], "tickvals")

  @assert isempty(work_layout["xaxis"])
  @assert isempty(work_layout["yaxis"])

  delete!(work_layout, "xaxis")
  delete!(work_layout, "yaxis")

  cur_canvas = BrailleCanvas(
    size_x, size_y, width=width, height=height,
    origin_x=origin_x, origin_y=origin_y
  )

  if is_linear_x
    if is_linear_y
      border = :solid
    else
      border = :ascii
    end
  else
    if is_linear_y
      border = :dashed
    else
      border = :dotted
    end
  end

  if !is_linear_x
    if strip(cur_xlabel) == ""
      cur_xlabel = "log"
    else
      cur_xlabel = "log " * cur_xlabel
    end
  end

  if haskey(work_layout, "title")
    cur_title = work_layout["title"]
    delete!(work_layout, "title")
  else
    cur_title = ""
  end

  cur_plot = UnicodePlots.Plot(
    cur_canvas; title = cur_title,
    xlabel = cur_xlabel, border=border
  )

  function _annotate!(varargs...; color=:light_black)
    UnicodePlots.annotate!(cur_plot, varargs..., color=color)
  end

  ylabel_index = Int(floor(nrows(cur_canvas)/2))

  if strip(cur_ylabel) == ""
    if !is_linear_y
      _annotate!(:l, ylabel_index, "log", color=:gray)
    end
  else
    _annotate!(:l, ylabel_index, cur_ylabel, color=:gray)

    if !is_linear_y
      log_label = "log"
      cur_diff = length(cur_ylabel) - length(log_label)

      log_label = lpad(log_label, length(log_label) + Int(floor(cur_diff/2)))
      log_label = rpad(log_label, length(cur_ylabel))

      _annotate!(:l, ylabel_index-1, log_label, color=:gray)
    end
  end

  function _clean_string(cur_value)
    rounded_value = round(cur_value, sigdigits=3)

    if ( abs(rounded_value) % 1 ) < ( 1e3 * eps() )
      string_value = string(Int(rounded_value))
    elseif 0.1 <= abs(rounded_value) < 1e4
      string_value = string(rounded_value)
    else
      tmp_coeff, tmp_exponent = split(@sprintf("%.2e", rounded_value), "e")
      tmp_coeff = rstrip(rstrip(tmp_coeff, '0'), '.')

      string_value = tmp_coeff * "e" * tmp_exponent
    end

    string_value
  end

  _annotate!(:l, nrows(cur_canvas), _clean_string(min_y))
  _annotate!(:l, 1, _clean_string(max_y))
  _annotate!(:bl, _clean_string(min_x))
  _annotate!(:br, _clean_string(max_x))

  return (cur_plot, has_legend, is_linear_x, is_linear_y, max_label_length)

end


text:
#=
This file tests the save and load feature of the Mapper.
=#
@testset "Testing Save/Load" begin
    # Create two maps that are the same
    m = place!(Example1.make_map(); move_attempts = 5000)
    sm = Mapper2.SAStruct(m)

    # Place struct sm
    Mapper2.SA.place!(sm; move_attempts = 5000)
    # Record the placed sm into map m
    Mapper2.SA.record(m, sm)
    route!(m)
    # Save the placement
    save(m, "tests")

    ##################################
    # Now create a new testmap.
    ##################################
    n = Example1.make_map()
    # Load the placement into "n"
    load(n, "tests")
    # Create a SA struct for the new map.
    sn = Mapper2.SAStruct(n)
    # Preplace "n" into "sn"
    Mapper2.SA.preplace(n, sn)
    # Assert the two costs are the same.
    @test Mapper2.SA.map_cost(sn) == Mapper2.SA.map_cost(sm)
    # Get the link histograms for the two maps. Make sure they match
    lm = MapperCore.global_link_histogram(m)
    ln = MapperCore.global_link_histogram(n)
    @test lm == ln

    rm("tests.jls")
end