text:
### Deprecated while v0.7 was latest

function ∇softmax(Δ, x; dims = 1)
    # This 2-arg version recomputes the forward pass, which is slow.
    # Removed from use in 0.7, but only prints a warning during 0.8:
    Base.depwarn("`∇softmax(Δ, x)` without `y = softmax(x)` argument is deprecated, as this is inefficient, please use `∇softmax_data(dy, y)`", :∇softmax)
    ∇softmax(Δ, x, softmax(x; dims); dims)
end
∇softmax!(Δ, x; dims = 1) = Δ .= ∇softmax(Δ, x; dims)
∇softmax!(out, Δ, x; dims = 1) = out .= ∇softmax(Δ, x; dims)

function ∇logsoftmax(Δ, x; dims = 1)
    Base.depwarn("`∇logsoftmax(Δ, x)` without `y = logsoftmax(x)` argument is deprecated, please use `∇logsoftmax_data(dy, y)`", :∇logsoftmax)
    ∇logsoftmax(Δ, x, logsoftmax(x; dims); dims)
end
∇logsoftmax!(Δ, x; dims = 1) = Δ .= ∇logsoftmax(Δ, x; dims)
∇logsoftmax!(out, Δ, x; dims = 1) = out .= ∇logsoftmax(Δ, x; dims)

### Deprecated while v0.8 was latest

export ∇softmax,
    ∇softmax!,
    logsoftmax,
    logsoftmax!,
    ∇logsoftmax,
    ∇logsoftmax!

function ∇softmax!(out::AbstractArray, Δ::AbstractArray,
                    x::AbstractArray, y::AbstractArray; dims = 1)
    Base.depwarn("`∇softmax!(dx, dy, x, y)` is deprecated, just use `∇softmax_data(dy, y)`", :∇softmax!)
    # Removed because using a mutating function blocks 2nd derivatives, and
    # the CUDA overload was slow anyway, https://github.com/FluxML/NNlibCUDA.jl/issues/30
    out .= Δ .* y
    out .= out .- y .* sum(out; dims)
end

function ∇logsoftmax!(out::AbstractArray, Δ::AbstractArray,
                    x::AbstractArray, y::AbstractArray; dims = 1)
    Base.depwarn("`∇logsoftmax!(dx, dy, x, y)` is deprecated, just use `∇logsoftmax_data(dy, y)`", :∇softmax!)
    out .= Δ .- sum(Δ; dims) .* exp.(y)
end

function ∇softmax(dy::AbstractArray{T}, x::AbstractArray, y::AbstractArray{S}; dims = 1) where {T,S}
    # Removed because there's no need to close over `x` here, that was done only to distinguish
    # this from `∇softmax(Δ, x; dims = 1)` which re-computed `y = softmax(x)`, which is slow.
    Base.depwarn("`∇softmax(dy, x, y)` should be replaced with `∇softmax_data(dy, y)`", :∇softmax)
    ∇softmax_data(dy, y)
end

function ∇logsoftmax(dy::AbstractArray, x::AbstractArray, y::AbstractArray; dims = 1)
    Base.depwarn("`∇logsoftmax(dy, x, y)` should be replaced with `∇logsoftmax_data(dy, y)`", :∇softmax)
    ∇logsoftmax_data(dy, y)
end


text:
export filter

import DataFrames: filter

function filter(df::AbstractDataFrame, arg; kwargs...)
    filter(arg, df; kwargs...)
end

if false
    using Pkg
    Pkg.activate("c:/git/DataConvenience")
    using DataFrames, DataConvenience
    df = DataFrame(a=1:3)

    filter(df, :a => ==(1))

    @> df begin
        filter(:a => ==(1))
    end
end

ext: jl

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
load("debug.jl")

module TestDebugEval
import Base.*
import Debug.*

macro assign(lhs, rhs)
    esc(quote
        ($lhs) = ($rhs)   # f_line - 5
    end)
end

const f_line = 15 # line number of the next line:
@debug function f(n) # this must be line number f_line
    x = 0     # line number f_line + 1
    for k=1:n                      # 2
        x += k                     # 3
    end                            # 4
    x = x*x                        # 5
    let y=2x                       # 6
        y                          # 7
    end                            # 8
    try                            # 9
        error("throw something")   #10
    catch e                        #11
        x                          #12
    end                            #13
    [begin x^2 end for x=1:3]      #14
    g = (x, y)->begin              #15
        x = -1                     #16
        a, b = 1, 2                #17
        x                          #18
    end                            #19
    g(n, n)                        #20
    @assign w x/2                  #21
    x                              #22
end

function debug_hook(line::Int, file, scope::Scope)
    line -= f_line
    print(line, ":")

    if (line ==  2) debug_eval(scope, :(x += 1)) end

    if (line !=  1) debug_eval(scope, :(print("\tx = ", x))) end
    if (line ==  3) debug_eval(scope, :(print("\tk = ", k))) end
    if (line ==  7) debug_eval(scope, :(print("\ty = ", y))) end
    if (line == 12) debug_eval(scope, :(print("\te = ", e))) end
    if (line >  15) debug_eval(scope, :(print("\tg = ", g))) end
    if (16 <= line <= 18) debug_eval(scope, :(print("\ty = ", y))) end
    if (line == 18) debug_eval(scope, :(print("\ta = ", a, "\tb = ", b, ))) end
#    if (line >  21) debug_eval(scope, :(print("\tw = ", w))) end
    println()
end

f(3)

end


text: nthroots(n::Integer) = [ cospi(2k/n)+sinpi(2k/n)im for k = 0:n-1 ]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@doc raw"""
    function layerUpdateParams!(
        model::Model,
        cLayer::FoB,
        cnt::Integer = -1;
        tMiniBatch::Integer = 1,
        kwargs...,
    ) where {FoB <: Union{FCLayer, BatchNorm}}

update trainable parameters for `FCLayer` and `BatchNorm` layers
"""
function layerUpdateParams!(
    model::Model,
    cLayer::FoB,
    cnt::Integer = -1;
    tMiniBatch::Integer = 1,
    kwargs...,
) where {FoB <: Union{FCLayer, BatchNorm}}

    optimizer = model.optimizer
    α = model.α
    β1, β2, ϵAdam = model.β1, model.β2, model.ϵAdam

    if cLayer.updateCount >= cnt
        return nothing
    end #if cLayer.updateCount >= cnt

    if !all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)
        return nothing
    end

    #initialize the needed variables to hold the corrected values
    #it is being init here cause these are not needed elsewhere
    VCorrected = Dict(:dw=>similar(cLayer.dW), :db=>similar(cLayer.dB))
    SCorrected = Dict(:dw=>similar(cLayer.dW), :db=>similar(cLayer.dB))
    if optimizer==:adam || optimizer==:momentum

        if tMiniBatch == 1
            cLayer.V[:dw] .= 0
            cLayer.V[:db] .= 0
        end
        cLayer.V[:dw] .= β1 .* cLayer.V[:dw] .+ (1-β1) .* cLayer.dW
        cLayer.V[:db] .= β1 .* cLayer.V[:db] .+ (1-β1) .* cLayer.dB

        ##correcting
        VCorrected[:dw] .= cLayer.V[:dw] ./ (1-β1^tMiniBatch)
        VCorrected[:db] .= cLayer.V[:db] ./ (1-β1^tMiniBatch)

        if optimizer==:adam

            if tMiniBatch == 1
                cLayer.S[:dw] .= 0
                cLayer.S[:db] .= 0
            end
            cLayer.S[:dw] .= β2 .* cLayer.S[:dw] .+ (1-β2) .* (cLayer.dW.^2)
            cLayer.S[:db] .= β2 .* cLayer.S[:db] .+ (1-β2) .* (cLayer.dB.^2)

            ##correcting
            SCorrected[:dw] .= cLayer.S[:dw] ./ (1-β2^tMiniBatch)
            SCorrected[:db] .= cLayer.S[:db] ./ (1-β2^tMiniBatch)

            ##update parameters with adam
            cLayer.W .-= (α .* (VCorrected[:dw] ./ (sqrt.(SCorrected[:dw]) .+ ϵAdam)))
            cLayer.B .-= (α .* (VCorrected[:db] ./ (sqrt.(SCorrected[:db]) .+ ϵAdam)))

        else#if optimizer==:momentum

            cLayer.W .-= (α .* VCorrected[:dw])
            cLayer.B .-= (α .* VCorrected[:db])

        end #if optimizer==:adam
    else
        cLayer.W .-= (α .* cLayer.dW)
        cLayer.B .-= (α .* cLayer.dB)
    end #if optimizer==:adam || optimizer==:momentum

    cLayer.updateCount += 1

    return nothing
end #updateParams!

function layerUpdateParams!(
                       model::Model,
                       cLayer::IoA,
                       cnt::Integer = -1;
                       tMiniBatch::Integer = 1,
                       kwargs...,
                       ) where {IoA <: Union{Input, Activation, AddLayer, Flatten, ConcatLayer}}

    # optimizer = model.optimizer
    # α = model.α
    # β1, β2, ϵAdam = model.β1, model.β2, model.ϵAdam

    if cLayer.updateCount >= cnt
        return nothing
    end #if cLayer.updateCount >= cnt

    if !all(i->(i.updateCount==cLayer.nextLayers[1].updateCount), cLayer.nextLayers)
        return nothing
    end

    cLayer.updateCount += 1

    return nothing
end #updateParams!

export layerUpdateParams!


text:
#=
Part of this work is derived from astropy/photutils. The relevant derivations
are considered under a BSD 3-clause license. =#

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

"""
    CircularAperture(x, y, r)
    CircularAperture(position, r)

A circular aperture.

A circular aperture with radius `r`. `r` must be greater than or equal to 0.

# Examples
```jldoctest
julia> ap = CircularAperture(0, 0, 10)
21×21 CircularAperture{Int64} with indices -10:10×-10:10:
 0.0        0.0       0.0         …  0.0         0.0       0.0
 0.0        0.0       0.0            0.0         0.0       0.0
 0.0        0.0       0.0            0.0         0.0       0.0
 0.0        0.0       0.00571026     0.00571026  0.0       0.0
 0.0        0.0       0.491844       0.491844    0.0       0.0
 0.0        0.170878  0.982952    …  0.982952    0.170878  0.0
 0.0        0.659735  1.0            1.0         0.659735  0.0
 0.0590655  0.975524  1.0            1.0         0.975524  0.0590655
 0.293527   1.0       1.0            1.0         1.0       0.293527
 0.445643   1.0       1.0            1.0         1.0       0.445643
 ⋮                                ⋱                        ⋮
 0.293527   1.0       1.0            1.0         1.0       0.293527
 0.0590655  0.975524  1.0            1.0         0.975524  0.0590655
 0.0        0.659735  1.0            1.0         0.659735  0.0
 0.0        0.170878  0.982952    …  0.982952    0.170878  0.0
 0.0        0.0       0.491844       0.491844    0.0       0.0
 0.0        0.0       0.00571026     0.00571026  0.0       0.0
 0.0        0.0       0.0            0.0         0.0       0.0
 0.0        0.0       0.0            0.0         0.0       0.0
 0.0        0.0       0.0         …  0.0         0.0       0.0
```
"""
struct CircularAperture{T<:Number} <: AbstractAperture{T}
    x::T
    y::T
    r::T
end

function overlap(ap::CircularAperture, i, j)
    dist = sqrt((i - ap.y)^2 + (j - ap.x)^2)
    dr = sqrt(2) / 2 # corner-center distance of pixel
    dist > ap.r + dr && return Outside
    dist < ap.r - dr && return Inside
    return Partial
end

partial(ap::CircularAperture, x, y) = circular_overlap_single_exact(x - 0.5, y - 0.5, x + 0.5, y + 0.5, ap.r)
partial(sub_ap::Subpixel{T,<:CircularAperture}, x, y) where {T} = circular_overlap_single_subpixel(x - 0.5, y - 0.5, x + 0.5, y + 0.5, sub_ap.ap.r, sub_ap.N)

CircularAperture(center, r) = CircularAperture(center..., r)
CircularAperture(x, y, r) = CircularAperture(promote(x, y, r)...)

function Base.show(io::IO, c::CircularAperture)
    print(io, "CircularAperture($(c.x), $(c.y), r=$(c.r))")
end

function bounds(c::CircularAperture)
    _xmin = c.x - c.r - 0.5
    _ymin = c.y - c.r - 0.5
    xmin = isinteger(_xmin) ? ceil(Int, _xmin) + 1 : ceil(Int, _xmin)
    ymin = isinteger(_ymin) ? ceil(Int, _ymin) + 1 : ceil(Int, _ymin)
    xmax = ceil(Int, c.x + c.r - 0.5)
    ymax = ceil(Int, c.y + c.r - 0.5)
    return (xmin, xmax, ymin, ymax)
end

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

"""
    CircularAnnulus(x, y, r_in, r_out)
    CircularAnnulus(position, r_in, r_out)

A circular annulus with inner radius `r_in` and outer radius `r_out`. 0 ≤ `r_in` ≤ `r_out`.

# Examples
```jldoctest
julia> ap = CircularAnnulus(0, 0, 5, 10)
21×21 CircularAnnulus{Int64} with indices -10:10×-10:10:
 0.0        0.0       0.0         …  0.0         0.0       0.0
 0.0        0.0       0.0            0.0         0.0       0.0
 0.0        0.0       0.0            0.0         0.0       0.0
 0.0        0.0       0.00571026     0.00571026  0.0       0.0
 0.0        0.0       0.491844       0.491844    0.0       0.0
 0.0        0.170878  0.982952    …  0.982952    0.170878  0.0
 0.0        0.659735  1.0            1.0         0.659735  0.0
 0.0590655  0.975524  1.0            1.0         0.975524  0.0590655
 0.293527   1.0       1.0            1.0         1.0       0.293527
 0.445643   1.0       1.0            1.0         1.0       0.445643
 ⋮                                ⋱                        ⋮
 0.293527   1.0       1.0            1.0         1.0       0.293527
 0.0590655  0.975524  1.0            1.0         0.975524  0.0590655
 0.0        0.659735  1.0            1.0         0.659735  0.0
 0.0        0.170878  0.982952    …  0.982952    0.170878  0.0
 0.0        0.0       0.491844       0.491844    0.0       0.0
 0.0        0.0       0.00571026     0.00571026  0.0       0.0
 0.0        0.0       0.0            0.0         0.0       0.0
 0.0        0.0       0.0            0.0         0.0       0.0
 0.0        0.0       0.0         …  0.0         0.0       0.0
```
"""
struct CircularAnnulus{T <: Number} <: AbstractAperture{T}
    x::T
    y::T
    r_in::T
    r_out::T
end

CircularAnnulus(center, r_in, r_out) = CircularAnnulus(center..., r_in, r_out)
CircularAnnulus(x, y, r_in, r_out) = CircularAnnulus(promote(x, y, r_in, r_out)...)

function Base.show(io::IO, c::CircularAnnulus)
    print(io, "CircularAnnulus($(c.x), $(c.y), r_in=$(c.r_in), r_out=$(c.r_out))")
end

function overlap(ap::CircularAnnulus, i, j)
    dist = sqrt((i - ap.y)^2 + (j - ap.x)^2)
    dr = sqrt(2) / 2 # corner-center distance of pixel
    ap.r_in - dr < dist < ap.r_out + dr || return Outside
    ap.r_in + dr < dist < ap.r_out - dr && return Inside
    return Partial
end

function partial(ap::CircularAnnulus, x, y)
    f1 = circular_overlap_single_exact(x - 0.5, y - 0.5, x + 0.5, y + 0.5, ap.r_out)
    f2 = circular_overlap_single_exact(x - 0.5, y - 0.5, x + 0.5, y + 0.5, ap.r_in)
    return f1 - f2
end

function partial(sub_ap::Subpixel{T,<:CircularAnnulus}, x, y) where T
    ap = sub_ap.ap
    f1 = circular_overlap_single_subpixel(x - 0.5, y - 0.5, x + 0.5, y + 0.5, ap.r_out, sub_ap.N)
    f2 = circular_overlap_single_subpixel(x - 0.5, y - 0.5, x + 0.5, y + 0.5, ap.r_in, sub_ap.N)
    return f1 - f2
end

function bounds(c::CircularAnnulus)
    xmin = ceil(Int, c.x - c.r_out - 0.5)
    ymin = ceil(Int, c.y - c.r_out - 0.5)
    xmax = ceil(Int, c.x + c.r_out - 0.5)
    ymax = ceil(Int, c.y + c.r_out - 0.5)
    return (xmin, xmax, ymin, ymax)
end


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
@cuda_include "cudaReduce.cu"

### Constructors

for f in [:sum, :minimum, :maximum]
  qf = QuoteNode(f)
  qf_i = QuoteNode(Symbol(f, "_i"))
  _priname = Symbol("_", f, "!")
  @eval begin
    Base.$f(x::ANode) = ANode($qf, (x,), (1,))

    function Base.$f(x::ANode, i::Integer)
      nx = size(x, i)
      ny = prod(Int[size(x,j) for j in 1:(i-1)])
      nz = prod(Int[size(x,j) for j in (i+1):ndims(x)])
      incx = ny
      incy = i == 1 ? size(x, 1): 1
      incz = ny*nx

      sz = [size(x)...]
      sz[i] = 1

      ANode($qf_i, (x,), (sz...), (nx, ny, nz, incx, incy, incz))
    end

    do_forward!(n::ANode{$qf}, out, in) = $_priname(length(in), out, in)
    do_forward!(n::ANode{$qf_i}, out, in) = $_priname(
        n.arg[1], n.arg[2], n.arg[3], out, in, n.arg[4], n.arg[5], n.arg[6])
  end
end

### sum

function _sum!{T}(n::Integer, y::Array{T}, x::Array{T}, a::Real=zero(T))
  @assert n == length(x) && 1 == length(y)
  y[1] = a == 0 ? 0 : y[1]*a
  for i in 1:n @inbounds y[1] += x[i] end
  nothing
end

function _sum!{T}(nx::Integer, ny::Integer, nz::Integer, y::Array{T}, x::Array{T}, incx::Integer, incy::Integer, incz::Integer, a::Real=zero(T))
  @assert ny*nz == length(y)

  for i in 1:nz
    for j in 1:ny
      @inbounds yi = j + ny*(i-1)
      xi = (i-1)*incz + (j-1)*incy + 1

      s = a == 0 ? zero(eltype(y)) : convert(eltype(y), y[yi]*a)
      for k in 1:nx
        s += x[xi]
        xi += incx
      end

      @inbounds y[yi] = s
    end
  end
end

@cuda begin
_sum!{T}(n::Integer, y::CudaArray{T}, x::CudaArray{T}, a::Real=zero(T)) = _sum!(n, 1, 1, y, x, 1, 1, 1, a)

const cuda_sum_32_float = CuFunction(md, "cuda_sum_32_float")
_sum!(nx::Integer, ny::Integer, nz::Integer, y::CudaArray{Float32}, x::CudaArray{Float32}, incx::Integer, incy::Integer, incz::Integer, a::Real=zero(Float32)) =
    cudacall(cuda_sum_32_float, (1,ny,nz), 32, (Csize_t, Ptr{Cfloat}, Ptr{Cfloat}, Cfloat, UInt64, UInt64, UInt64), nx, x, y, a, incx, incy, incz; shmem_bytes=32*4, stream=STREAM[1])

const cuda_sum_32_double = CuFunction(md, "cuda_sum_32_double")
_sum!(nx::Integer, ny::Integer, nz::Integer, y::CudaArray{Float64}, x::CudaArray{Float64}, incx::Integer, incy::Integer, incz::Integer, a::Real=zero(Float64)) =
    cudacall(cuda_sum_32_double, (1,ny,nz), 32, (Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, UInt64, UInt64, UInt64), nx, x, y, a, incx, incy, incz; shmem_bytes=32*8, stream=STREAM[1])
end

function add_constant{T}(n::Csize_t, a::Array{T}, b::Array{T})
    @assert n == length(a) && 1 == length(b)
    a[:] += b[1]
end

@cuda_gsl add_constant{T<:Real}(n::Csize_t, a::&T, b::&T) "a[i] += b[0]"

gradient_node(n::ANode{:sum}, wrt::ANode, b::ANode) = broadcast_node(b, size(wrt))
gradient_node(n::ANode{:sum_i}, wrt::ANode, b::ANode) = broadcast_node(b, size(wrt))

### minimum

function _minimum!{T}(n::Integer, y::Array{T}, x::Array{T})
  @assert n == length(x) && 1 == length(y)
  y[1] = x[1]
  for i in 2:n
    @inbounds if x[i] < y[1] y[1] = x[i] end
  end
  nothing
end

function _minimum!{T}(nx::Integer, ny::Integer, nz::Integer, y::Array{T}, x::Array{T}, incx::Integer, incy::Integer, incz::Integer)
  @assert ny*nz == length(y)

  for i in 1:nz
    for j in 1:ny
      @inbounds yi = j + ny*(i-1)
      xi = (i-1)*incz + (j-1)*incy + 1

      s = x[xi]
      xi += incx
      for k in 2:nx
        if x[xi] < s s = x[xi] end
        xi += incx
      end

      @inbounds y[yi] = s
    end
  end
end

@cuda begin
_minimum!{T}(n::Integer, y::CudaArray{T}, x::CudaArray{T}) = _minimum!(n, 1, 1, y, x, 1, 1, 1)

const cuda_minimum_32_float = CuFunction(md, "cuda_minimum_32_float")
_minimum!(nx::Integer, ny::Integer, nz::Integer, y::CudaArray{Float32}, x::CudaArray{Float32}, incx::Integer, incy::Integer, incz::Integer) =
    cudacall(cuda_minimum_32_float, (1,ny,nz), 32, (Csize_t, Ptr{Cfloat}, Ptr{Cfloat}, UInt64, UInt64, UInt64), nx, x, y, incx, incy, incz; shmem_bytes=32*4, stream=STREAM[1])

const cuda_minimum_32_double = CuFunction(md, "cuda_minimum_32_double")
_minimum!(nx::Integer, ny::Integer, nz::Integer, y::CudaArray{Float64}, x::CudaArray{Float64}, incx::Integer, incy::Integer, incz::Integer) =
    cudacall(cuda_minimum_32_double, (1,ny,nz), 32, (Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, UInt64, UInt64, UInt64), nx, x, y, incx, incy, incz; shmem_bytes=32*8, stream=STREAM[1])
end

gradient_node(n::ANode{:minimum}, wrt::ANode, b::ANode) = b.*(n .== n.input[1])
gradient_node(n::ANode{:minimum_i}, wrt::ANode, b::ANode) = b.*(n .== n.input[1])

### maximum

function _maximum!{T}(n::Integer, y::Array{T}, x::Array{T})
  @assert n == length(x) && 1 == length(y)
  y[1] = x[1]
  for i in 2:n
    @inbounds if x[i] > y[1] y[1] = x[i] end
  end
  nothing
end

function _maximum!{T}(nx::Integer, ny::Integer, nz::Integer, y::Array{T}, x::Array{T}, incx::Integer, incy::Integer, incz::Integer)
  @assert ny*nz == length(y)

  for i in 1:nz
    for j in 1:ny
      @inbounds yi = j + ny*(i-1)
      xi = (i-1)*incz + (j-1)*incy + 1

      s = x[xi]
      xi += incx
      for k in 2:nx
        if x[xi] > s s = x[xi] end
        xi += incx
      end

      @inbounds y[yi] = s
    end
  end
end

@cuda begin
_maximum!{T}(n::Integer, y::CudaArray{T}, x::CudaArray{T}) = _maximum!(n, 1, 1, y, x, 1, 1, 1)

const cuda_maximum_32_float = CuFunction(md, "cuda_maximum_32_float")
_maximum!(nx::Integer, ny::Integer, nz::Integer, y::CudaArray{Float32}, x::CudaArray{Float32}, incx::Integer, incy::Integer, incz::Integer) =
    cudacall(cuda_maximum_32_float, (1,ny,nz), 32, (Csize_t, Ptr{Cfloat}, Ptr{Cfloat}, UInt64, UInt64, UInt64), nx, x, y, incx, incy, incz; shmem_bytes=32*4, stream=STREAM[1])

const cuda_maximum_32_double = CuFunction(md, "cuda_maximum_32_double")
_maximum!(nx::Integer, ny::Integer, nz::Integer, y::CudaArray{Float64}, x::CudaArray{Float64}, incx::Integer, incy::Integer, incz::Integer) =
    cudacall(cuda_maximum_32_double, (1,ny,nz), 32, (Csize_t, Ptr{Cdouble}, Ptr{Cdouble}, UInt64, UInt64, UInt64), nx, x, y, incx, incy, incz; shmem_bytes=32*8, stream=STREAM[1])
end

gradient_node(n::ANode{:maximum}, wrt::ANode, b::ANode) = b.*(n .== n.input[1])
gradient_node(n::ANode{:maximum_i}, wrt::ANode, b::ANode) = b.*(n .== n.input[1])


text:
# Characters
const SUITCHAR=collect("CDHSN")
const CARDCHAR=collect("23456789TJQKA")
const HEXCHAR=collect("0123456789ABCDEF")
const PLAYERCHAR=collect("WNES")
hex2int(c::Char)=(c <= '9' ? c - '0' : c - '7')
int2hex(i::Integer)=HEXCHAR[i+1]

# Trump suits
const CLUBS=1
const DIAMONDS=2
const HEARTS=3
const SPADES=4
const NOTRUMP=5

# Regular bids are represented by integers 1:35
const TRUMPBIDS=35
makebid(level,suit)=(5*(level-1)+suit)
bidlevel(bid)=1+div(bid-1,5)
bidtrump(bid)=mod1(bid,5)
bidtrumpchar(bid)=SUITCHAR[bidtrump(bid)]
bidlevelchar(bid)='0'+bidlevel(bid)

# Cards
makecard(suit,value)=(value + 13*(suit-1))
cardvalue(card)=mod1(card,13)
cardsuit(card)=1+div(card-1,13)
cardvaluechar(card)=CARDCHAR[cardvalue(card)]
cardsuitchar(card)=SUITCHAR[cardsuit(card)]

# Three extra bids
const PASS=36
const DOUBLE=37
const REDOUBLE=38
const NUMBIDS=38

# Players
const WEST=1
const NORTH=2
const EAST=3
const SOUTH=4

# Double dummy data
# using StaticArrays
# struct DoubleDummy2
#     hands::SVector{4,UInt64} # Using 52 bits of UInt64 to encode a hand
#     results::SMatrix{4,5,UInt8}  # (4,5) array of NS tricks
# end

max_stars_repo_name: JuliaTagBot/play

max_issues_repo_name: JuliaTagBot/play

max_forks_repo_name: JuliaTagBot/play


text:
####################################################
"""
    plt = plotstrack(track::VisClaw.Track, index::AbstractVector=1:length(track.lon); kwargs...)

    plotstrack!(plt::Plots.Plot, track::VisClaw.Track, index::AbstractVector=1:length(track.lon); kwargs...)

plot a typhoon/hurricane track using Plots
"""
function plotstrack!(plt, track::VisClaw.Track, index::AbstractVector=1:length(track.lon); kwargs...)
    # parse keyword args
    kwdict = KWARG(kwargs)
	label, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_label, "")

	# plot
	plt = Plots.plot!(plt, track.lon[index], track.lat[index]; axis_ratio=:equal, label=label, kwdict...)
	return plt
end
####################################################
"""
$(@doc plotstrack!)
"""
plotstrack(track::VisClaw.Track, index::AbstractVector=1:length(track.lon); kwargs...) =
plotstrack!(Plots.plot(), track, index; kwargs...)
####################################################


text:
export
    Importance,
    Rejection,
    LW,
    rejection_proposal,
    lw_proposal,
    make_custom_proposal

"""
    Importance <: InstantAlgorithm

Representation of an importance sampling algorithm.

# arguments
- proposal_function Specifies how the algorithm should make proposals. This is a function
that takes a runtime and an instance and returns a proposer.
The proposer takes parent values and proposes a value for the instance along with a log
score.
- num_particles The number of completed particles to use. This is not necessarily the
number attempted. If there are rejections, the algorithm will continue to create particles
until `num_particles` have been completed. Warning: With impossible evidence, the process
will not terminate.
"""
struct Importance <: InstantAlgorithm
    proposal_function :: Function
    num_particles :: Int
end

function answer(pf::ProbFunction, ::Importance, runtime::Runtime, instances::Vector{VariableInstance})
    particles = get_state(runtime, :particles)
    score = 0.0
    total = 0.0
    for (s, lw) in zip(particles.samples, particles.log_weights)
        x = [s[get_name(inst)] for inst in instances]
        w = exp(lw)
        if pf.fn(x)
            score += w
        end
        total += w
    end
    return score / total
end

function answer(::Marginal, ::Importance, runtime::Runtime, instance::VariableInstance)
    O = output_type(get_node(instance))
    dict = Dict{O, Float64}()
    particles = get_state(runtime, :particles)
    for (s, lw) in zip(particles.samples, particles.log_weights)
        x = s[get_name(instance)]
        w = exp(lw)
        dict[x] = get(dict, x, 0.0) + w
    end
    z = sum(values(dict))
    vs = O[]
    ps = Float64[]
    for (v, p) in dict
        push!(vs, v)
        push!(ps, p/z)
    end
    return Cat(vs, ps)
end

const Reject = ErrorException("Reject")

"""
    rejection_proposal(::Runtime, instance::VariableInstance, parent_values::Tuple)

Return a proposer and scorer to implement standard rejection sampling from the prior.
It proposes a value for the `instance` from its sfunc, and scores it by the evidence,
if any. If the score is -Infinity, it throws a Reject exception.
"""
function rejection_proposal(runtime::Runtime, instance::VariableInstance)
    proposer(parent_values) = (sample(get_sfunc(instance), parent_values), 0.0)
end

function _get_hard_evidence(runtime, instance)::Union{HardScore, Nothing}
    if has_evidence(runtime, instance)
        ev = get_evidence(runtime, instance)
        if ev isa HardScore
            return ev
        end
    end
    return nothing
end

"""
    lw_proposal(runtime::Runtime, instance::VariableInstance, parent_values::Tuple)

Return a proposer and scorer to implement likelihood weighting.

This proposal scheme is the same as the prior proposal unless a variable has hard evidence.
In the case of hard evidence, the proposer sets the value of the variable to the evidence
value and scores it by the log conditional probability of the evidence given the parent
values.
"""
function lw_proposal(runtime::Runtime, instance::VariableInstance)
    evidence = _get_hard_evidence(runtime, instance)
    if !isnothing(evidence)
        return _hard_evidence_proposal(evidence, instance)
    else
        return rejection_proposal(runtime, instance)
    end
end

function _hard_evidence_proposal(evidence, instance)
    sf = get_sfunc(instance)
    x = evidence.value
    proposer(parent_values) = (x, logcpdf(sf, parent_values, x))
    return proposer
end

"""
    make_custom_proposal(custom_sfs::Dict{Symbol, SFunc})

Create a proposal function for a custom proposal scheme.

Returns a proposal function that can be provided to the Importance constructor.
Evidence is handled similar to `lw`, except that the custom proposal is used for soft
evidence.

# Arguments
- custom_sfs A dictionary mapping variable names to custom sfuncs used for their proposal.
Need not be complete; if a variable is not in this dictionary, its standard sfunc will be
used.
"""
function make_custom_proposal(custom_sfs::Dict{Symbol, SFunc})
    function proposal(runtime, instance)
        evidence = _get_hard_evidence(runtime, instance)
        if !isnothing(evidence)
            return _hard_evidence_proposal(evidence, instance)
        else
            name = get_name(instance)
            if name in keys(custom_sfs)
                prior_sf = get_sfunc(instance)
                proposal_sf = custom_sfs[name]
                function proposer(parent_values)
                    x = sample(proposal_sf, parent_values)
                    l = logcpdf(prior_sf, parent_values, x) -
                            logcpdf(proposal_sf, parent_values, x)
                    return (x, l)
                end
                return proposer
            else
                return rejection_proposal(runtime, instance)
            end
        end
    end
    return proposal
end

Rejection(num_particles) = Importance(rejection_proposal, num_particles)

LW(num_particles) = Importance(lw_proposal, num_particles)

function _importance(runtime::Runtime, num_samples::Int, proposal_function::Function,
    samples::Vector{Dict{Symbol, Any}}, lws)
    net = runtime.network
    nodes = topsort(get_initial_graph(net))
    proposers = Function[]
    evidences = Score[]
    interventions = Union{Dist,Nothing}[]
    for v in nodes
        if v isa Variable
            inst = current_instance(runtime, v)
            push!(proposers, proposal_function(runtime, inst))
            if has_evidence(runtime, inst)
                push!(evidences, get_evidence(runtime, inst))
            else
                push!(evidences, FunctionalScore{output_type(v)}(x -> 1.0))
            end
            if has_intervention(runtime, inst)
                push!(interventions, get_intervention(runtime, inst))
            else
                push!(interventions, nothing)
            end
        end
    end
    s = 1
    while s <= num_samples
        try
            vnum = 1
            for v in nodes
                if v isa Variable
                    inst = current_instance(runtime, v)
                    if !isnothing(interventions[vnum])
                        iv = interventions[vnum]
                        samples[s][v.name] = sample(iv, ())
                    else
                        proposer = proposers[vnum]
                        pars = get_initial_parents(net, v)
                        parvals = tuple([samples[s][p.name] for p in pars]...)
                        (x,lw) = proposer(parvals)
                        pe = get_log_score(evidences[vnum], x)
                        if !isfinite(pe)
                            throw(Reject)
                        end
                        samples[s][v.name] = x
                        lws[s] += lw + pe
                    end
                    vnum += 1
                end
            end
            s += 1
        catch e
            if e != Reject
                throw(e)
            end
        end
    end
    for v in nodes
        if v isa Variable
            ps = Dict{output_type(v), Float64}()
            for i in 1:length(samples)
                x = samples[i][v.name]
                ps[x] = get(ps, x, 0.0) + exp(lws[i])
            end
            i = current_instance(runtime, v)
            set_value!(runtime, i, :belief, Cat(ps))
        end
    end

    log_prob_evidence = logsumexp(lws) - log(num_samples)
    set_state!(runtime, :log_prob_evidence, log_prob_evidence)
    set_state!(runtime, :particles, Particles(samples, lws))
end

function infer(algorithm::Importance, runtime::InstantRuntime,
    evidence::Dict{Symbol, Score} = Dict{Symbol, Score}(),
    interventions::Dict{Symbol, Dist} = Dict{Symbol, Dist}(),
    placeholder_beliefs::Dict{Symbol, Dist} = get_placeholder_beliefs(runtime))
    net = get_network(runtime)
    nodes = get_nodes(net)
    ensure_all!(runtime)
    # See comment in instantalgorithm.jl
    # During particle filtering, importance needs joint samples of placeholders, not marginals
    # So we first check to see if particles already exists, and only add placeholder beliefs if
    # they're not found in existing particles
    samples = Dict{Symbol, Any}[]
    lws = Float64[]
    if has_state(runtime, :particles)
        particles = get_state(runtime, :particles)
        # Copy over only values of nodes in the current runtime
        for i in 1:algorithm.num_particles
            if length(particles.samples) > 0
                index_into_particles = (i-1) % length(particles.samples) + 1
                newsample = Dict{Symbol, Any}()
                oldsample = particles.samples[index_into_particles]
                for n in nodes
                    if n.name in keys(oldsample)
                        newsample[n.name] = oldsample[n.name]
                    end
                end
                push!(samples, newsample)
                push!(lws, particles.log_weights[index_into_particles])
            else
                push!(samples, Dict{Symbol, Any}())
                push!(lws, 0.0)
            end
        end
    else
        for i in 1:algorithm.num_particles
            push!(samples, Dict{Symbol, Any}())
        end
        lws = zeros(Float64, algorithm.num_particles)
    end

    placeholders = get_placeholders(net)
    for ph in placeholders
        pi = placeholder_beliefs[ph.name]
        for i in 1:algorithm.num_particles
            if !(ph.name in keys(samples[i]))
                samples[i][ph.name] = sample(pi, ())
            end
        end
    end
    for (n,e) in evidence
        v = get_node(net, n)
        inst = current_instance(runtime, v)
        post_evidence!(runtime, inst, e)
    end
    for (n,i) in interventions
        v = get_node(net, n)
        inst = current_instance(runtime, v)
        post_intervention!(runtime, inst, i)
    end
    _importance(runtime, algorithm.num_particles, algorithm.proposal_function, samples, lws)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 18f8736c-8754-11eb-3a0c-c729d6c6b0d8
begin
  using Pkg; Pkg.activate("MLJ_env", shared=true)
  using PlutoUI

    using CategoricalArrays
    # import MLJLinearModels
    import DataFrames
    import CSV
    # import DecisionTree
    using MLJ
    # import MLJLinearModels
    # import MultivariateStats
    # import MLJFlux
    using Plots
end

# ╔═╡ 5aa90fa0-8762-11eb-18a8-a5a45edfd5a7
using Test

# ╔═╡ e3664c92-8753-11eb-325b-95ca5bcab046
md"""
## A workshop introducing the machine learning toolbox MLJ - Part 1 - Data Representation

  - source:[MLJ](https://alan-turing-institute.github.io/MLJ.jl/stable/)
  - based on [Machine Learning in Julia using MLJ, JuliaCon2020](https://github.com/ablaom/MachineLearningInJulia2020)
"""

# ╔═╡ 18e124b4-8754-11eb-07d4-a9eee0b8ba3b
md"""
##### General resources

  - [List of methods introduced in this tutorial](methods.md)
  - [MLJ Cheatsheet](https://alan-turing-institute.github.io/MLJ.jl/dev/mlj_cheatsheet/)
  - [Common MLJ Workflows](https://alan-turing-institute.github.io/MLJ.jl/dev/common_mlj_workflows/)
  - [MLJ manual](https://alan-turing-institute.github.io/MLJ.jl/dev/)
  - [Data Science Tutorials in Julia](https://alan-turing-institute.github.io/DataScienceTutorials.jl/)
"""

# ╔═╡ 992cc480-8785-11eb-0572-c33c1c988efc
html"""

"""

# ╔═╡ 53bc7eb6-8764-11eb-326d-41138236d179
md"""
##### TOC
  - [Data Representation](#data-representation)
  - [Two-dimensional data](#2-dimensional)
  - [Fixing scientific types in tabular data](#fixing-scitype-tabular-data)
  - [Resources or this part](#resources)
  - [Exercises for Part 1](#exercises)
"""

# ╔═╡ 9dbf597a-8764-11eb-01e6-67bac4f668b5
html"""

  back to TOC

"""

# ╔═╡ 18c455fa-8754-11eb-2683-9dfe2784e339
md"""
#### Data Representation

To help you focus on the intended *purpose* or *interpretation* of
data, MLJ models specify data requirements using *scientific types*,
instead of machine types. An example of a scientific type is
 `OrderedFactor`. The other basic "scalar" scientific types are
 illustrated below:
"""

# ╔═╡ 18acb3c8-8754-11eb-2d0a-5d82e2b47551
PlutoUI.LocalResource("assets/scitypes.png")

# ╔═╡ 185c6c10-8754-11eb-16a3-4b2439aa1d84
md"""
 A scientific type is an ordinary Julia type (so it can be used for
 method dispatch, for example) but it usually has no instances. The
 `scitype` function is used to articulate MLJ's convention about how
 different machine types will be interpreted by MLJ models:
"""

# ╔═╡ 7f3e28ec-8754-11eb-3c7e-e71dc9b2ca4e
begin
  time = [2.3, 4.5, 4.2, 1.8, 7.1]
  scitype(time)
end

# ╔═╡ 841c8930-8754-11eb-1596-0535d54e7558
begin
  _height = [185, 153, 163, 114, 180]
  scitype(_height)
end

# ╔═╡ 88b066ce-8754-11eb-1ec8-61fc466493d5
md"""
To fix data which MLJ is interpreting incorrectly, we use the `coerce` method:
"""

# ╔═╡ 8edfcc74-8754-11eb-3ac9-913d6db5d0c7
height = coerce(_height, Continuous)

# ╔═╡ 8ec4a37c-8754-11eb-05c6-e9d5aa55479f
md"""
Here's an example of data we would want interpreted as  `OrderedFactor` but isn't:
"""

# ╔═╡ 8e9259d0-8754-11eb-01ee-cd92c6663059
begin
  _exam_mark = ["rotten", "great", "bla",  missing, "great"]
  scitype(_exam_mark)
end

# ╔═╡ a152aad4-8754-11eb-33d1-cf4df8cd76d2
begin
  exam_mark_ = coerce(_exam_mark, OrderedFactor)
  levels(exam_mark_)
end

# ╔═╡ a134e1fc-8754-11eb-1664-1d96b05f7ca4
md"""
Use `levels!` to put the classes in the right order:
"""

# ╔═╡ ac1420c4-8754-11eb-3751-4df7eff8ada4
begin
  levels!(exam_mark_, ["rotten", "bla", "great"])
  exam_mark_[1] < exam_mark_[2]
end

# ╔═╡ abfabfc6-8754-11eb-0a32-93a7e254063d
levels(exam_mark_)

# ╔═╡ bd7c20be-8754-11eb-07a1-e182ca10e932
md"""
When sub-sampling, no levels are lost:
"""

# ╔═╡ abdefbf4-8754-11eb-1bc5-4950210910d7
levels(exam_mark_[1:2])

# ╔═╡ abc65bbe-8754-11eb-325f-d98618a30eb8
md"""
**Note on binary data.** There is no separate scientific type for binary data. Binary data is `OrderedFactor{2}` or  `Multiclass{2}`. If a binary measure like `truepositive` is a  applied to `OrderedFactor{2}` then the "positive" class is assumed to appear *second* in the ordering. If such a measure is applied to  `Multiclass{2}` data, a warning is issued. $(html"
")
A single `OrderedFactor` can be coerced to a single `Continuous` variable, for models that require this, while a `Multiclass` variable can only be one-hot encoded.

"""

# ╔═╡ b80f22e2-8764-11eb-1965-c13d87724698
html"""

  back to TOC

"""

# ╔═╡ a1223a20-8754-11eb-01c3-f5ae2dab2d10
md"""
### Two-dimensional data

Whenever it makes sense, MLJ Models generally expect two-dimensional data to be *tabular*. All the tabular formats implementing the  [Tables.jl API](https://juliadata.github.io/Tables.jl/stable/) (see this
 [list](https://github.com/JuliaData/Tables.jl/blob/master/INTEGRATIONS.md))
 have a scientific type of `Table` and can be used with such models. $(html"
")
The simplest example of a table is the julia native *column table*, which is just a named tuple of equal-length vectors:

"""

# ╔═╡ a1017da8-8754-11eb-1439-8da18a5b4fb2
begin
  column_table = (h=height, e=exam_mark_, t=time)
  #-
  scitype(column_table)
end

# ╔═╡ 8e76b9fa-8754-11eb-1d19-37ef9df3cb02
md"""
Notice the `Table{K}` type parameter `K` encodes the scientific types of the columns. (This is useful when comparing table scitypes with `<:`). To inspect the individual column scitypes, we use the `schema` method instead:

"""

# ╔═╡ e5835d02-8754-11eb-31e3-117515bbf480
schema(column_table)

# ╔═╡ e56c4c7a-8754-11eb-276d-27b1b96cfcc2
begin
  row_table = [(a=1, b=3.4),
    (a=2, b=4.5),
    (a=3, b=5.6)]
  schema(row_table)
end

# ╔═╡ e5501170-8754-11eb-036e-e346dd918f26
begin
  # using DataFrames
  df = DataFrames.DataFrame(column_table)
  #-
  schema(df)
end

# ╔═╡ e53b6164-8754-11eb-1ff9-67fff6d1aa66
begin
  # using CSV
  const DIR = "."
  file = CSV.File(joinpath(DIR, "data", "horse.csv"));
  schema(file) # (triggers a file read)
end

# ╔═╡ e5207a7a-8754-11eb-23a4-5547a267fcf1
begin
  matrix_table = MLJ.table(rand(2,3))
  schema(matrix_table)
  # The matrix is *not* copied, only wrapped.
end

# ╔═╡ c7fbf07c-8764-11eb-2c7c-577d1599693a
html"""

  back to TOC

"""

# ╔═╡ e5075f04-8754-11eb-25e7-0d6208bc1994
md"""
### Fixing scientific types in tabular data

To show how we can correct the scientific types of data in tables, we introduce a cleaned up version of the UCI Horse Colic Data Set (the cleaning work-flow is described [here](https://alan-turing-institute.github.io/DataScienceTutorials.jl/end-to-end/horse/#dealing_with_missing_values))
"""

# ╔═╡ 05353384-8762-11eb-00f0-b107afd1e3c2
begin
  hfile = CSV.File(joinpath(DIR, "data", "horse.csv"));
  horse = DataFrames.DataFrame(hfile); # convert to data frame without copying columns
  first(horse, 4)
end

# ╔═╡ 051e501a-8762-11eb-38f6-2bf36faca359
md"""
From [the UCI
docs](http://archive.ics.uci.edu/ml/datasets/Horse+Colic) we can
surmise how each variable ought to be interpreted (a step in our
work-flow that cannot reliably be left to the computer):

variable                    | scientific type (interpretation)
----------------------------|-----------------------------------
`:surgery`                  | Multiclass
`:age`                      | Multiclass
`:rectal_temperature`       | Continuous
`:pulse`                    | Continuous
`:respiratory_rate`         | Continuous
`:temperature_extremities`  | OrderedFactor
`:mucous_membranes`         | Multiclass
`:capillary_refill_time`    | Multiclass
`:pain`                     | OrderedFactor
`:peristalsis`              | OrderedFactor
`:abdominal_distension`     | OrderedFactor
`:packed_cell_volume`       | Continuous
`:total_protein`            | Continuous
`:outcome`                  | Multiclass
`:surgical_lesion`          | OrderedFactor
`:cp_data`                  | Multiclass
"""

# ╔═╡ 0502deca-8762-11eb-151a-29f6fa8801ef
md"""
Let's see how MLJ will actually interpret the data, as it is currently encoded:
"""

# ╔═╡ 04ec04d4-8762-11eb-126d-a90a7d326b72
schema(horse)

# ╔═╡ 030776da-8762-11eb-05ef-59c20464000a
md"""
As a first correction step, we can get MLJ to "guess" the appropriate fix, using the `autotype` method:
"""

# ╔═╡ 23e4b28c-8762-11eb-326b-dba661ca1e8d
autotype(horse)

# ╔═╡ 23ca0dec-8762-11eb-2bc6-d7d4e43255ec
md"""
Okay, this is not perfect, but a step in the right direction, which we implement like this:
"""

# ╔═╡ 23af76a8-8762-11eb-3ab4-03edc83b8310
begin
  coerce!(horse, autotype(horse));
  schema(horse)
end

# ╔═╡ 23952e56-8762-11eb-11ab-9b0f41f746fc
md"""
All remaining `Count` data should be `Continuous`:
"""

# ╔═╡ 46679cfc-8762-11eb-001b-bfd5fbb0d55b
begin
  coerce!(horse, Count => Continuous);
  schema(horse)
end

# ╔═╡ 4654e77e-8762-11eb-1c4f-5594e9a167cc
md"""
 We'll correct the remaining truant entries manually:
"""

# ╔═╡ 4636460c-8762-11eb-3612-4550488a4826
begin
  coerce!(horse,
        :surgery               => Multiclass,
        :age                   => Multiclass,
        :mucous_membranes      => Multiclass,
        :capillary_refill_time => Multiclass,
        :outcome               => Multiclass,
        :cp_data               => Multiclass);
  schema(horse)
end

# ╔═╡ de71ad74-8764-11eb-3d40-99f69a188626
html"""

  back to TOC

"""

# ╔═╡ 23819eec-8762-11eb-2aa3-65dd52774eb6
md"""
### Resources for Part 1

- From the MLJ manual:
   - [A preview of data type specification in
  MLJ](https://alan-turing-institute.github.io/MLJ.jl/dev/getting_started/#A-preview-of-data-type-specification-in-MLJ-1)
   - [Data containers and scientific types](https://alan-turing-institute.github.io/MLJ.jl/dev/getting_started/#Data-containers-and-scientific-types-1)
   - [Working with Categorical Data](https://alan-turing-institute.github.io/MLJ.jl/dev/working_with_categorical_data/)
- [Summary](https://alan-turing-institute.github.io/MLJScientificTypes.jl/dev/#Summary-of-the-MLJ-convention-1) of the MLJ convention for representing scientific types
- [MLJScientificTypes.jl](https://alan-turing-institute.github.io/MLJScientificTypes.jl/dev/)
- From Data Science Tutorials:
    - [Data interpretation: Scientific Types](https://alan-turing-institute.github.io/DataScienceTutorials.jl/data/scitype/)
    - [Horse colic data](https://alan-turing-institute.github.io/DataScienceTutorials.jl/end-to-end/horse/)
- [UCI Horse Colic Data Set](http://archive.ics.uci.edu/ml/datasets/Horse+Colic)
"""

# ╔═╡ c037d4ce-8763-11eb-07c9-53b58dc26075
html"""

  back to TOC

"""

# ╔═╡ 5ac50c98-8762-11eb-0577-ad8d0739b47b
md"""
### Exercises for Part 1

#### Exercise 1
Try to guess how each code snippet below will evaluate:
"""

# ╔═╡ 5a913224-8762-11eb-2655-35a1ef0a1080
@test  scitype(42) == Count

# ╔═╡ 5a74975e-8762-11eb-0dce-59bb2cc49162
begin
  questions = ["who", "why", "what", "when"]
  @test scitype(questions) == AbstractVector{Textual}
  @test elscitype(questions) == Textual
end

# ╔═╡ 7677f504-8762-11eb-2b07-cfb36a99343b
begin
  A = rand(2, 3)
  @test scitype(A) == AbstractMatrix{Continuous}
  @test elscitype(A) == Continuous
end

# ╔═╡ 7661f966-8762-11eb-28fd-d36d8e4ad3d2
begin
  using SparseArrays
  Asparse = sparse(A)
  @test scitype(A) == AbstractMatrix{Continuous}
end

# ╔═╡ 76461336-8762-11eb-3074-e9fdbc9b1c5a
begin
  # using CategoricalArrays
  C1 = categorical(A)

  @test scitype(C1) == AbstractMatrix{Multiclass{6}}
  @test elscitype(C1) == Multiclass{6}
end

# ╔═╡ 762e61f8-8762-11eb-34b3-377a8a390bc6
begin
  v = [1, 2, missing, 4]
  @test scitype(v) == AbstractVector{Union{Missing, Count}}
  @test elscitype(v) == Union{Missing, Count}
  @test scitype(v[1:2]) == AbstractVector{Union{Missing, Count}}
end

# ╔═╡ 7612ee66-8762-11eb-2c31-d9e4afb204b9
md"""
#### Exercise 2

Coerce the following vector to make MLJ recognize it as a vector of  ordered factors (with an appropriate ordering):
"""

# ╔═╡ 8edd43ec-8762-11eb-28a0-2f52edaaa3a6
quality = ["good", "poor", "poor", "excellent", missing, "good", "excellent"]

# ╔═╡ 8ec5f76e-8762-11eb-0d91-6ddc349e67b5
begin
  nquality = coerce(quality, OrderedFactor)
  levels!(nquality, ["poor", "good", "excellent"])

  @test scitype(nquality) == AbstractVector{Union{Missing, OrderedFactor{3}}}
  @test nquality[1] > nquality[2] # as nquality[1] == "good" and nquality[2] == "poor"
end

# ╔═╡ 8eab3c6c-8762-11eb-1b3a-134e39a0621a
md"""
#### Exercise 3 (fixing scitypes in a table)

Fix the scitypes for the [House Prices in King
County](https://mlr3gallery.mlr-org.com/posts/2020-01-30-house-prices-in-king-county/)
dataset:
"""

# ╔═╡ 8e8fa024-8762-11eb-13f1-7b9f34905d87
begin
  hofile = CSV.File(joinpath(DIR, "data", "house.csv"));
  house = DataFrames.DataFrame(hofile); # convert to data frame without copying columns
  first(house, 4)
end

# ╔═╡ e200acee-8762-11eb-2ad3-0f7693b1377d
md"""
(Two features in the original data set have been deemed uninformative
and dropped, namely `:id` and `:date`. The original feature
`:yr_renovated` has been replaced by the `Bool` feature `is_renovated`.)
"""

# ╔═╡ e402aa38-8762-11eb-021e-e30f60a98bfa
@show schema(house)

# ╔═╡ fd5bc7bc-8762-11eb-19fb-9f28925654c6
autotype(house)

# ╔═╡ Cell order:
# ╟─e3664c92-8753-11eb-325b-95ca5bcab046
# ╠═18f8736c-8754-11eb-3a0c-c729d6c6b0d8
# ╟─18e124b4-8754-11eb-07d4-a9eee0b8ba3b
# ╟─992cc480-8785-11eb-0572-c33c1c988efc
# ╟─53bc7eb6-8764-11eb-326d-41138236d179
# ╟─9dbf597a-8764-11eb-01e6-67bac4f668b5
# ╟─18c455fa-8754-11eb-2683-9dfe2784e339
# ╠═18acb3c8-8754-11eb-2d0a-5d82e2b47551
# ╟─185c6c10-8754-11eb-16a3-4b2439aa1d84
# ╠═7f3e28ec-8754-11eb-3c7e-e71dc9b2ca4e
# ╠═841c8930-8754-11eb-1596-0535d54e7558
# ╟─88b066ce-8754-11eb-1ec8-61fc466493d5
# ╠═8edfcc74-8754-11eb-3ac9-913d6db5d0c7
# ╟─8ec4a37c-8754-11eb-05c6-e9d5aa55479f
# ╠═8e9259d0-8754-11eb-01ee-cd92c6663059
# ╠═a152aad4-8754-11eb-33d1-cf4df8cd76d2
# ╟─a134e1fc-8754-11eb-1664-1d96b05f7ca4
# ╠═ac1420c4-8754-11eb-3751-4df7eff8ada4
# ╠═abfabfc6-8754-11eb-0a32-93a7e254063d
# ╟─bd7c20be-8754-11eb-07a1-e182ca10e932
# ╠═abdefbf4-8754-11eb-1bc5-4950210910d7
# ╟─abc65bbe-8754-11eb-325f-d98618a30eb8
# ╟─b80f22e2-8764-11eb-1965-c13d87724698
# ╟─a1223a20-8754-11eb-01c3-f5ae2dab2d10
# ╠═a1017da8-8754-11eb-1439-8da18a5b4fb2
# ╟─8e76b9fa-8754-11eb-1d19-37ef9df3cb02
# ╠═e5835d02-8754-11eb-31e3-117515bbf480
# ╠═e56c4c7a-8754-11eb-276d-27b1b96cfcc2
# ╠═e5501170-8754-11eb-036e-e346dd918f26
# ╠═e53b6164-8754-11eb-1ff9-67fff6d1aa66
# ╠═e5207a7a-8754-11eb-23a4-5547a267fcf1
# ╟─c7fbf07c-8764-11eb-2c7c-577d1599693a
# ╟─e5075f04-8754-11eb-25e7-0d6208bc1994
# ╠═05353384-8762-11eb-00f0-b107afd1e3c2
# ╟─051e501a-8762-11eb-38f6-2bf36faca359
# ╟─0502deca-8762-11eb-151a-29f6fa8801ef
# ╠═04ec04d4-8762-11eb-126d-a90a7d326b72
# ╟─030776da-8762-11eb-05ef-59c20464000a
# ╠═23e4b28c-8762-11eb-326b-dba661ca1e8d
# ╟─23ca0dec-8762-11eb-2bc6-d7d4e43255ec
# ╠═23af76a8-8762-11eb-3ab4-03edc83b8310
# ╟─23952e56-8762-11eb-11ab-9b0f41f746fc
# ╠═46679cfc-8762-11eb-001b-bfd5fbb0d55b
# ╟─4654e77e-8762-11eb-1c4f-5594e9a167cc
# ╠═4636460c-8762-11eb-3612-4550488a4826
# ╠═de71ad74-8764-11eb-3d40-99f69a188626
# ╟─23819eec-8762-11eb-2aa3-65dd52774eb6
# ╟─c037d4ce-8763-11eb-07c9-53b58dc26075
# ╟─5ac50c98-8762-11eb-0577-ad8d0739b47b
# ╠═5aa90fa0-8762-11eb-18a8-a5a45edfd5a7
# ╠═5a913224-8762-11eb-2655-35a1ef0a1080
# ╠═5a74975e-8762-11eb-0dce-59bb2cc49162
# ╠═7677f504-8762-11eb-2b07-cfb36a99343b
# ╠═7661f966-8762-11eb-28fd-d36d8e4ad3d2
# ╠═76461336-8762-11eb-3074-e9fdbc9b1c5a
# ╠═762e61f8-8762-11eb-34b3-377a8a390bc6
# ╟─7612ee66-8762-11eb-2c31-d9e4afb204b9
# ╠═8edd43ec-8762-11eb-28a0-2f52edaaa3a6
# ╠═8ec5f76e-8762-11eb-0d91-6ddc349e67b5
# ╟─8eab3c6c-8762-11eb-1b3a-134e39a0621a
# ╠═8e8fa024-8762-11eb-13f1-7b9f34905d87
# ╟─e200acee-8762-11eb-2ad3-0f7693b1377d
# ╠═e402aa38-8762-11eb-021e-e30f60a98bfa
# ╠═fd5bc7bc-8762-11eb-19fb-9f28925654c6

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:
__precompile__(true)

"Helper to generate new julia packages."
module NewPkg

export Generate

include("generate.jl")

doc"""
    generate(pkg, description)

Generate a new package named `pkg` with `MIT "Expat" License`.
Generate creates a git repo at `Pkg.dir(pkg)` for the package with an initial file structure.

Arguments:

* `pkg` - a name for the new julia package.
* `description` -  a short description of the new julia package.

Keyword parameters:

* `authors` - a string or array of author names, the final default value will be the `package GitRepo user.name`.
* `authors_url` - a webpage or email address related to the copyright or authors, the final default value will be
  `https://github.com/package GitRepo user.name/`. To skip it set it to `"NONE"`.
* `path` - a location where the package will be generated, the default location is `Pkg.dir()`
* `github_name` - github organisation or user name, the final default value will be `package GitRepo user.name`. To
    skip it set it to `"NONE"`. This is used to set any `package repo remote github url` and any mkdocs: `site_url`,
    `repo_url` etc...
* `github_ssh` - if true sets the package repo remote github url to "ssh version" otherwise the "https", the default
    value is `true`.
* `mkdocs` - enables generation of a `mkdocs.yml` configuration file as well as an initial `$pkg/docs` folder, the
    default value is `true`. <br />
    **Note:** initial `site_url`, `repo_url` ... are generated for github.
"""
generate(pkg::AbstractString, description::AbstractString;
            authors::Union{AbstractString, Array} = [],
            authors_url::AbstractString           = "",
            path::AbstractString                  = Pkg.Dir.path(),
            github_name::AbstractString           = "",
            github_ssh::Bool                      = true,
            mkdocs::Bool                          = true
        ) = Generate.package(pkg, description;
                                authors     = authors,
                                authors_url = authors_url,
                                path        = path,
                                github_name = github_name,
                                github_ssh  = github_ssh,
                                mkdocs      = mkdocs
                            )

"""
    config(force::Bool=false)
Interactive configuration of the development environment.

NewPkg operations require `git` minimum configuration that keeps user signature (user.name & user.email).
"""
function config(force::Bool=false)
    # setup global git configuration
    cfg = LibGit2.GitConfig(LibGit2.Consts.CONFIG_LEVEL_GLOBAL)
    try
        println("Julia NewPkg configuration:")

        username = LibGit2.get(cfg, "user.name", "")
        if isempty(username) || force
            username = LibGit2.prompt("Enter user name", default=username)
            LibGit2.set!(cfg, "user.name", username)
        else
            println("User name: $username")
        end

        useremail = LibGit2.get(cfg, "user.email", "")
        if isempty(useremail) || force
            useremail = LibGit2.prompt("Enter user email", default=useremail)
            LibGit2.set!(cfg, "user.email", useremail)
        else
            println("User email: $useremail")
        end
    finally
        finalize(cfg)
    end
    lowercase(LibGit2.prompt("Do you want to change this configuration?", default="N")) == "y" && config(true)
    return
end

function __init__()
    # Check if git configuration exists
    cfg = LibGit2.GitConfig(LibGit2.Consts.CONFIG_LEVEL_GLOBAL)
    try
        username = LibGit2.get(cfg, "user.name", "")
        if isempty(username)
            warn("Julia NewPkg is not configured. Please, run `NewPkg.config()` before performing any operations.")
        end
    finally
        finalize(cfg)
    end
end

end # module


text:
using FeedbackParticleFilters, Distributions, LinearAlgebra, Random, PDMats, FillArrays

println("Testing DiffusionStateModel.jl:")

@testset "DiffusionStateModel.jl" begin

    f(x::AbstractVector) = [x[1]*x[2], x[1]-x[3], x[1]+x[2]*x[3]]
    g(x::AbstractVector) = [x[1] x[2]; x[2] x[3]; 0. 1.]
    init                 = MvNormal(LinearAlgebra.Diagonal(FillArrays.Fill(1., 3)))

    print("  inner constructor")
    model = DiffusionStateModel(f, g, init)
    print("-")
    @test_throws BoundsError DiffusionStateModel(f, g, [1., 2.]) # wrong length of initial condition
    print("-")
    @test_throws ErrorException DiffusionStateModel(f, g, [1., 2., 3., 4.]) # wrong length of initial condition
    print("-")
    g2(x::AbstractVector) = [x[1] x[2]; x[2] x[3]]
    @test_throws ErrorException DiffusionStateModel(f, g2, init) # wrong size of diffusion term
    print("-")
    @test_throws ErrorException DiffusionStateModel(f, g, [0, 0, 0]) # wrong type of initial condition
    println("DONE")

    print("  method drift")
    print("-")
    @test drift(model)([1., 2., 3.]) == [2., -2., 7.]
    print("-")
    @test hasmethod(drift(model), Tuple{AbstractMatrix})
    println("DONE")

    print("  method diffusion")
    print("-")
    @test diffusion(model)([1., 2., 3.]) == [1. 2.; 2. 3.; 0. 1.]
    print("-")
    @test hasmethod(diffusion(model), Tuple{AbstractMatrix})
    println("DONE")

    print("  method initial_condition")
    print("-")
    @test initial_condition(model) == model.init
    println("DONE")

    print("  method state_dim")
    print("-")
    @test state_dim(model) == 3
    println("DONE")

    print("  method noise_dim")
    print("-")
    @test noise_dim(model) == 2
    println("DONE")

    print("  method initialize")
    print("-")
    print("-")
    model2 = DiffusionStateModel(f, g, [1., 2., 3.])
    @test initialize(model2) == [1., 2., 3.]
    println("DONE")

    print("  method drift_function")
    print("-")
    @test drift_function(model)([1., 2., 3.]) == drift(model)([1., 2., 3.])
    println("DONE")

    print("  method diffusion_function")
    print("-")
    @test diffusion_function(model)([1., 2., 3.]) == diffusion(model)([1., 2., 3.])
    println("DONE")

    print("  calling instance of DiffusionStateModel")
    print("-")
    x0 = ones(3)
    @test model(x0, 0.01) != x0
    print("-")
    xx = hcat(ones(3), zeros(3))
    @test model(xx, 0.01) != xx
    println("DONE")

    print("  method propagate!")
    print("-")
    x0 = ones(3)
    x1 = copy(x0)
    propagate!(x1, model, 0.01)
    @test x1 != x0
    print("-")
    xx0 = hcat(ones(3), zeros(3))
    xx1 = copy(xx0)
    propagate!(xx1, model, 0.01)
    @test xx1 != xx0
    println("DONE")

end; #DiffusionStateModel.jl

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
doc"""
    second_derivative(f::Function, x::AbstractFloat)

# Description

Evaluate the second derivative of `f` at `x` using finite differences. In
mathematical notation, we calculate,

```math
\frac{f(x + \epsilon) - 2 f(x) + f(x - \epsilon)}{\epsilon^2},
```

where ``\epsilon`` is chosen to be small enough to approximate the second
derivative, but not so small as to suffer from extreme numerical inaccuracy.

# Arguments

* `f::Function`: The function to be differentiated.
* `x::AbstractFloat`: The point at which to evaluate the derivative of `f`. Its
    type must implement `eps`.

# Returns

* `y::Real`: The second derivative of `f` evaluated at `x`.

# Examples

```jl
import FiniteDiff: second_derivative
y = second_derivative(sin, 1.0)
```
"""
@inline function second_derivative(f::Function, x::AbstractFloat)::Real
    ϵ = step_size(HessianMode(), x)
    return (f(x + ϵ) - 2 * f(x) + f(x - ϵ)) / (ϵ * ϵ)
end

doc"""
    second_derivative!(
        output::AbstractArray,
        f::Function,
        x::AbstractFloat,
    )

# Description

Evaluate the second derivative of `f` at `x` using finite differences. In
mathematical notation, we calculate,

```math
\frac{f(x + \epsilon) - 2 f(x) + f(x - \epsilon)}{\epsilon^2},
```

where ``\epsilon`` is chosen to be small enough to approximate the second
derivative, but not so small as to suffer from extreme numerical inaccuracy.

The first argument, `output`, will be mutated so that the value of the second
derivative is its first element.

# Arguments

* `output::AbstractArray`: An array whose first element will be mutated.
* `f::Function`: The function to be differentiated.
* `x::AbstractFloat`: The point at which to evaluate the derivative of `f`. Its
    type must implement `eps`.

# Returns

* `void::Nothing`: This function is called for its side-effects.

# Examples

```jl
import FiniteDiff: second_derivative
y = Array(Float64, 1)
second_derivative!(y, sin, 1.0)
```
"""
function second_derivative!(
    output::AbstractArray,
    f::Function,
    x::AbstractFloat,
)::Void
    output[1] = second_derivative(f, x)
    return
end

"""
    second_derivative(f::Function; mutates::Bool=false)

# Description

Construct a new function that will evaluate the second derivative of `f` at
any point `x`. A keyword argument specifies whether the resulting function
should be mutating or non-mutating.

# Arguments

* `f::Function`: The function to be differentiated.

# Keyword Arguments

* `mutates::Bool = false`: Determine whether the resulting function will mutate
    its inputs or will be a pure function.

# Returns

* `void::Nothing`: This function is called for its side-effects.

# Examples

```jl
import FiniteDiff: second_derivative
f′′ = second_derivative(sin)
f′′(1.0)
```
"""
function second_derivative(f::Function; mutates::Bool=false)::Function
    if mutates
        return (output, x) -> second_derivative!(output, f, x)
    else
        return (x, ) -> second_derivative(f, x)
    end
end


text:
function sample(model, input::Array{Float32,2};
    outdim::Int = 5)::Array{Float32,2}

    # get the size of the input -
    (R,C) = size(input)
    output_array = Array{Float32,2}(undef, R, outdim)

    for i ∈ 1:R
        output_sample = model(input[i,:])
        for j ∈ 1:outdim
            output_array[i,j] = output_sample[j]
        end
    end

    return output_array
end

function sample(model, input::Vector{Float32})::Vector{Float32}
    return model(input)
end

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
module PatientAllocationResults

using DataFrames
using LinearAlgebra
using Dates
using Distributions

function results_all(
		sent::Array{<:Real,3},
		beds::Array{<:Real,1},
		initial_patients::Array{<:Real,1},
		discharged_patients::Array{<:Real,2},
		admitted_patients::Array{<:Real,2},
		locations::Array{String,1},
		start_date::Date,
		los,
)
	N, _, T = size(sent)

	summary = results_summary(sent, beds, initial_patients, discharged_patients, admitted_patients, locations, start_date, los)
	complete = results_complete(sent, beds, initial_patients, discharged_patients, admitted_patients, locations, start_date, los)

	sent_matrix_table = results_sentmatrix_table(sent, locations)
	sent_matrix_vis = results_sentmatrix_vis(sent, initial_patients, admitted_patients, locations)

	sent_to = Dict(locations[i] => locations[row] for (i,row) in enumerate(eachrow(sum(sent, dims=3)[:,:,1] .> 0)))

	active, active_null = results_active_patients(sent, initial_patients, discharged_patients, admitted_patients, los)

	overflow = [max(0, active[i,t] - beds[i]) for i in 1:N, t in 1:T]
	load = [active[i,t] / beds[i] for i in 1:N, t in 1:T]

	overflow_null = [max(0, active_null[i,t] - beds[i]) for i in 1:N, t in 1:T]
	load_null = [active_null[i,t] / beds[i] for i in 1:N, t in 1:T]

	total_overflow = sum(overflow_null)
	avgload = sum(active_null) / sum(beds)

	return (
		total_overflow=total_overflow,
		average_load=avgload,
		summary_table=summary,
		complete_table=complete,
		sent_matrix_table=sent_matrix_table,
		sent_matrix_vis=sent_matrix_vis,
		sent_to=sent_to,
		active_patients=active,
		overflow=overflow,
		load=load,
		active_patients_nosent=active_null,
		overflow_nosent=overflow_null,
		load_nosent=load_null,
	)
end

function results_summary(
		sent::Array{<:Real,3},
		beds::Array{<:Real,1},
		initial_patients::Array{<:Real,1},
		discharged_patients::Array{<:Real,2},
		admitted_patients::Array{<:Real,2},
		locations::Array{String,1},
		start_date::Date,
		los,
)
	N, _, T = size(sent)

	active, active_null = results_active_patients(sent, initial_patients, discharged_patients, admitted_patients, los)

	overflow = [sum(max(0, active[i,t] - beds[i]) for t in 1:T) for i in 1:N]
	load = [sum(active[i,t] / beds[i] for t in 1:T)/T for i in 1:N]

	overflow_null = [sum(max(0, active_null[i,t] - beds[i]) for t in 1:T) for i in 1:N]
	load_null = [sum(active_null[i,t] / beds[i] for t in 1:T)/T for i in 1:N]

	summary = DataFrame(
		location=locations,
		total_sent=sum(sent, dims=[2,3])[:],
		total_received=sum(sent, dims=[1,3])[:],
		overflow=overflow,
		overall_load=load,
		overflow_nosent=overflow_null,
		overall_load_nosent=load_null,
	)

	return summary
end

function results_complete(
		sent::Array{<:Real,3},
		beds::Array{<:Real,1},
		initial_patients::Array{<:Real,1},
		discharged_patients::Array{<:Real,2},
		admitted_patients::Array{<:Real,2},
		locations::Array{String,1},
		start_date::Date,
		los,
)
	N, _, T = size(sent)

	active_patients, active_patients_null = results_active_patients(sent, initial_patients, discharged_patients, admitted_patients, los)

	overflow = [max(0, active_patients[i,t] - beds[i]) for i in 1:N, t in 1:T]
	load = [(active_patients[i,t] / beds[i]) for i in 1:N, t in 1:T]

	overflow_null = [max(0, active_patients_null[i,t] - beds[i]) for i in 1:N, t in 1:T]
	load_null = [(active_patients_null[i,t] / beds[i]) for i in 1:N, t in 1:T]

	outcomes = DataFrame()
	for (i,s) in enumerate(locations)
		single_state_outcome = DataFrame(
			location=fill(s, T),
			date=start_date .+ Day.(0:T-1),
			sent=sum(sent[i,:,:], dims=1)[:],
			received=sum(sent[:,i,:], dims=1)[:],
			new_patients=admitted_patients[i,:],
			active_patients=active_patients[i,:],
			active_patients_nosent=active_patients_null[i,:],
			capacity=fill(beds[i], T),
			overflow=overflow[i,:],
			load=load[i,:],
			overflow_nosent=overflow_null[i,:],
			load_nosent=load_null[i,:],
			sent_to=[sum(sent[i,:,t])>0 ? collect(zip(locations[sent[i,:,t] .> 0], sent[i,sent[i,:,t].>0,t])) : "[]" for t in 1:T],
			sent_from=[sum(sent[:,i,t])>0 ? collect(zip(locations[sent[:,i,t] .> 0], sent[sent[:,i,t].>0,i,t])) : "[]" for t in 1:T],
		)
		outcomes = vcat(outcomes, single_state_outcome)
	end

	return outcomes
end

function results_sentmatrix_table(sent::Array{<:Real,3}, locations::Array{String,1})
	sent_matrix = DataFrame(sum(sent, dims=3)[:,:,1])
	rename!(sent_matrix, Symbol.(locations))
	insertcols!(sent_matrix, 1, :state => locations)
	return sent_matrix
end

function results_sentmatrix_vis(sent::Array{<:Real,3}, initial_patients::Array{<:Real,1},
		admitted_patients::Array{<:Real,2}, locations::Array{String,1})
	selfedges = initial_patients + sum(admitted_patients, dims=2)[:] - sum(sent, dims=[2,3])[:]
	sent_vis_matrix = sum(sent, dims=3)[:,:,1] + diagm(selfedges)
	sent_vis_matrix = DataFrame(sent_vis_matrix)
	rename!(sent_vis_matrix, Symbol.(locations))
	return sent_vis_matrix
end

function results_active_patients(
		sent::Array{<:Real,3},
		initial_patients::Array{<:Real,1},
		discharged_patients::Array{<:Real,2},
		admitted_patients::Array{<:Real,2},
		los,
)
	N, _, T = size(sent)

	L = nothing
	if isa(los, Int)
		L = vcat(ones(Int, los), zeros(Int, T-los))
	elseif isa(los, Array{<:Real,1})
		if length(los) >= T
			L = los
		else
			L = vcat(los, zeros(Float64, T-length(los)))
		end
	elseif isa(los, Distribution)
		L = 1.0 .- cdf.(los, 0:T)
	else
		error("Invalid length of stay distribution")
	end

	active = [(
			initial_patients[i]
			- sum(discharged_patients[i,1:t])
			+ sum(L[t-t₁+1] * (
				admitted_patients[i,t₁]
				- sum(sent[i,:,t₁])
				+ sum(sent[:,i,t₁])
			) for t₁ in 1:t)
			+ sum(sent[i,:,t])
		) for i in 1:N, t in 1:T
	]
	active_null = [(
			initial_patients[i]
			- sum(discharged_patients[i,1:t])
			+ sum(L[t-t₁+1] * admitted_patients[i,t₁] for t₁ in 1:t)
		) for i in 1:N, t in 1:T
	]

	return active, active_null
end

end

max_stars_repo_name: flixpar/covid-resource-allocation

max_issues_repo_name: flixpar/covid-resource-allocation

max_forks_repo_name: flixpar/covid-resource-allocation


text:
# ---------------------------------------------------------------------------- #
#
#   writeTecplot.jl
#
#   Several functions to write Tecplot output
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    writeTecplot( flname::String, prob::Elas, mesh::Mesh3D,
                  uh::Array{Float64,3}, σh::Array{Float64,3};
                  σmises = Array{Float64,3}(0,0,0) )

Outputs tecplot data file to `flname` for 3D structural elasticity problems.
"""
function writeTecplotP1( flname::String, prob::luteos.Elas, mesh::Mesh3D,
                         uh::Matrix{Float64}, σh::Matrix{Float64}; σmises = Matrix{Float64}(0,0) )

    nelem = size( mesh.nodes, 3 )
    ndat  = 2*mesh.dim + mesh.dim^2
    if size(σmises,1) > 0
        ndat += 1
    end

    # Open file
    fid = open( flname, "w" )

    ### Write header
    @printf( fid, "TITLE = \"%s\"\n", prob.name )
    @printf( fid, "VARIABLES = \"X\", \"Y\", \"Z\"" )
    for jj in 1:mesh.dim
        @printf( fid, ", \"U%i\"",jj )
    end
    for ii in 1:mesh.dim, jj in 1:mesh.dim
        @printf( fid, ", \"<greek>s</greek>%i%i\"", ii, jj )
    end
    if size( σmises, 1 ) > 0
        @printf( fid, ", \"<greek>s</greek>mises\"" )
    end
    @printf( fid, "\n" )
    @printf( fid, "ZONE, DATAPACKING=POINT, NODES=%i, ELEMENTS=%i, ZONETYPE=FETETRAHEDRON\n", mesh.n, nelem )

    ### Write values
    for pp in 1:mesh.n
        if size( σmises, 1 ) == 0
            wrtval = [mesh.p[pp,:]; uh[pp,:]; σh[pp,:] ]
        else
            wrtval = [mesh.p[pp,:]; uh[pp,:]; σh[pp,:]; σmises[pp,:] ]
        end

    for ii in 1:ndat
      @printf(fid, "%16.15e\t", wrtval[ii] )
    end
    @printf(fid,"\n")

    end

    ### Write connectivity
    for ee in 1:nelem

        for kk in 1:(mesh.dim+1)
            @printf(fid, "%i\t", mesh.t[ee,kk] )
        end
        @printf(fid,"\n")

    end

    close(fid)

end

"""
    writeTecplot( flname::String, prob::Elas, mesh::Mesh3D,
                  uh::Array{Float64,3}, σh::Array{Float64,3};
                  σmises = Array{Float64,3}(0,0,0) )

Outputs tecplot data file to `flname` for 3D structural elasticity problems.
"""
function writeTecplotP1( flname::String, prob::luteos.Problem, mesh::Mesh3D, solution::Matrix{Float64} )

    nelem   = size( mesh.nodes, 3 )
    nnodes  = size( mesh.nodes, 1 )
    ndat    = size( solution, 2 )

    @assert size(solution,1) == mesh.n

    # Open file
    fid = open( flname, "w" )

    ### Write header
    @printf( fid, "TITLE = \"%s\"\n", prob.name )
    @printf( fid, "VARIABLES = \"X\", \"Y\", \"Z\"" )
    for jj in 1:ndat
        @printf( fid, ", \"U%i\"",jj )
    end
    @printf( fid, "\n" )
    @printf( fid, "ZONE, DATAPACKING=POINT, NODES=%i, ELEMENTS=%i, ZONETYPE=FETETRAHEDRON\n", mesh.n, nelem )

    ### Write values
    for pp in 1:mesh.n

        wrtval = [mesh.p[pp,:] solution[pp,:] ]

        for ii in 1:ndat+3
            @printf(fid, "%16.15e\t", wrtval[ii] )
        end
        @printf(fid,"\n")

    end

    ### Write connectivity
    for ee in 1:nelem

        for kk in 1:(mesh.dim+1)
            @printf(fid, "%i\t", mesh.t[ee,kk] )
        end
        @printf(fid,"\n")

    end

    close(fid)

end

"""
    writeTecplot( flname::String, prob::Elas, mesh::Mesh3D,
                  uh::Array{Float64,3}, σh::Array{Float64,3};
                  σmises = Array{Float64,3}(0,0,0) )

Outputs tecplot data file to `flname` for 3D structural elasticity problems.
"""
function writeTecplotP1Defl( flname::String, prob::luteos.Problem, mesh::Mesh3D, solution::Matrix{Float64} )

  nelem   = size( mesh.nodes, 3 )
  nnodes  = size( mesh.nodes, 1 )
  ndat    = size( solution, 2 )

  @assert size(solution,1) == mesh.n

  # defl = solution ./ maximum(solution) * maximum( mesh.p ) * 0.25
  defl = solution

  # Open file
  fid = open( flname, "w" )

  ### Write header
  @printf( fid, "TITLE = \"%s\"\n", prob.name )
  @printf( fid, "VARIABLES = \"X\", \"Y\", \"Z\"" )
  for jj in 1:ndat
    @printf( fid, ", \"U%i\"",jj )
  end
  @printf( fid, "\n" )
  @printf( fid, "ZONE, DATAPACKING=POINT, NODES=%i, ELEMENTS=%i, ZONETYPE=FETETRAHEDRON\n",
    mesh.n, nelem )

  ### Write values
  for pp in 1:mesh.n

    wrtval = [mesh.p[pp,:]+defl[pp,:] solution[pp,:] ]

    for ii in 1:ndat+3
      @printf(fid, "%16.15e\t", wrtval[ii] )
    end
    @printf(fid,"\n")

  end

  ### Write connectivity
  for ee in 1:nelem

    for kk in 1:(mesh.dim+1)
      @printf(fid, "%i\t", mesh.t[ee,kk] )
    end
    @printf(fid,"\n")

  end

  close(fid)

end

function convertToNodal( mesh::Mesh3D, prob::luteos.Elas, uh::Array{Float64,3}, σh::Array{Float64,3} )

    @assert size(mesh.nodes,1) == 4

    nelem = size( mesh.nodes, 3 )

    ncounts = fill( 0, mesh.n )

    σhmises = luteos.compMises3D( σh )

    uhvec    = fill( 0., mesh.n, mesh.dim )
    σhvec    = fill( 0., mesh.n, mesh.dim^2 )
    σhmisvec = fill( 0., mesh.n, 1 )

    for kk in 1:nelem

        indsp = mesh.t[kk,:]

        ncounts[indsp] .+= 1

        for ii in 1:(mesh.dim+1)
            for jj in 1:mesh.dim
                uhvec[indsp[ii],jj] += uh[ii,jj,kk]
            end
            for jj in 1:mesh.dim^2
                σhvec[indsp[ii],jj] += σh[ii,jj,kk]
            end
            σhmisvec[indsp[ii],1] += σhmises[ii,1,kk]
        end

    end

    for jj in 1:mesh.n
        uhvec[jj,:]   ./= ncounts[jj]
        σhvec[jj,:]   ./= ncounts[jj]
        σhmisvec[jj,:] ./= ncounts[jj]
    end

    return uhvec, σhvec, σhmisvec

end

# TODO: Need to add surface output


text:
mutable struct FiniteSignalType{T}
    s0::Vector{T} # default signal.
    a::Float64
    b::Float64

    s::Vector{T} # buffer.
    #a::Float64
    #b::Float64
    a_ind
    b_ind
end

function FiniteSignalType(s, a, b)
    return FiniteSignalType(copy(s), a, b, copy(s), 1, length(s))
end

function resetsignal!(S::FiniteSignalType{T}) where T

    S.a_ind = 1
    S.b_ind = length(S.s0)
    S.s[:] = S.s0

    return nothing
end

# positive is delay, or shift towards more positive.
function applyintegershift!(S::FiniteSignalType{T}, d::Int) where T
    #
    S.a_ind = S.a_ind-d
    S.b_ind = S.b_ind-d

    return nothing
end

function evalfinitesignal(s::Vector{T}, i::Int)::T where T
    if 1 <= i <= length(s)
        return s[i]
    end

    return zero(T)
end

function evalsignal(S::FiniteSignalType{T}, n::Int)::T where T

    #ind_t_Float = convertcompactdomain(t, S.a, S.b, S.a_ind, S.b_ind)
    #ind_t_Float = round(Int, ind_t_Float)
    #@assert 1 <= n <= length(S.s)

    return evalfinitesignal(S.s, S.a_ind + n-1)
end

#  https://lcondat.github.io/publis/condat_icip08_rotations.pdf
function getthiriancoeffs!(b::Vector{T}, τ::T) where T <: Real

    #resize!(b, N)
    N = length(b)
    for k = 1:N
        b[k] = (-1)^k * binomial(N,k) * prod( (τ-n)/(τ-n-k) for n = 0:N)
    end

    return nothing
end

# https://lcondat.github.io/publis/condat_icip08_rotations.pdf
# τ ∈ [0, 0.5]
# mutates s and b_buffer.
function applydelaythirian!(b_buffer::Vector{T}, s::Vector{T}, τ::T) where T <: Real

    getthiriancoeffs!(b_buffer, τ)

    N = length(b_buffer)
    M = length(s)

    for i = M-1:-1:0
        for k = 1:N

            #s[i+1] += b_buffer[k]*(s[i-k+1]-s[i+k+1])

            tmp = evalfinitesignal(s, i-k+1) - evalfinitesignal(s, i+k+1)

            s[i+1] += b_buffer[k]*tmp
        end
    end

    return nothing
end

function applyadvancethirian!(b_buffer::Vector{T}, s::Vector{T}, τ::T) where T <: Real

    getthiriancoeffs!(b_buffer, τ)

    N = length(b_buffer)
    M = length(s)

    #reverse!(s)
    for i = M-1:-1:0
        for k = 1:N

            # tmp = evalfinitesignal(s, i-k+1) - evalfinitesignal(s, i+k+1)
            # s[i+1] += b_buffer[k]*tmp

            # reverse: if we have the ORDERED lists:
            # U = {1, 2, ..., M}, V = {M, M-1, ... , 1}, then
            # j = M-i+1, where i is a given element from U, and j is the
            #   corresponding element in V.

            j = M-i #j = M-(i+1)+1
            j1 = j+k #j1 = M-(i-k+1)+1
            j2 = j-k #j2 = M-(i+k+1)+1
            tmp = evalfinitesignal(s, j1) - evalfinitesignal(s, j2)
            s[j] += b_buffer[k]*tmp
        end
    end
    #reverse!(s)

    return nothing
end

function applythirian!(b_buffer)

    if τ > 0.5
        return nothing # shift, then apply thirian.
    end

    return nothing # applythirian!()
end

function applyshift!(S::FiniteSignalType{T}, shift::T, b_thirian) where T

    d = round(Int, shift)
    τ = shift-d
    println(τ)

    applyintegershift!(S, d)
    if τ > 0
        applydelaythirian!(b_thirian, S.s, τ)
    else
        applyadvancethirian!(b_thirian, S.s, abs(τ))
    end

    return nothing
end

# other interesting articles about fractional delay: https://www.dsprelated.com/showarticle/1327.php


text:
using Test

function max_(arr)
    if length(arr) == 1
        return arr[1]
    end
    sub_max = max_(arr[2:end])
    return arr[1] > sub_max ? arr[1] : sub_max
end

@test max_([1,5,8,20,9,10]) == 20


text:
"""    PyABOD(n_neighbors = 5,
              method = "fast")
$(make_docs_link("abod"))
"""
@pymodel mutable struct PyABOD <: UnsupervisedDetector
    n_neighbors::Integer = 5::(_ > 0)
    method::String = "fast"::(_ in ("default", "fast"))
end

"""    PyCBLOF(n_clusters = 8,
               alpha = 0.9,
               beta = 5,
               use_weights = false,
               random_state = nothing,
               n_jobs = 1)
$(make_docs_link("cblof"))
"""
@pymodel mutable struct PyCBLOF <: UnsupervisedDetector
    n_clusters::Integer = 8::(_ > 1)
    alpha::Real = 0.9::(0.5 < _ < 1)
    beta::Real = 5::(_ > 1)
    use_weights::Bool = false
    random_state::Union{Nothing, Integer} = nothing
    n_jobs::Integer = 1::(_ >= -1)
end

"""    PyCOF(n_neighbors = 5)
$(make_docs_link("cof"))
"""
@pymodel mutable struct PyCOF <: UnsupervisedDetector
    n_neighbors::Integer = 5::(_ > 0)
end

"""    PyCOPOD()
$(make_docs_link("copod"))
"""
@pymodel mutable struct PyCOPOD <: UnsupervisedDetector end

"""    PyHBOS(n_bins = 10,
              alpha = 0.1,
              tol = 0.5)
$(make_docs_link("hbos"))
"""
@pymodel mutable struct PyHBOS <: UnsupervisedDetector
    n_bins::Integer = 10::(_ > 1)
    alpha::Real = 0.1::(0 < _ < 1)
    tol::Real = 0.5::(0 < _ < 1)
end

"""    PyIForest(n_estimators = 100,
                 max_samples = "auto",
                 max_features = 1.0
                 bootstrap = false,
                 behaviour = "new",
                 random_state = nothing,
                 verbose = 0,
                 n_jobs = 1)
$(make_docs_link("iforest"))
"""
@pymodel mutable struct PyIForest <: UnsupervisedDetector
    n_estimators::Integer = 100::(_ > 0)
    max_samples::Union{String, Real} = "auto"
    max_features::Real = 1.0
    bootstrap::Bool = false
    behaviour::String = "new"
    random_state::Union{Nothing, Integer} = nothing
    verbose::Integer = 0::(0 <= _ <= 2)
    n_jobs::Integer = 1::(_ >= -1)
end

"""    PyKNN(n_neighbors = 5,
             method = "largest",
             radius = 1.0,
             algorithm = "auto",
             leaf_size = 30,
             metric = "minkowski",
             p = 2,
             metric_params = nothing,
             n_jobs = 1)
$(make_docs_link("knn"))
"""
@pymodel mutable struct PyKNN <: UnsupervisedDetector
    n_neighbors::Integer = 5::(_ > 0)
    method::String = "largest"::(_ in ("largest", "mean", "median"))
    radius::Real = 1.0
    algorithm::String = "auto"::(_ in ("auto", "ball_tree", "kd_tree", "brute"))
    leaf_size::Integer = 30::(_ > 0)
    metric::String = "minkowski"::(_ in ("cityblock", "cosine", "euclidean", "l1", "l2", "manhatten", "braycurtis", "canberra", "chebyshev", "correlation", "dice", "hamming", "jaccard", "kulsinski", "mahalanobis", "matching", "minkowski", "rogerstanimoto", "russellrao", "seuclidean", "sokalmichener", "sokalsneath", "sqeuclidean", "yule"))
    p::Union{Nothing, Integer} = 2
    metric_params::Union{Nothing, Any} = nothing
    n_jobs::Integer = 1::(_ >= -1)
end

"""    PyLMDD(n_iter = 50,
              dis_measure = "aad",
              random_state = nothing)
$(make_docs_link("lmdd"))
"""
@pymodel mutable struct PyLMDD <: UnsupervisedDetector
    n_iter::Integer = 50::(_ > 0)
    dis_measure::String = "aad"::(_ in ("aad", "var", "iqr"))
    random_state::Union{Nothing, Integer} = nothing
end

"""    PyLODA(n_bins = 10,
              n_random_cuts = 100)
$(make_docs_link("loda"))
"""
@pymodel mutable struct PyLODA <: UnsupervisedDetector
    n_bins::Integer = 10::(_ > 1)
    n_random_cuts::Integer = 100::(_ > 0)
end

"""    PyLOF(n_neighbors = 5,
             method = "largest",
             algorithm = "auto",
             leaf_size = 30,
             metric = "minkowski",
             p = 2,
             metric_params = nothing,
             n_jobs = 1)
$(make_docs_link("lof"))
"""
@pymodel mutable struct PyLOF <: UnsupervisedDetector
    n_neighbors::Integer = 5::(_ > 0)
    algorithm::String = "auto"::(_ in ("auto", "ball_tree", "kd_tree", "brute"))
    leaf_size::Integer = 30::(_ > 0)
    metric::String = "minkowski"::(_ in ("cityblock", "cosine", "euclidean", "l1", "l2", "manhatten", "braycurtis", "canberra", "chebyshev", "correlation", "dice", "hamming", "jaccard", "kulsinski", "mahalanobis", "matching", "minkowski", "rogerstanimoto", "russellrao", "seuclidean", "sokalmichener", "sokalsneath", "sqeuclidean", "yule"))
    p::Union{Nothing, Integer} = 2
    metric_params::Union{Nothing, Any} = nothing
    n_jobs::Integer = 1::(_ >= -1)
end

"""    PyLOCI(alpha = 0.5,
              k = 3)
$(make_docs_link("loci"))
"""
@pymodel mutable struct PyLOCI <: UnsupervisedDetector
    alpha::Real = 0.5::(0 < _ < 1)
    k::Real = 3::(_ > 0)
end

"""    PyMCD(store_precision = true,
             assume_centered = false,
             support_fraction = nothing,
             random_state = nothing)
$(make_docs_link("mcd"))
"""
@pymodel mutable struct PyMCD <: UnsupervisedDetector
    store_precision::Bool = true
    assume_centered::Bool = false
    support_fraction::Union{Nothing, Real} = nothing
    random_state::Union{Nothing, Integer} = nothing
end

"""    PyOCSVM(kernel = "rbf",
               degree = 3,
               gamma = "auto",
               coef0 = 0.0,
               tol = 0.001,
               nu = 0.5,
               shrinking = true,
               cache_size = 200,
               verbose = false,
               max_iter = -1)
$(make_docs_link("ocsvm"))
"""
@pymodel mutable struct PyOCSVM <: UnsupervisedDetector
    kernel::String = "rbf"::(_ in ("linear", "poly", "rbf", "sigmoid", "precomputed"))
    degree::Integer = 3::(_ > 1)
    gamma::Union{String, Real} = "auto"
    coef0::Real = 0.0
    tol::Real = 0.001
    nu::Real = 0.5::(0 < _ <= 1)
    shrinking::Bool = true
    cache_size::Integer = 200::(_ > 0)
    verbose::Bool = false
    max_iter::Integer = -1
end

"""    PyPCA(n_components = nothing,
             n_selected_components = nothing,
             copy = true,
             whiten = false,
             svd_solver = "auto",
             tol = 0.0
             iterated_power = "auto",
             standardization = true,
             weighted = true,
             random_state = nothing)
$(make_docs_link("pca"))
"""
@pymodel mutable struct PyPCA <: UnsupervisedDetector
    n_components::Union{Nothing, Real} = nothing
    n_selected_components::Union{Nothing, Integer} = nothing
    copy::Bool = true
    whiten::Bool = false
    svd_solver::String = "auto"::(_ in ("auto", "full", "arpack", "randomized"))
    tol::Real = 0.0
    iterated_power::Union{String, Integer} = "auto"
    standardization::Bool = true
    weighted::Bool = true
    random_state::Union{Nothing, Integer} = nothing
end

"""    PyROD(parallel_execution = false)
$(make_docs_link("rod"))
"""
@pymodel mutable struct PyROD <: UnsupervisedDetector
    parallel_execution::Bool = false
end

"""    PySOD(n_neighbors = 5,
             ref_set = 10,
             alpha = 0.8)
$(make_docs_link("sod"))
"""
@pymodel mutable struct PySOD <: UnsupervisedDetector
    n_neighbors::Integer = 20::(_ > 0)
    ref_set::Integer = 10::(_ > 0)
    alpha::Real = 0.8::(0 < _ < 1)
end

"""    PySOS(perplexity = 4.5,
             metric = "minkowski",
             eps = 1e-5)
$(make_docs_link("sos"))
"""
@pymodel mutable struct PySOS <: UnsupervisedDetector
    perplexity::Real = 4.5::(_ > 0)
    metric::String = "minkowski"::(_ in ("cityblock", "cosine", "euclidean", "l1", "l2", "manhatten", "braycurtis", "canberra", "chebyshev", "correlation", "dice", "hamming", "jaccard", "kulsinski", "mahalanobis", "matching", "minkowski", "rogerstanimoto", "russellrao", "seuclidean", "sokalmichener", "sokalsneath", "sqeuclidean", "yule"))
    eps::Real = 1e-5::(_ > 0)
end


text:
# series of deconvolution tools

export selectFactorType
export approxDeconv, deconvSolveKey

## Initial version of selecting the dimension of a factor -- will be consolidated with existing infrastructure later

"""
    $SIGNATURES

Inverse solve of predicted noise value and returns tuple of (newly predicted, and known "measured"
noise) values.

Notes
- Only works for first value in measurement::Tuple at this stage.
- "measured" is used as starting point for the "predicted" values solve.
- Not all factor evaluation cases are support yet.
- NOTE only works on `.threadid()==1` at present, see #1094
- This function is still part of the initial implementation and needs a lot of generalization improvements.

DevNotes
- TODO Test for various cases with multiple variables.
- TODO make multithread-safe, and able, see #1094
- TODO Test for cases with `nullhypo`
- FIXME FactorMetadata object for all use-cases, not just empty object.
- TODO resolve #1096 (multihypo)
  - TODO Test cases for `multihypo`.
- TODO figure out if there is a way to consolidate with evalFactor and approxConv?
  - basically how to do deconv for just one sample with unique values (wrt TAF)
- TODO N should not be hardcoded to 100

Related

[`approxDeconv`](@ref), [`_solveCCWNumeric!`](@ref)
"""
function approxDeconv(fcto::DFGFactor,
                      ccw::CommonConvWrapper = _getCCW(fcto);
                      N::Int=100,
                      measurement::Tuple=sampleFactor(ccw, N),
                      retries::Int=3  )
  #
  # but what if this is a partial factor -- is that important for general cases in deconv?
  _setCCWDecisionDimsConv!(ccw)

  # FIXME This does not incorporate multihypo??
  varsyms = getVariableOrder(fcto)
  # vars = getPoints.(getBelief.(dfg, varsyms, solveKey) )
  fcttype = getFactorType(fcto)

  # get measurement dimension
  zDim = _getZDim(fcto)
  # TODO consider using ccw.cpt[thrid].res # likely needs resizing
  res_ = zeros(zDim)
  # TODO, consolidate fmd with getSample/sampleFactor and _buildLambda
  fctSmpls = deepcopy(measurement[1])
  fmd = _getFMdThread(ccw)

  # TODO assuming vector on only first container in measurement::Tuple
  makeTarget = (i) -> view(measurement[1], :, i)

  # NOTE
  # build a lambda that incorporates the multihypo selections
  # set these first
  # ccw.cpt[].activehypo / .p / .params  # params should already be set from construction
  certainidx, allelements, activehypo, mhidx = assembleHypothesesElements!(nothing, N, 0, length(varsyms))
  # only doing the current active hypo
  @assert activehypo[2][1] == 1 "deconv was expecting hypothesis nr == (1, 1:d)"

  lent = length(makeTarget(1))
  islen1 = size(makeTarget(1),1) == 1

  for idx in 1:N
    # towards each particle in their own thread (not 100% ready yet, factors should be separate memory)
    thrid = Threads.threadid()
    cpt_ = ccw.cpt[thrid]
    targeti_ = makeTarget(idx)

    # TODO must first resolve hypothesis selection before unrolling them -- deferred #1096
    cpt_.activehypo = activehypo[2][2]

    onehypo!, _ = _buildCalcFactorLambdaSample( ccw,
                                                idx,
                                                cpt_,
                                                targeti_,
                                                measurement  )
    #

    # lambda with which to find best measurement values
    hypoObj = (tgt) -> (targeti_.=tgt; onehypo!() )

    # find solution via SubArray view pointing to original memory location
    targeti_ .= _solveLambdaNumeric(fcttype, hypoObj, res_, measurement[1][:,idx], islen1)
  end

  # return (deconv-prediction-result, independent-measurement)
  return measurement[1], fctSmpls
end

"""
    $SIGNATURES

Generalized deconvolution to find the predicted measurement values of the factor `fctsym` in `dfg`.
Inverse solve of predicted noise value and returns tuple of (newly predicted, and known "measured" noise) values.

Notes
- Opposite operation contained in `approxConv`.
- For more notes see [`solveFactorMeasurements`](@ref).

Related

[`approxConv`](@ref), `deconvSolveKey`
"""
function approxDeconv(dfg::AbstractDFG,
                      fctsym::Symbol,
                      solveKey::Symbol=:default;
                      retries::Int=3  )
  #

  # which factor
  fct = getFactor(dfg, fctsym)
  N = size(getPoints(getBelief(dfg,getVariableOrder(fct)[1])), 2)
  pts = approxDeconv(fct, N=N, retries=retries )
  return pts
end

"""
    $SIGNATURES

Calculate the relative difference between two variables and across solveKeys.

Example

```julia
fg = generateCanonicalFG_Kaess()
solveTree!(fg, storeOld=true)
# calculate the relative motion induced by the solver from init to solve.
pts = deconvSolveKey(fg, :x1, :default, :x1, :graphinit)
```

Notes
- Can pass user `tfg::AbstractDFG` for better in-place operations.

DevNotes
- TODO use dfg, rather than building new tfg internally.

Related

[`approxDeconv`](@ref), [`mmd`](@ref)
"""
function deconvSolveKey(dfg::AbstractDFG,
                        refSym::Symbol,
                        refKey::Symbol,
                        tstSym::Symbol,
                        tstKey::Symbol;
                        tfg = initfg()  )
  #
  # create a new temporary factor graph for calculations

  # add the first "reference" variable
  Xref = getBelief(dfg, refSym, refKey)
  refSym_ = Symbol(refSym, "_ref")
  refVarType = getVariableType(dfg, refSym)
  if !exists(tfg, refSym_)
    addVariable!(tfg, refSym_, refVarType)
  end
  initManual!(tfg, refSym_, Xref)

  # add the second "test" variable
  tstVarType = getVariableType(dfg, tstSym)
  Xtst = getBelief(dfg, tstSym, tstKey)
  tstSym_ = Symbol(tstSym, "_tst")
  if !exists(tfg, tstSym_)
    addVariable!(tfg, tstSym_, tstVarType)
  end
  initManual!(tfg, tstSym_, Xtst)

  # add the new dummy factor with default manifold for computations
  fctType = selectFactorType(refVarType, tstVarType)
  nf = addFactor!(tfg, [refSym_; tstSym_], fctType())

  # TODO connect from dfg all other data that might form part of FactorMetadata in tfg
  pts = approxDeconv(tfg, nf.label)

  # assuming tfg was passed in by the user
  deleteFactor!(tfg, nf.label)

  # return result
  return pts, fctType
end

#


text:
function compute_grad_num(n::Int64,h::Float64,t0::Float64,tmax::Float64,elements::Array{Float64,2},tt2::Array{Float64,2},
   count2::Array{Int64,1},rstar::Float64,tobs_tot::Array{Float64,1},sobs_tot::Array{Float64,1},lgndof::Float64,V1exp2nuinv::Float64,
   indx::Array{Int64,1},nparam::Int64,nplanet::Int64,nvary::Int64)
# compute_grad_num(n,h,t0,tmax,elements,tt2,count2,rstar,tobs_tot,sobs_tot,lgndof,V1exp2nuinv,indx,nparam,nplanet,nvary)
# Computes gradient numerically.
# First, compute model in BigFloat precision:
hbig = big(h); t0big = big(t0); tmaxbig = big(tmax); elements_big = big.(elements); rstarbig = big(rstar)
tt2big = big.(tt2)
zbig = zero(BigFloat)
@time dq0big = ttv_elements!(n,t0big,hbig,tmaxbig,elements_big,tt2big,count2,zbig,0,0,rstarbig)
lgndof_big = big(lgndof); V1exp2nuinv_big = big(V1exp2nuinv); tobs_tot_big = big.(tobs_tot); sobs_tot_big = big.(sobs_tot)
chi_ttv = zero(BigFloat)
for j=1:ntrans
  chi_ttv += lnprob([tobs_tot_big[j],tt2big[iplanet[j],indx[j]],sobs_tot_big[j],lgndof_big,V1exp2nuinv_big])
end
# Add in Student's t normalization:
chi_ttv += ntrans*log_norm_students_t([lgndof_big,V1exp2nuinv_big])
gradf_num = zeros(BigFloat,nparam)
i = 0
# Now perturb each element and compute numerical derivative:
println("Computing numerical gradient")
dq = big(1e-15)
for ip=1:nplanet, ivary=1:nvary
  tstart = time()
  i +=1
  # If varying mass, this is at end of dtdq0 parameter list:
  if ivary == 1; iv = 7; else; iv= ivary-1; end
  # Perturb the parameter associated with ivary:
  lgndof_1 = big(0.0); lgndof_2 = big(0.0); V1exp2nuinv_1 = big(0.0); V1exp2nuinv_2 = big(0.0)
  if ivary < nparam-1
    # Compute finite difference derivative with respect to elements[ip+1,ivary]
    elements_big_1 = copy(elements_big); tt2big_1 = copy(tt2big)
    elements_big_1[ip+1,ivary] = elements_big[ip+1,ivary]*(1-dq)
    elements_big_2 = copy(elements_big); tt2big_2 = copy(tt2big)
    elements_big_2[ip+1,ivary] = elements_big[ip+1,ivary]*(1+dq)
    dq0big1 = ttv_elements!(n,t0big,hbig,tmaxbig,elements_big_1,tt2big_1,count2,zbig,0,0,rstarbig)
    dq0big2 = ttv_elements!(n,t0big,hbig,tmaxbig,elements_big_2,tt2big_2,count2,zbig,0,0,rstarbig)
  end
  # Gradient of "chi-square" with respect to each varying parameter:
  for it=1:ntrans
    # Set up parameter vector for computing likelihood of this data point:
    if ivary < nparam-1
      p1 = [tobs_tot_big[it],tt2big_1[iplanet[it],indx[it]],sobs_tot_big[it],lgndof_big,V1exp2nuinv_big]
      p2 = [tobs_tot_big[it],tt2big_2[iplanet[it],indx[it]],sobs_tot_big[it],lgndof_big,V1exp2nuinv_big]
      gradf_num[i] += (lnprob(p2)-lnprob(p1))/(2dq*elements_big[ip+1,ivary])
    end
    # gradf_num[i] += grad_tmp[2]*dtdelements[iplanet[it],indx[it],iv,ip+1]
    # Compute with finite difference
    if use_student
      if i == 1
      # Derivative wrt ndof:
        # gradf_num[nparam-1] += grad_tmp[4]
        lgndof_1 = lgndof_big *(1-dq)
        lgndof_2 = lgndof_big *(1+dq)
        p1 = [tobs_tot_big[it],tt2big[iplanet[it],indx[it]],sobs_tot_big[it],lgndof_1,V1exp2nuinv_big]
        p2 = [tobs_tot_big[it],tt2big[iplanet[it],indx[it]],sobs_tot_big[it],lgndof_2,V1exp2nuinv_big]
        deriv_num = (lnprob(p2)-lnprob(p1))/(2dq*lgndof_big)
        if isnan(deriv_num)
          println("NaN derivative.  p1: ",p1," p2: ",p2," dq: ",dq," lgndof_big: ",lgndof_big," lnprob(p1): ",lnprob(p1)," lnprob(p2): ",lnprob(p2))
        else
          gradf_num[nparam-1] += deriv_num
        end
        # Derivative wrt V1exp2nuinv:
        # gradf_num[nparam  ] += grad_tmp[5]
        V1exp2nuinv_1 = V1exp2nuinv_big * (1-dq)
        V1exp2nuinv_2 = V1exp2nuinv_big * (1+dq)
        p1 = [tobs_tot_big[it],tt2big[iplanet[it],indx[it]],sobs_tot_big[it],lgndof_big,V1exp2nuinv_1]
        p2 = [tobs_tot_big[it],tt2big[iplanet[it],indx[it]],sobs_tot_big[it],lgndof_big,V1exp2nuinv_2]
#        gradf_num[nparam] +=
        deriv_num = (lnprob(p2)-lnprob(p1))/(2dq*V1exp2nuinv_big)
        if isnan(deriv_num)
          println("NaN derivative.  p1: ",p1," p2: ",p2," dq: ",dq," V1exp2nuinv_big: ",V1exp2nuinv_big," lnprob(p1): ",lnprob(p1)," lnprob(p2): ",lnprob(p2))
        else
          gradf_num[nparam] += deriv_num
        end
      end
    end
  end
  telapse = time()- tstart
  println("Done with component: ",ip,", ",ivary," grad: ",gradf_num[i]," time: ",telapse/60," min" )
end
# Now add in derivative of normalization:
lgndof_1 = lgndof_big *(1-dq)
lgndof_2 = lgndof_big *(1+dq)
norm_1 = log_norm_students_t([lgndof_1,V1exp2nuinv_big])
norm_2 = log_norm_students_t([lgndof_2,V1exp2nuinv_big])
deriv_num = ntrans*(norm_2-norm_1)/(2dq*lgndof_big)
#gradf_num[nparam-1] += ntrans*(norm_2-norm_1)/(2*lgndof_big*dq)
if isnan(deriv_num)
  println("NaN derivative.  lgndof_1: ",lgndof_1," lgndof_2: ",lgndof_2," dq: ",dq," V1exp2nuinv_big: ",V1exp2nuinv_big," norm_1: ",norm_1," norm_2: ",norm_2," ntrans: ",ntrans)
else
  gradf_num[nparam-1] += deriv_num
end
V1exp2nuinv_1 = V1exp2nuinv_big * (1-dq)
V1exp2nuinv_2 = V1exp2nuinv_big * (1+dq)
norm_1 = log_norm_students_t([lgndof_big,V1exp2nuinv_1])
norm_2 = log_norm_students_t([lgndof_big,V1exp2nuinv_2])
deriv_num = ntrans*(norm_2-norm_1)/(2dq*V1exp2nuinv_big)
#gradf_num[nparam  ] += ntrans*(norm_2-norm_1)/(2*V1exp2nuinv_big*dq)
if isnan(deriv_num)
  println("NaN derivative.  V1exp2nuinv_1: ",V1exp2nuinv_1," V1exp2nuinv_2: ",V1exp2nuinv_2," dq: ",dq," lgndof_big: ",lgndof_big," norm_1: ",norm_1," norm_2: ",norm_2," ntrans: ",ntrans)
else
  gradf_num[nparam] += deriv_num
end
return gradf_num
end


text:
export Operator

"""
    abstract type Operator{T<:Number} <: AbstractMatrix{Number} end

A matrix representing a linear operator that acts upon a complex-valued Hilbert space.
The `Operator` is an abstract type that parents all linear operator in quantum mechanics.
These matrix types can represent quantum states, evolution, and measurements, each with
their individual constraint.
These constraints are place upon the children of th `Operator` abstract type.
"""
abstract type Operator{T<:Number} <: AbstractMatrix{Number} end
Base.size(O :: Operator) = size(O.M)
Base.getindex(O :: Operator, id :: Vararg{Int,2}) = getindex(O.M, id...)
Base.setindex!(O :: Operator, val, id :: Vararg{Int,2}) = (O.M[id...] = val)

"""
Operator Multiplication:

```julia
*(operators :: Vararg{Operator}) :: Matrix
```
"""
*(operators :: Vararg{Operator}) :: Matrix = *(map(O -> O.M, operators)...)

"""
Kronecker Product:

```julia
kron(operators :: Vararg{Operator}) :: Matrix
```
"""
kron(operators :: Vararg{Operator}) :: Matrix =  kron(map(O -> O.M, operators)...)

"""
Matrix Rank:

```julia
rank(O :: Operator) :: Int64
```
"""
rank(O :: Operator) :: Int64 = rank(O.M, atol=O.atol)

"""
Matrix Square Root:

```julia
sqrt(O :: Operator) :: Matrix
```
"""
sqrt(O :: Operator) :: Matrix = sqrt(O.M)

"""
[`Operator`](@ref) types can multiply [`Bra`](@ref) and [`Ket`](@ref) types.

``O|\\psi\\rangle = | \\psi' \\rangle``:

```julia
*(O :: Operator, ket :: Ket) :: Vector
```

``\\langle \\psi |O = \\langle \\psi'|``:

```julia
*(bra :: Bra, O :: Operator) :: Adjoint{T, Vector{T}} where T <: Number
```

Inner product, ``\\langle \\psi |O|\\sigma\\rangle = \\langle \\psi|\\sigma'\\rangle``:

```julia
*(bra :: Bra, O :: Operator, ket :: Ket) :: Number
```
"""
*(O :: Operator, ket :: Ket) :: Vector = O.M * ket.ψ
*(bra :: Bra, O :: Operator) :: Adjoint{T, Vector{T}} where T <: Number = bra.ψ * O.M
*(bra :: Bra, O :: Operator, ket :: Ket) :: Number = bra.ψ * O.M * ket.ψ

*(O :: Operator, M :: Matrix{<:Number}) :: Matrix = O.M * M
*(M :: Matrix{<:Number}, O :: Operator) :: Matrix = M * O.M


text:
## ----------------------------------------------------------------------------

const _KBASE_CONFIG_KEY = "knowledgebase"

function _knowledgebase()::String
    config = read_configfile()
    kb = get(config, _KBASE_CONFIG_KEY, "")
    return kb
end


text:
a = 1
b = 2
println("a = $(a), b = $(b)")
a ⊻= b
b ⊻= a
a ⊻= b
println("a = $(a), b = $(b)")


text:
import YAML
import Random
using ArgParse
using Bayesian_Survival_Analysis
using CSV
using DataFrames
# using Revise

load_config(path::String) = YAML.load_file(path; dicttype=Dict{Symbol, Any})
set_seed(conf::Dict) = Random.seed!(conf[:seed])

function generate_dummy_df(conf::Dict)
  set_seed(conf)

  df = DataFrame(
    # 入社した時の年齢
    old = cat([generateold(;conf[:generateold][i]...) for i in 1:3]..., dims=1),
    # 性別
    gender = cat([generategender(;conf[:generategender][i]...) for i in 1:3]..., dims=1),
    # 役職(1: 新卒, 2: 中途, 3: 役員)
    workertype = cat(ones(Int, conf[:proper_num]), ones(Int, conf[:mid_career_num]) * 2, ones(Int, conf[:high_career_num]) * 3, dims=1)
  )
  # 学位(0: 高卒, 1: 学士, 2: 修士, 3: 博士)
  insertcols!(df, 1, :gakui => generategakui(;df=df, stats=conf[:generategakui][:stats]))

  # 入社年(何年前に入社したか?)
  insertcols!(df, 1, :yearofjoin => generateyearofjoin(;df=df, stats=conf[:generateyearofjoin][:stats]))

  # 勤続年数
  insertcols!(df, 1, :seniority => generateseniority(;df=df))

  # 打ち切りフラグ
  insertcols!(df, 1, :censored => generatecensored(;df=df))
end

output_csv(df::DataFrame, conf::Dict) = CSV.write(conf[:output_path], df, delim=',', writeheader=true)

function parse_commandline()
  s = ArgParseSettings()

  @add_arg_table s begin
    "--conf", "-c"
      help = "config path"
      arg_type = String
      default = "./configs/generate_data_config.yaml"
  end

  return parse_args(s; as_symbols=true)
end

function main()
  parsed_args = parse_commandline()
  @show parsed_args
  conf = load_config(parsed_args[:conf])
  df = generate_dummy_df(conf)
  show(df)
  output_csv(df, conf)
end

main()