text:
using Combinatorics: combinations
using StatsBase: sample

@doc raw"""
    get_contours(fit::AbstractFit, χsq, parameters_combination::Vector{Int};
                 npts::Int=20, limits=true, sigma = 1.0)

For a given fit `fit` and the ``\chi^2`` function `χsq`, gives an array of parameter arrays,
with each array corresponding to a set of parameters obtained from calculating the `MINOS` ``1σ`` contour
(try to find `npts` points in the contour) for the two parameters in `parameters_combination`.

`parameters_combination` is an `Int` array of the numbers of that two parameters, e.g. it is `[1, 2]` for the first
two parameters and `[2, 3]` or the second and third parameters.

If `limits` is `true`, then fix one parameter to its bounds from `MINOS` of the best fit and get the values for
the other parameters; this runs over all parameters.
"""
function get_contours(fit::Fit, χsq, parameters_combination::Vector{Int}; npts::Int=20, limits=true, sigma = 1.0)
    fit.migrad_ok() || migrad(fit);        # if migrad has not been run then run it first
    length(fit.merrors) == 0 && minos(fit) # if minos has not been run then run it first
    fmin_1σ = fit.fval + 1.0;   # χsq from the previous best fit + 1.0
    tol = 0.05  # tolerance allowing χ² to be ≤  fmin_1σ + tol
    # setting the parameters to the the best ones from the previous fit
    kwdarg = Dict{Symbol, Any}(Symbol(k) => v for (k,v) in fit.fitarg)
    _args = Symbol.(fit.parameters) #  func_argnames(χsq)
    @views for i in 1: length(_args)    # reset the limits of parameters to be the bounds from minos
        if !kwdarg[Symbol(:fix_, _args[i])]
            kwdarg[Symbol(:limit_, _args[i])] = get(fit.args, i-1) .+
            (fit.merrors[String(_args[i])][3], fit.merrors[String(_args[i])][4])
             # (fit.merrors[String(_args[i]), -1], fit.merrors[String(_args[i]), 1])
        end
    end

    container = Vector{Vector{Float64}}()
    # consider the upper and lower bounds for each parameter
    if limits
        @views for a in _args
            if !kwdarg[Symbol(:fix_, a)]
                _fit1 = Minuit(χsq; kwdarg..., a => kwdarg[Symbol(:limit_, a)][1], Symbol("fix_", a) => true)
                _fit2 = Minuit(χsq; kwdarg..., a => kwdarg[Symbol(:limit_, a)][2], Symbol("fix_", a) => true)
                _fit1.strategy = 1; migrad(_fit1)
                _fit2.strategy = 1; migrad(_fit2)
                _fit1.fval  < fmin_1σ + tol && push!(container, vcat(_fit1.fval, args(_fit1)) )
                _fit2.fval  < fmin_1σ + tol && push!(container, vcat(_fit2.fval, args(_fit2)) )
            end
        end
    end

    # choose a pair of parameters, and try to compute the MINOS contour of npts points
    para1, para2 = _args[parameters_combination[1]], _args[parameters_combination[2]]
    contour_parameters = fit.mncontour(para1, para2, sigma = sigma, numpoints=npts )[3]

    # for each contour point, get the values of the other parameters
    for pa in contour_parameters
        dict = Dict(zip((para1, para2), pa) ) # using dict to overwrite those in kwdarg
        _fit1 = Minuit(χsq; kwdarg..., dict...,  #errordef=1,
                  Symbol("fix_", para1) => true, Symbol("fix_", para2) => true)
        _fit1.strategy = 1; migrad(_fit1)
        # filtering only those with χ² in 1σ
        _fit1.fval  < fmin_1σ + tol && push!(container, vcat(_fit1.fval, args(_fit1)) )
    end

    return container
end

function get_contours(fit::ArrayFit, χsq, parameters_combination::Vector{Int}; npts::Int=20, limits=true, sigma = 1.0)
    fit.migrad_ok() || migrad(fit);        # if migrad has not been run then run it first
    length(fit.merrors) == 0 && minos(fit) # if minos has not been run then run it first
    fmin_1σ = fit.fval + 1.0;   # χsq from the previous best fit + 1.0
    tol = 0.05  # tolerance allowing χ² to be ≤  fmin_1σ + tol
    # setting the parameters to the the best ones from the previous fit
    start_array = args(fit) #PyVector{Real}(fit.args) # does not support assignment
    kwdarg = Dict{Symbol, Any}(Symbol(k) => v for (k,v) in fit.fitarg)
    _args = Symbol.(fit.parameters) # get names of parameters
    nargs = length(_args)
    @views for i in 1: nargs    # reset the limits of parameters to be the bounds from minos
        if !kwdarg[Symbol(:fix_, _args[i])]
            kwdarg[Symbol(:limit_, _args[i])] = get(fit.args, i-1) .+
            (fit.merrors[String(_args[i])][3], fit.merrors[String(_args[i])][4])
        end
    end

    # igrad ? gradfun(par) = gradient(χsq, par) : gradfun = nothing

    container = Vector{Vector{Float64}}()
    # consider the upper and lower bounds for each parameter
    if limits
        @views for i in 1:nargs # for a in _args
            if !kwdarg[Symbol(:fix_, _args[i])]
                a = _args[i]
                start_array[i] = kwdarg[Symbol(:limit_, a)][1]
                _fit1 = Minuit(χsq, start_array; name = _args, kwdarg..., Symbol("fix_", a) => true)
                start_array[i] = kwdarg[Symbol(:limit_, a)][2]
                _fit2 = Minuit(χsq, start_array; name = _args, kwdarg..., Symbol("fix_", a) => true)
                _fit1.strategy = 1; migrad(_fit1)
                _fit2.strategy = 1; migrad(_fit2)
                _fit1.fval  < fmin_1σ + tol && push!(container, vcat(_fit1.fval, args(_fit1)) )
                _fit2.fval  < fmin_1σ + tol && push!(container, vcat(_fit2.fval, args(_fit2)) )
            end
        end
    end

    # choose a pair of parameters, and try to compute the MINOS contour of npts points
    para1, para2 = _args[parameters_combination[1]], _args[parameters_combination[2]]
    contour_parameters = fit.mncontour(para1, para2, sigma = sigma, numpoints=npts )[3]

    # for each contour point, get the values of the other parameters
    for pa in contour_parameters
        start_array[parameters_combination[1]], start_array[parameters_combination[2]] = pa
        # start_array[parameters_combination[2]] = pa[2]
        _fit1 = Minuit(χsq, start_array; name = _args, kwdarg..., # dict...,  #errordef=1,
                  Symbol("fix_", para1) => true, Symbol("fix_", para2) => true)
        _fit1.strategy = 1; migrad(_fit1)
        # filtering only those with χ² in 1σ
        _fit1.fval  < fmin_1σ + tol && push!(container, vcat(_fit1.fval, args(_fit1)) )
    end

    return container
end

@doc raw"""
    get_contours_all(fit::AbstractFit, χsq; npts=20, limits=true, sigma = 1.0)

For a given fit `fit` and the ``\chi^2`` function `χsq`, gives a list of parameters sets which are at the edge of
``1σ`` `MINOS` contours for all combinations of varying parameters. The case of `limits` being `true` runs only once.
"""
function get_contours_all(fit::AbstractFit, χsq; npts=20, limits=true, sigma = 1.0)
    npara = length(func_argnames(χsq))
    container = Vector()
    push!(container, vcat(fit.fval, args(fit)) ) # the first row set to tbe best fit
    # npara = fit.narg
    free_pars_indices = findall(x-> x in fit.list_of_vary_param(), fit.parameters)
    for arr in combinations(free_pars_indices, 2)  # (1:npara, 2)
        push!(container, get_contours(fit, χsq, arr, npts=npts, limits=limits, sigma = sigma)...)
        limits = false  # run limits only once
    end
    return container
end

"""
    contour_df(fit::AbstractFit, χsq; npts=20, limits=true, sigma = 1.0)

parameters in the form of a dataframe.
"""
function contour_df(fit::AbstractFit, χsq; npts=20, limits=true, sigma = 1.0)
    parameters_1sigma =  get_contours_all(fit, χsq, npts=npts, limits=limits, sigma = sigma)
    df0 = DataFrame(parameters_1sigma)
    argnames = Array{Symbol,1}(undef, fit.narg)
    @. argnames = Symbol(fit.parameters)
    try
        DataFrame( collect.(eachrow(df0)), vcat(:chisq, argnames) )
    catch
        @warn "No parameter sets were found."
    end
end

"""
    get_contours_given_parameter(fit::AbstractFit, χsq, para::T, range)
        where {T <: Union{Symbol, String}}

gives parameter sets in one sigma for a given parameter constrained in a range.
If `fit` is an `ArrayFit` and no user-defined names have been given to the parameters,
then `para` is `"x0"` or `:x0` for the 1st parameter, `"x1"` or `:x1` for the 2nd parameter, ...
"""
function get_contours_given_parameter(fit::Fit, χsq, para::T, range) where {T <: Union{Symbol, String}}
    fit.migrad_ok() || migrad(fit);        # if migrad has not been run then run it first
    length(fit.merrors) == 0 && minos(fit) # if minos has not been run then run it first
    fmin_1σ = fit.fval + 1.0;   # χsq from the previous best fit + 1.0
    # setting the parameters to the the best ones from the previous fit
    kwdarg = Dict{Symbol, Any}(Symbol(k) => v for (k,v) in fit.fitarg)
    _args = Symbol.(fit.parameters) #  func_argnames(χsq)
    @views for i in 1: length(_args)    # reset the limits of parameters to be the bounds from minos
        if !kwdarg[Symbol(:fix_, _args[i])]
            kwdarg[Symbol(:limit_, _args[i])] = get(fit.args, i-1) .+
             (fit.merrors[String(_args[i])][3], fit.merrors[String(_args[i])][4])
              # (fit.merrors[String(_args[i]), -1], fit.merrors[String(_args[i]), 1])
        end
    end

    container = Vector{Vector{Float64}}()
    for a in range
        _fit1 = Minuit(χsq; kwdarg..., para => a, Symbol(:fix_, para) => true)
        _fit1.strategy = 2; _fit1.migrad()
        _fit1.fval  ≤ fmin_1σ && push!(container, vcat(_fit1.fval, args(_fit1)) )
    end

    return container
end

function get_contours_given_parameter(fit::ArrayFit, χsq, para::T, range) where {T <: Union{Symbol, String}}
    fit.migrad_ok() || migrad(fit);        # if migrad has not been run then run it first
    length(fit.merrors) == 0 && minos(fit) # if minos has not been run then run it first
    fmin_1σ = fit.fval + 1.0;   # χsq from the previous best fit + 1.0
    # setting the parameters to the the best ones from the previous fit
    start_array = args(fit) #PyVector{Real}(fit.args) # does not support assignment
    kwdarg = Dict{Symbol, Any}(Symbol(k) => v for (k,v) in fit.fitarg)
    _args = Symbol.(fit.parameters)
    nargs = length(_args)
    @views for i in 1:nargs    # reset the limits of parameters to be the bounds from minos
        if !kwdarg[Symbol(:fix_, _args[i])]
            kwdarg[Symbol(:limit_, _args[i])] = get(fit.args, i-1) .+
             (fit.merrors[String(_args[i])][3], fit.merrors[String(_args[i])][4])
        end
    end

    container = Vector{Vector{Float64}}()
    for a in range
        start_array[findfirst(x-> x == Symbol(para), _args)] = a
        _fit1 = Minuit(χsq, start_array; name = _args, kwdarg..., Symbol(:fix_, para) => true)
        _fit1.strategy = 2; migrad(_fit1)
        _fit1.fval  ≤ fmin_1σ && push!(container, vcat(_fit1.fval, args(_fit1)) )
    end

    return container
end

@doc """
    contour_df_given_parameter(fit::AbstractFit, χsq, para::T, range; limits = true)
        where {T <: Union{Symbol, String}}

return parameter sets in one sigma for a given parameter constrained in a range as a `DataFrame`.

See also [`get_contours_given_parameter`](@ref).
"""
function contour_df_given_parameter(fit::AbstractFit, χsq, para::T, range; limits = true) where {T <: Union{Symbol, String}}
    parameters =  get_contours_given_parameter(fit, χsq, para, range)
    df0 = DataFrame(parameters)
    argnames = Array{Symbol,1}(undef, fit.narg)
    @. argnames = Symbol(fit.parameters)
    try
        DataFrame( collect.(eachrow(df0)), vcat(:chisq, argnames) )
    catch
        @warn "No parameter sets were found in the given range of $para. Try another range."
    end
end

@doc raw"""
    get_contours_samples(fit::AbstractFit, χsq, paras, ranges; nsamples = 100, MNbounds = true)

return 1σ parameter sets as an `Array` (the latter returns a `DataFrame`) for given parameters constrained in `ranges`:
* if `paras` is a single parameter, then take equally spaced `nsamples` in `ranges` given in the form of `(min, max)`;
* if `paras` contain more (``\geq 2``) parameters, then `paras` should be of the form `(:para1, :para2)`, `ranges` should be of the form `((min1, max1), (min2, max2))`,
and values for the parameters given in `paras` are randomly sampled in the given `ranges`;
* if `MNbounds` is true, then constrain the parameters in the range provided by `MINOS` no matter whether that is valid or not (to be improved by checking the validity)
* if `igrad` is true, then use `ForwardDiff.gradient` to compute the gradient.
* For using array parameters, if no user-defined names have been given to the parameters,
`paras` should be given such that `"x0"` or `:x0` for the 1st parameter, `"x1"` or `:x1` for the 2nd parameter, ...
"""
function get_contours_samples(fit::Fit, χsq, paras, ranges; nsamples = 100, MNbounds = true, igrad = false)
    fit.migrad_ok() || migrad(fit);        # if migrad has not been run then run it first
    length(fit.merrors) == 0 && minos(fit) # if minos has not been run then run it first
    fmin_1σ = fit.fval + 1.0;   # χsq from the previous best fit + 1.0
    # setting the parameters to the the best ones from the previous fit
    kwdarg = Dict{Symbol, Any}(Symbol(k) => v for (k,v) in fit.fitarg)
    _args = Symbol.(fit.parameters)
    if MNbounds
        @views for i in 1: length(_args)    # reset the limits of parameters to be the bounds from minos
            if !kwdarg[Symbol(:fix_, _args[i])]
            kwdarg[Symbol(:limit_, _args[i])] = get(fit.args, i-1) .+
                 (fit.merrors[String(_args[i])][3], fit.merrors[String(_args[i])][4])
                  # (fit.merrors[String(_args[i]), -1], fit.merrors[String(_args[i]), 1])
            end
        end
    end

    container = Vector{Vector{Float64}}()
    push!(container, vcat(fit.fval, args(fit)) ) # the first row set to tbe best fit

    # make range for each parameter with length len; if one only 1 parameter
    @views if isa(ranges[1], Number)
        sam = LinRange(ranges[1], ranges[end], nsamples)
        pfix = Symbol(:fix_, paras) => true
    else
        r = map(x -> LinRange(x[1], x[end], nsamples), ranges)
        sam = zip(sample.(r, nsamples)...)
        pfix = map( x -> Symbol(:fix_, x) => true, paras)
    end

    igrad ? gradf(a...) = gradient(x->χsq(x...), [a...]) : gradf = nothing

    @inbounds @simd for p in unique(sam)
        @views if isa(ranges[1], Number)
            _fit1 = Minuit(χsq; kwdarg..., paras => p, Symbol(:fix_, paras) => true, grad = gradf)
        else
            pdict = Dict(zip(paras, p))
            _fit1 = Minuit(χsq; kwdarg..., pdict..., pfix..., grad = gradf)
        end
        _fit1.strategy = 1; migrad(_fit1)
        _fit1.fval  ≤ fmin_1σ && push!(container, vcat(_fit1.fval, args(_fit1)) )
    end

    return container
end

function get_contours_samples(fit::ArrayFit, χsq, paras, ranges; nsamples = 100, MNbounds = true, igrad = false)
    fit.migrad_ok() || migrad(fit);        # if migrad has not been run then run it first
    length(fit.merrors) == 0 && minos(fit) # if minos has not been run then run it first
    fmin_1σ = fit.fval + 1.0;   # χsq from the previous best fit + 1.0
    # setting the parameters to the the best ones from the previous fit
    start_array = args(fit)
    kwdarg = Dict{Symbol, Any}(Symbol(k) => v for (k,v) in fit.fitarg)
    _args = Symbol.(fit.parameters)
    if MNbounds
        @views for i in 1: length(_args)    # reset the limits of parameters to be the bounds from minos
            if !kwdarg[Symbol(:fix_, _args[i])]
                kwdarg[Symbol(:limit_, _args[i])] = get(fit.args, i-1) .+
                 (fit.merrors[String(_args[i])][3], fit.merrors[String(_args[i])][4])
            end
        end
    end

    container = Vector{Vector{Float64}}()
    push!(container, vcat(fit.fval, args(fit)) ) # the first row set to tbe best fit

    # make range for each parameter with length len; if one only 1 parameter
    @views if isa(ranges[1], Number)
        sam = LinRange(ranges[1], ranges[2], nsamples)
        pfix = Symbol(:fix_, paras) => true
    else
        r = map(x -> LinRange(x[1], x[2], nsamples), ranges)
        sam = zip(sample.(r, nsamples)...)
        pfix = map( x -> Symbol(:fix_, x) => true, paras)
    end

    igrad ? gradf(a...) = gradient(x->χsq(x...), [a...]) : gradf = nothing

# if using @inbounds, then crashes julia
    for p in unique(sam)
        @views if isa(ranges[1], Number)
            start_array[findfirst(x-> x == Symbol(paras), _args)] = p
            _fit1 = Minuit(χsq, start_array; name = _args, kwdarg..., Symbol(:fix_, paras) => true, grad = gradf)
        else
            @views for i in eachindex(paras)
                start_array[findfirst(x-> x == Symbol(paras[i]), _args)] = p[i]
            end
            _fit1 = Minuit(χsq, start_array; name = _args, kwdarg..., pfix..., grad = gradf)
        end
        _fit1.strategy = 1; migrad(_fit1)
        _fit1.fval  ≤ fmin_1σ && push!(container, vcat(_fit1.fval, args(_fit1)) )
    end

    return container
end

@doc raw"""
    contour_df_samples(fit::AbstractFit, χsq, paras, ranges; nsamples = 100, MNbounds=true)

gives  1σ parameter sets as a `DataFrame` for given parameters constrained in `ranges`:
* if `paras` is a single parameter, then take equally spaced `nsamples` in `ranges` given in the form of `(min, max)`;
* if `paras` contain more parameters, then `paras` should be of the form `(:para1, :para2)`, `ranges` should be of the form `((min1, max1), (min2, max2))`;
* `paras` can be more than 2. Values for the parameters given in `paras` are randomly sampled in the given `ranges`.
* is `MNbounds` is true, then constrain the parameters in the range provided by `MINOS` no matter whether that is valid or not (to be improved by checking the validity)
* if `igrad` is true, then use `ForwardDiff.gradient` to compute the gradient
"""
function contour_df_samples(fit::AbstractFit, χsq, paras, ranges; nsamples = 100, MNbounds=true, igrad = false)
    parameters = get_contours_samples(fit, χsq, paras, ranges, nsamples = nsamples, MNbounds = MNbounds, igrad = igrad)
    df0 = DataFrame(parameters)
    argnames = Array{Symbol,1}(undef, fit.narg)
    @. argnames = Symbol(fit.parameters)
    try
        DataFrame( collect.(eachrow(df0)), vcat(:chisq, argnames) )
    catch
        @warn "No parameter sets were found in the given ranges of $para. Try another ranges."
    end
end


text:
@testset "problems" begin
    is_ci = get(ENV, "CI", nothing) == "true"
    if is_ci
        test_files = filter(isfile, joinpath.(@__DIR__, readdir(@__DIR__)))
    else
        repo = LibGit2.GitRepo(root)
        problem_dir = relpath(@__DIR__, root) # "test/problems"

        test_files = filter(diff_files(repo, "HEAD~1", "")) do f
            startswith(f, problem_dir) && endswith(f, ".jl")
        end
    end
    # don't recursively include "problem.jl" itself
    test_files = filter(x -> x != "problems.jl", basename.(test_files))

    @info "testing $(length(test_files)) problems" CI = is_ci
    for f in test_files
        # run tests in a sandbox
        m = Module(gensym())
        # modules created with Module() does not have include defined
        Core.eval(m, :(include(x) = Base.include($m, x)))
        Core.eval(m, :(include("../requirements.jl")))
        Core.eval(m, :(include($f)))
    end
end

ext: jl


text:
# Automatically generated using Clang.jl

const XKB_KEYCODE_INVALID = Float32(0x0fffffff)
const XKB_LAYOUT_INVALID = Float32(0x0fffffff)
const XKB_LEVEL_INVALID = Float32(0x0fffffff)
const XKB_MOD_INVALID = Float32(0x0fffffff)
const XKB_LED_INVALID = Float32(0x0fffffff)
const XKB_KEYCODE_MAX = Float32(0x0fffffff) - 1

# Skipping MacroDefinition: xkb_keycode_is_legal_ext ( key ) ( key <= XKB_KEYCODE_MAX )
# Skipping MacroDefinition: xkb_keycode_is_legal_x11 ( key ) ( key >= 8 && key <= 255 )
# Skipping MacroDefinition: XKB_KEYMAP_USE_ORIGINAL_FORMAT ( ( enum xkb_keymap_format ) - 1 )

const xkb_context = Cvoid
const xkb_keymap = Cvoid
const xkb_state = Cvoid
const xkb_keycode_t = UInt32
const xkb_keysym_t = UInt32
const xkb_layout_index_t = UInt32
const xkb_layout_mask_t = UInt32
const xkb_level_index_t = UInt32
const xkb_mod_index_t = UInt32
const xkb_mod_mask_t = UInt32
const xkb_led_index_t = UInt32
const xkb_led_mask_t = UInt32

struct xkb_rule_names
    rules::Cstring
    model::Cstring
    layout::Cstring
    variant::Cstring
    options::Cstring
end

@cenum xkb_keysym_flags::UInt32 begin
    XKB_KEYSYM_NO_FLAGS = 0
    XKB_KEYSYM_CASE_INSENSITIVE = 1
end

@cenum xkb_context_flags::UInt32 begin
    XKB_CONTEXT_NO_FLAGS = 0
    XKB_CONTEXT_NO_DEFAULT_INCLUDES = 1
    XKB_CONTEXT_NO_ENVIRONMENT_NAMES = 2
end

@cenum xkb_log_level::UInt32 begin
    XKB_LOG_LEVEL_CRITICAL = 10
    XKB_LOG_LEVEL_ERROR = 20
    XKB_LOG_LEVEL_WARNING = 30
    XKB_LOG_LEVEL_INFO = 40
    XKB_LOG_LEVEL_DEBUG = 50
end

@cenum xkb_keymap_compile_flags::UInt32 begin
    XKB_KEYMAP_COMPILE_NO_FLAGS = 0
end

@cenum xkb_keymap_format::UInt32 begin
    XKB_KEYMAP_FORMAT_TEXT_V1 = 1
end

const xkb_keymap_key_iter_t = Ptr{Cvoid}

@cenum xkb_key_direction::UInt32 begin
    XKB_KEY_UP = 0
    XKB_KEY_DOWN = 1
end

@cenum xkb_state_component::UInt32 begin
    XKB_STATE_MODS_DEPRESSED = 1
    XKB_STATE_MODS_LATCHED = 2
    XKB_STATE_MODS_LOCKED = 4
    XKB_STATE_MODS_EFFECTIVE = 8
    XKB_STATE_LAYOUT_DEPRESSED = 16
    XKB_STATE_LAYOUT_LATCHED = 32
    XKB_STATE_LAYOUT_LOCKED = 64
    XKB_STATE_LAYOUT_EFFECTIVE = 128
    XKB_STATE_LEDS = 256
end

@cenum xkb_state_match::UInt32 begin
    XKB_STATE_MATCH_ANY = 1
    XKB_STATE_MATCH_ALL = 2
    XKB_STATE_MATCH_NON_EXCLUSIVE = 65536
end

@cenum xkb_consumed_mode::UInt32 begin
    XKB_CONSUMED_MODE_XKB = 0
    XKB_CONSUMED_MODE_GTK = 1
end

const XKB_X11_MIN_MAJOR_XKB_VERSION = 1
const XKB_X11_MIN_MINOR_XKB_VERSION = 0

@cenum xkb_x11_setup_xkb_extension_flags::UInt32 begin
    XKB_X11_SETUP_XKB_EXTENSION_NO_FLAGS = 0
end


text:
"""
Feature extractor.

Performs Shi-Tomasi feature extraction algorithm with support
for masking-out regions in which to avoid feature detection.
"""
struct Extractor
    """
    Maximum number of points to detect.
    """
    max_points::Int
    descriptor::BRIEF
    """
    Radius of the circle that is used to create a mask with regions
    in which to avoid feature detection.
    """
    radius::Int
    grid_resolution::Tuple{Int, Int}
    cell_size::Int
end
Extractor(max_points, radius, grid_resolution, cell_size) =
    Extractor(max_points, BRIEF(;size=256), radius, grid_resolution, cell_size)

function _shi_tomasi(image, n_keypoints::Int; min_response::Real = 1e-4)
    corners = fill(false, size(image)...)
    responses = shi_tomasi(image)
    maxima = findlocalmaxima(responses)

    max_responses = [responses[mx] for mx in maxima]
    P = sortperm(max_responses; lt=(x, y) -> x > y)
    length(P) > n_keypoints && (P = P[1:n_keypoints];)

    nb_detected = 0
    for p in P
        maxima_corner = maxima[p]
        responses[maxima_corner] < min_response && continue
        corners[maxima_corner] = true
        nb_detected += 1
    end

    corners, responses, nb_detected
end

"""
Detect keypoints in the `image`.

# Arguments:

- `current_points`:
    Vector of points, which define circular regions, where to avoid
    detecting new features. This can be used to decrease the amount of
    similar keypoints. Pass empty vector to skip this step.
- `σ::Real`:
    Standard deviation for the gaussian blur.
    This is used to smooth out circles in the mask to reduce the chance
    of detecting aliased circle corners. Default value is `3`.
    Pass `0` to skip this step.

# Returns:

Keypoints in the `(y, x)` format.
"""
function detect(e::Extractor, image, current_points; σ_mask = 3)
    length(current_points) ≥ e.max_points && return CartesianIndex{2}[]

    if !isempty(current_points)
        mask = get_mask(image, current_points, e.radius)
        # Smooth out mask, to avoid detecting features on the circle edges.
        σ_mask ≉ 0 && (mask = imfilter(mask, Kernel.gaussian(σ_mask));)
        image = image .* mask
    end

    height, width = size(image)
    n_cells = e.grid_resolution[1] * e.grid_resolution[2]
    n_detect = e.max_points - length(current_points)
    n_cell_detect = ceil(Int, n_detect / n_cells)

    features = CartesianIndex{2}[]
    sizehint!(features, n_detect)

    for y in 0:(e.grid_resolution[1] - 1), x in 0:(e.grid_resolution[2] - 1)
        y_shift, x_shift = y * e.cell_size, x * e.cell_size
        y_range = (y_shift + 1):min(height, (y + 1) * e.cell_size)
        x_range = (x_shift + 1):min(width, (x + 1) * e.cell_size)

        sub_features, _, n_detected = _shi_tomasi(
            @view(image[y_range, x_range]), n_cell_detect)
        n_detected == 0 && continue
        for sk in Keypoints(sub_features)
            push!(features, CartesianIndex{2}(sk[1] + y_shift, sk[2] + x_shift))
        end
    end

    features
end

"""
# Returns:

Vector of descriptors of `BitVector` type
and vector of keypoints' coordinates of type `CartesianIndex{2}`.
"""
function describe(e::Extractor, image, keypoints)
    create_descriptor(image, keypoints, e.descriptor)
end

"""
Create a mask from a set of points.

Each point defines a circular region that is used to avoid
detecting features in that region.

Circular regions are filled with `zero(T)` values,
while the rest of the mask is has `ones(T)` values.
"""
function get_mask(image::Matrix{T}, points, radius) where T
    mask = ones(T, size(image))
    for point in points
        draw!(mask, CirclePointRadius(to_cartesian(point), radius), zero(T))
    end
    mask
end


text:
using OrdinaryDiffEq, ParameterizedFunctions, Plots, ODE, ODEInterfaceDiffEq, LSODA, DiffEqDevTools, Sundials
using LinearAlgebra
LinearAlgebra.BLAS.set_num_threads(1)

gr() #gr(fmt=:png)

f = @ode_def Hires begin
  dy1 = -1.71*y1 + 0.43*y2 + 8.32*y3 + 0.0007
  dy2 = 1.71*y1 - 8.75*y2
  dy3 = -10.03*y3 + 0.43*y4 + 0.035*y5
  dy4 = 8.32*y2 + 1.71*y3 - 1.12*y4
  dy5 = -1.745*y5 + 0.43*y6 + 0.43*y7
  dy6 = -280.0*y6*y8 + 0.69*y4 + 1.71*y5 -
           0.43*y6 + 0.69*y7
  dy7 = 280.0*y6*y8 - 1.81*y7
  dy8 = -280.0*y6*y8 + 1.81*y7
end

u0 = zeros(8)
u0[1] = 1
u0[8] = 0.0057
prob = ODEProblem(f,u0,(0.0,321.8122))

sol = solve(prob,Rodas5(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)
abstols = 1.0 ./ 10.0 .^ (4:11)
reltols = 1.0 ./ 10.0 .^ (1:8);

plot(sol)

plot(sol,tspan=(0.0,5.0))

#sol = solve(prob,ode23s()); println("Total ODE.jl steps: $(length(sol))")
#using GeometricIntegratorsDiffEq
#try
#    sol = solve(prob,GIRadIIA3(),dt=1/10)
#catch e
#    println(e)
#end

setups = [
#Dict(:alg=>ROCK2()),
#Dict(:alg=>ROCK4())
#Dict(:alg=>ESERK5())
]

abstols = 1.0 ./ 10.0 .^ (5:8)
reltols = 1.0 ./ 10.0 .^ (1:4);
setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>Rodas3()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>rodas()),
          Dict(:alg=>radau()),
          Dict(:alg=>RadauIIA5()),
          Dict(:alg=>ROS34PW1a()),
          Dict(:alg=>lsoda()),
          ]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;dense = false,verbose=false,
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>Kvaerno3()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>KenCarp3()),
    # Dict(:alg=>SDIRK2()), # Removed because it's bad
          Dict(:alg=>radau())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;dense = false,verbose=false,
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>KenCarp5()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>KenCarp3()),
          Dict(:alg=>ARKODE(order=5)),
          Dict(:alg=>ARKODE()),
          Dict(:alg=>ARKODE(order=3))]
names = ["Rosenbrock23" "KenCarp5" "KenCarp4" "KenCarp3" "ARKODE5" "ARKODE4" "ARKODE3"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      names=names,save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;dense = false,verbose=false,
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>ImplicitEulerExtrapolation()),
          #Dict(:alg=>ImplicitDeuflhardExtrapolation()), # Diverges
          #Dict(:alg=>ImplicitHairerWannerExtrapolation()), # Diverges
          Dict(:alg=>ABDF2()),
          Dict(:alg=>QNDF()),
          Dict(:alg=>Exprb43()),
          Dict(:alg=>Exprb32()),
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)

abstols = 1.0 ./ 10.0 .^ (7:13)
reltols = 1.0 ./ 10.0 .^ (4:10)

setups = [Dict(:alg=>GRK4A()),
          Dict(:alg=>Rodas4P()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>ddebdf()),
          Dict(:alg=>Rodas5()),
          Dict(:alg=>rodas()),
          Dict(:alg=>radau()),
          Dict(:alg=>lsoda()),
          Dict(:alg=>RadauIIA5()),
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
                      dense=false,appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2)
plot(wp)

setups = [
          Dict(:alg=>Rodas5()),
          Dict(:alg=>Kvaerno5()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>Rodas4()),
          Dict(:alg=>radau())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
                      dense=false,appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>KenCarp5()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>KenCarp3()),
          Dict(:alg=>ARKODE(order=5)),
          Dict(:alg=>ARKODE()),
          Dict(:alg=>ARKODE(order=3))]
names = ["Rosenbrock23" "KenCarp5" "KenCarp4" "KenCarp3" "ARKODE5" "ARKODE4" "ARKODE3"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                      names=names,save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
                      dense=false,appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)

#setups = [#Dict(:alg=>Hairer4()),
          #Dict(:alg=>Hairer42()),
          #Dict(:alg=>Rodas3()),
          #Dict(:alg=>Kvaerno4()),
          #Dict(:alg=>KenCarp5()),
          #Dict(:alg=>Cash4())
#]
#wp = WorkPrecisionSet(prob,abstols,reltols,setups;
#                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
#plot(wp)

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
using ActuaryUtilities

cfs = [-7000,11.33137357858057, 11.167072336683004, 14.73988904111101, 14.543375011974934, 14.787026189603953, 12.220165229291128, 10.62228303189419, 10.2054439968351, 12.38830768146916, 11.55662975000198, 10.328454158492287, 11.12038681607141, 10.324480820201908, 10.198382859550456, 11.622094119543757, 12.273321426327541, 13.968817009905088, 14.633381040913772, 14.653117743106158, 14.75415942102587, 11.115621653400536, 12.491301011838484, 11.081360378007831, 12.332972598070072, 10.880589146045264, 10.235136200262108, 12.505550229750812, 11.250354442977768, 13.007688935886225, 10.949658164170005, 10.327509073165091, 12.081824327395134, 11.095672179619607, 12.879864014703418, 11.686291918815805, 10.487486433970597, 10.782904021875865, 10.828899614147883, 13.740734146503733, 13.070147160401934, 11.281533610255876, 14.608845949786723, 11.101609237726503, 13.6672383716922, 11.991725804610223, 14.223338394270746, 12.54775932654576, 11.669406707832616, 11.172531062630963, 14.797739369742395, 14.742002039355507, 11.26059357755014, 10.005760731173147, 13.166395327031115, 10.953541497765157, 13.85589184547591, 11.938616564984416, 12.602780339084227, 10.190660874910725, 12.62253014530776, 13.25618155360879, 11.909328696772635, 13.645695694839475, 10.771692737595489, 10.019159704029537, 12.748198733268266, 14.949945651485546, 10.471989762364004, 11.694659010236517, 10.840097149970063, 13.81959526563658, 13.725232362048477, 14.696971575480223, 14.061440252466008, 11.894097203077104, 11.623559529246583, 11.422764616449735, 11.911568770804685, 14.257873178733874, 11.908786512287042, 11.704867143116388, 13.547102659825594, 10.231510523063077, 12.580425576504725, 10.204318593593095, 10.839329475371066, 14.06382079409852, 14.741551103697574, 14.125100696255025, 10.534750441149983, 12.047335974087941, 12.040901005285527, 13.061055670703222, 12.033900288651994, 14.64914110907408, 14.74150289957578, 14.18198082748481, 12.64569352877869, 10.795441413087191, 12.32677569400135, 14.201284184243757, 12.198731964230793, 13.829237277802761, 12.88556220895728, 13.681945290090994, 11.310817183394416, 10.642344563381029, 12.001334393089806, 13.940994758512655, 14.588579514241378, 11.784013809882609, 14.857525454933082, 11.890304205798119, 14.163978789672555, 12.802006519405316, 10.453613727321805, 10.545421710663241, 13.860350461267206, 13.011623286218281, 12.529578764153623, 12.0682203481126, 12.222973006731433, 10.782217823800455, 14.815112126486351, 14.021264432689495, 13.907504381640877, 13.263889558252384, 13.886523931472032, 14.92463305217065, 10.39804588193537, 11.717727995106285, 13.475830384363098, 13.96526870747185, 12.537688951669926, 10.713851417257366, 12.651783539839993, 14.94145019017787, 11.298228878804299, 11.238462723076646, 12.54816626881003, 12.991279893628143, 12.938650960410705, 14.509841908929733, 13.061627992570637, 13.532302333363827, 10.87538595312562, 13.401986714313082, 14.452077231881011, 10.62775390563043, 10.259969609401132, 10.630066121172158, 10.67746807962661, 14.301509125704097, 12.248186770125802, 13.529935347104352, 12.214180545525151, 13.737238621655871, 14.935337045307609, 10.715008490421852, 13.792764627232366, 11.566835806502482, 13.381623410224446, 11.759551394829062, 14.46722368188435, 14.820860467933276, 14.367067474219198, 12.785099227552985, 12.125131138279556, 13.59004374972493, 11.035277329107453, 13.552102778902684, 11.167577469887052, 12.079972349457684, 14.78489299740463, 11.117948297108185, 11.81118834164748, 13.306337834666962, 14.61353141178894, 14.277193746539648, 12.166555805663473, 11.265160163436654, 10.362543406467916, 11.17487038814809, 13.221576830876275, 11.124927016257411, 11.778833269702968, 12.884165622590684, 11.70614082816003, 13.732507964152038, 13.562497475153835, 14.365673914585269, 11.218499437025232, 14.313161434674605, 10.84553935633242, 13.892594863691276, 12.905232048487619, 14.526707951361043, 11.151370775543134, 11.202538557418025, 14.351557380104524, 14.618321803616304, 13.450749033215105, 13.971103464192872, 12.453728166838555, 11.203648664519335, 13.354317055813262, 14.559195884990523, 13.853813999405409, 10.265888555997426, 11.740493314590367, 13.565808716032713, 14.465381469251845, 11.700055034584512, 13.666543210894279, 10.05828768839079, 13.793893081997751, 11.840407775390236, 14.292054186783561, 13.327715352294016, 13.176691551807002, 14.385175692391908, 14.385014460873183, 11.407277293369141, 14.231388088664055, 12.250059232161178, 14.86235229497673, 13.77507647192779, 12.896049896482687, 10.60002301502453, 13.447070421155189, 13.358321797846632, 11.21916910048809, 13.30255271341071, 13.493253069737648, 10.948637495146045, 10.001718613821161, 13.91551521108937, 11.581689906353411, 14.977887214372007, 14.65693440845071, 14.665474268673329, 12.691023502305779, 11.636579739013367, 13.145065301442465, 14.054029882788013, 12.872698520853412, 14.603129068297749, 14.939572657418722, 11.170995199160744, 11.847096639046411, 13.883460710012764, 10.126422540987893, 13.407267339989394, 10.421297015199778, 14.793530503055097, 14.009521709723277, 10.551203664546083, 13.606121157647152, 13.893115513324952, 12.387229455207336, 10.152726271083313, 10.857363093853683, 10.553311047158111, 14.48112608836534, 12.009660541892252, 10.562740549914857, 10.185345168077752, 10.610848424790655, 12.382686150676019, 12.70009488541087, 11.751661090276412, 11.18298540002945, 10.938505981462802, 14.11088111794641, 14.40865642716125, 12.87387852504945, 11.345425461857499, 13.450458614928777, 14.558969761971529, 11.784359157402896, 11.528844481965896, 12.841302594496753, 10.516217100119583, 12.034691758196473, 10.631794257379834, 14.71733657041019, 14.820479028358028, 14.780542929577692, 13.052248672850174, 13.437337644194365, 12.7734425699407, 12.709862021282632, 10.034620588718868, 10.850848229793126, 14.536407975610807, 11.155751652337214, 12.193467364682604, 13.089842912540833, 13.093755843694804, 13.963076758897023, 10.782746529071671, 10.57634443951439, 10.585652732773541, 12.388022899491169, 11.083235905164493, 11.807420712958251, 11.670764084958904, 14.343548642031877, 14.089822252264673, 11.146774006554228, 14.83903097303459, 14.644353374427286, 13.244758181211573, 12.410509038164601, 11.643494403195634, 14.092682065844572, 14.367586298974796, 14.810106489623035, 11.945040477569586, 14.401686777111255, 11.887093935272377, 14.462898789116835, 10.226154690178888, 13.016393381958249, 13.453915585803575, 13.578554072849135, 14.335352027320045, 12.265010406231234, 10.757058969921816, 11.876531697073972, 12.0491828922183, 11.092453908216743, 12.109105293560898, 12.938024069917983, 14.208843596502678, 11.466720564074699, 12.308443242742026, 13.1609372489227, 11.65441215195797, 13.728800043792663, 13.847762282971136, 14.87752275142158, 11.149184723288272, 13.640462888204215, 14.237264786466763, 13.405122090695556, 12.80566689508301, 14.483299214044472, 11.767922175746978, 13.53510927127212, 13.684872748818286, 12.553957482155372, 11.004089018712035, 14.520634498250335, 13.83809738781293, 10.672326920794527, 11.419480324115304, 14.460878992784673, 11.511794268770446, 14.47105838704902, 14.935713252855354, 13.024901111795435, 12.111956580271531, 10.646741917366109, 14.072054873809176, 12.094227297375376, 13.981650144730342, 10.648895767446431, 11.571433949013413, 14.661395278086946, 11.868831686617465, 12.875727107393082, 10.771518408030383, 13.631011657679027, 13.11336224401068, 12.145472806863014, 11.991080661902533, 12.130009173958825, 13.515256179224657, 11.513137707370925, 13.754468506719768, 14.913179336936315, 13.58791955747601, 10.134569430875423, 13.134677899479813, 12.183005498224883, 14.656674018446617, 13.657413280225123, 11.142466962252714, 11.105613647573488, 10.890935137974486, 13.660316392058583, 11.992634871732069, 14.039474798844504, 10.10273673903627, 12.512433682080895, 14.791392396799488, 13.952028496884394, 10.68072002921823, 14.148714257563924, 14.160473668158666, 10.990107339106551, 10.34171121730688, 13.669669138346311, 12.290862020072826, 10.272777155525883, 14.966454869056701, 10.300288608313576, 10.737680728398, 10.630967159879477, 13.809878958904225, 14.447266368588924, 14.456981802538502, 10.906800282148177, 13.631250425690478, 10.65625820412889, 12.765564689360023, 13.229159965330654, 11.490266136800688, 13.964187960877286, 13.923896694765084, 13.29754185129391, 13.337423023172482, 10.920180637381083, 14.69321208846982, 12.737081695738025, 12.086001267790907, 14.63420753638875, 12.995937177677611, 14.845264104376177, 13.172941798235495, 14.00450081411807, 14.486240101939611, 12.575356268048141, 10.572728182169437, 11.817972894521603, 11.020289503530373, 12.33588646221233, 14.38354993054751, 12.48657836304595, 14.045752344056417, 14.291881431404809, 11.449852919176939, 12.285269949280199, 13.329686561081914, 12.392607450843016, 11.421898690830588, 13.033623164214463, 13.740312988118928, 14.86808204262957, 14.052280858899593, 11.857022479258601, 13.344298820910902, 12.815321121685948, 12.542919644519769, 10.221542848537679, 12.825606712087758, 13.80917375726888, 11.873087750288265, 13.192060469426856, 10.513049332061176, 13.71129267069243, 10.669450949394452, 14.380201616255327, 13.516014153886436, 10.239185746011342, 14.183752673275308, 11.003573730151318, 10.824874699517318, 14.3140211605076, 10.735839396867446, 11.829061949174482, 12.444084226653683, 11.660344546867854, 14.613122657676357, 11.117541787276522, 13.71389482775391, 11.844802667361172, 13.809064011020476, 12.545980413077045, 10.37221231321253, 11.099185606018425, 12.069367607148795, 14.753997650830298, 14.860064867206237, 11.451965081962701, 12.408066998699887, 10.09219363945688, 14.804046243343, 14.462425372901544, 10.396509478003066, 13.465502211253604, 12.193256544832277, 10.419849455572237, 12.335938167870554, 14.077332364067836, 12.04238111937329, 11.219586684431405, 11.191981905319334, 13.990789915195073, 13.321355746701418, 13.867997027534281, 11.34094094573737, 11.894860840187725, 10.014643905455108, 13.9450805524248, 11.369872605476774, 13.823701043903593, 11.632450588220253, 14.402129467913422, 10.851529252537635, 13.446334384773026, 12.881948367117989, 14.301072054230792, 13.901034543377902, 14.702417172046475, 13.483519738663922, 12.257856089368287, 12.736794695148546, 11.151630919820665, 11.870685024867045, 13.667266606426, 12.666123528665413, 10.58046880567575, 12.669966690049995, 12.850748247740196, 11.685094024655385, 10.84736035947657, 10.879550786172686, 14.583912720522537, 11.556576292636521, 14.446763862373635, 12.2747770021459, 12.397640643976004, 10.138635782380568, 10.19699771111878, 14.33115980578789, 13.069369589587058, 11.654396130820546, 14.32588536103658, 11.234253581651313, 10.886330102178094, 13.905356191291201, 11.31211292200017, 14.591797165769588, 11.57906280643942, 12.687157601724977, 14.231068598718757, 11.533307917760986, 11.00278786406948, 11.791374923466336, 12.355367201523464, 13.412638962511421, 14.059784911142893, 14.929574141464332, 13.985148455547298, 13.478658632335467, 12.088803020571913, 11.076902880948095, 14.85282953089796, 13.557554489783813, 14.09449869180156, 11.058526150649634, 11.871607439321343, 13.416916873614177, 11.896264445342519, 12.140871542725174, 13.967040028548759, 13.92387024388108, 12.557100140784577, 13.709689153237505, 11.455617154148374, 11.983798170601936, 12.593738595063888, 10.852961452174824, 10.891042640728921, 11.018980720091063, 10.661661557734767, 11.968388540103959, 14.575380826864798, 10.310664116233681, 10.79793225618651, 13.175129120149073, 11.868037287978872, 12.428601245011638, 12.297281024758835, 11.75532988335384, 11.935023303980994, 11.445578716193893, 11.635791018694775, 12.229231977915058, 14.947200554894952, 10.584576916849892, 13.543547215152838, 10.948422942905365, 10.085858606874586, 11.40453743267752, 14.389263873891274, 14.270628905078073, 10.31802746298273, 11.057013004827654, 12.128662408389431, 13.19152229812298, 11.84574262866243, 12.066273890508961, 12.7038048452766, 10.95981080906325, 10.588705330230361, 12.890457480879702, 14.110412467286706, 14.63310013295482, 13.421448571071592, 11.873600277480264, 11.823406873514843, 14.649020650603521, 11.535855502475416, 14.986616581388084, 10.854693292973192, 11.738250952256035, 14.731546956588033, 10.867566954581275, 12.862704113181548, 11.185985864878273, 12.661047000348495, 10.959402002641353, 13.85593485994154, 12.727773370880131, 13.270188110208368, 12.155893411161436, 11.251995770990543, 10.734028201769819, 13.341862536135448, 12.56663295008742, 14.829248071219354, 12.63146009202443, 12.087984402629713, 12.687332073795645, 13.024881500305945, 14.636625631024831, 14.330939418774697, 14.489727988032627, 10.894975953362461, 14.063442584619581, 13.335748722823846, 14.54081960859539, 13.094259093921789, 13.25815397434542, 11.534888981037271, 11.52402384179464, 14.566248514539671, 14.896673567000065, 11.976985722167944, 10.091588801010063, 13.181522254786389, 14.03106487342248, 12.406847857913657, 12.528017683354857, 11.269268754233366, 12.577593443628423, 14.062655963542348, 14.364172468278504, 12.919707007699827, 14.460732124630018, 14.421854029205928, 14.342538815876301, 13.715570088998515, 11.952513510277187, 10.29676956281684, 10.149818255362092, 13.253896447552147, 14.423345225432362, 12.268963697091056, 14.286774548972177, 11.905175204077148, 13.493744141319985, 13.637722653600923, 14.541889104088948, 14.379843837582854, 11.916040077750571, 11.358478642311324, 14.07020030698621, 10.887539603342036, 14.906750067630346, 13.68689407511227, 11.578510345740794, 10.204256090759499, 12.512830804014698, 13.182658964969256, 12.784048294940208, 10.54838424521611, 13.516015586421368, 11.533499652630915, 11.455410070235674, 10.200545340237175, 13.422586005304904, 11.917659998933392, 13.338723186085131, 12.707389546708027, 13.372196887397864, 14.15110846956745]

@benchmark irr($cfs)

max_stars_repo_name: JeannotJeannot/Learn

max_issues_repo_name: JeannotJeannot/Learn

max_forks_repo_name: JeannotJeannot/Learn


text:
function tenrnet(X::AbstractMatrix{T}; donorm::Bool=true) where T<:Real
    ℊ,𝒸=size(X)
    if donorm
        lbsz=sum(X,dims=1)
        # X=(X./lbsz)*median(lbsz)
        X=(X./lbsz)*1e4
        X=log.(X.+1)
    end
    A=zeros(Float16, ℊ, ℊ, NLAYERS)
    for k=1:NLAYERS
        println("network ... $k")
        𝕩=X[:,randperm(𝒸)][:,1:NCELLS]        # jackknife (m-out-of-n)
        # 𝕩=X[:,rand(1:𝒸,NCELLS)];            # bootstrapping (m-out-of-n)
        𝕩ᵀ=transpose(𝕩)
        a=pcnet(𝕩ᵀ,NCOMP1)
        a[abs.(a).<quantile(vec(abs.(a)),0.95)].=0.0
        @inbounds A[:,:,k]=sparse(a)
    end
    Z=tensordecomp(A,NCOMP2)
    return Z,A
end

max_stars_repo_name: jamesjcai/ScTenifold

max_issues_repo_name: jamesjcai/ScTenifold

max_forks_repo_name: jamesjcai/ScTenifold


text:
using LoopVectorization
using Tullio

"""
    naive_lpd(log_likelihood::AbstractArray{<:Real}, chain_index::Vector}<:Int)
Calculate the naive (in-sample) estimate of the expected log probability density, otherwise
known as the in-sample Bayes score. Not recommended for most uses.

# Arguments
  - $LIKELIHOOD_ARRAY_ARG
  - $CHAIN_INDEX_DOC
"""
function naive_lpd(log_likelihood::AbstractArray{<:Real, 3})
    @info "We advise against using `naive_lpd`, as it gives inconsistent and strongly " *
    "biased estimates. Use `psis_loo` instead."
    return _naive_lpd(log_likelihood)
end

function naive_lpd(
    log_likelihood::AbstractMatrix{<:Real},
    chain_index::AbstractVector{<:Integer} = _assume_one_chain(log_likelihood)
)
    @nospecialize(chain_index)
    log_likelihood = _convert_to_array(log_likelihood, chain_index)
    return _naive_lpd(log_likelihood)
end

function _naive_lpd(log_likelihood::AbstractArray{<:Real, 3})
    dims = size(log_likelihood)
    data_size = dims[1]
    mcmc_count = dims[2] * dims[3]  # total number of samples from posterior
    log_count = log(mcmc_count)

    pointwise_naive = similar(log_likelihood, data_size)
    @tullio pointwise_naive[i] = exp(log_likelihood[i, j, k] - log_count) |> log
    return @tullio naive := pointwise_naive[i]
end


text:
using .ShessMark, .StatBack

function mark_chess(robot)
    init(isodd(sum(move_to_start!((Sud,West)))) # isodd - возвращает true, если аргумент нечетный
    #УТВ: Робот - в юго-западном углу и глобальная переменная модуля инициализирована

    side = Ost
    mark_chess(robot,side)
    while isborder(Nord)==false
        move!(robot,Nord)
        side = inverse(side)
        mark_chess(robot,side)
    end
    #УТВ: Робот - у северной границы поля И маркеры расставлены в шахматном порядке

    for side in (West,Sud) movements!(robot,side) end
    #УТВ: Робот - в юго-западном углу

    move_to_back(robot)
    #УТВ: Робот - в исходном положении
end

function mark_chess(r,side)
    ShessMark.putmarker!(r) # здесь вызывается не стандартная команда Робота, а одноименная функция из модуля
    while isborder(r,side)==false
        move!(r,side)
        ShessMark.putmarker!(r)
    end
end


text:
VERSION >= v"0.4.0-dev+6641" && __precompile__()
module StrPack

export @struct
export pack, unpack, sizeof
export DataAlign
export align_default, align_packed, align_packmax, align_structpack, align_table
export align_x86_pc_linux_gnu, align_native
export show_struct_layout

using Base.Meta
using Compat

import Base.read, Base.write
import Base.isequal

bswap(c::Char) = identity(c) # white lie which won't work for multibyte characters in UTF-16 or UTF-32

# Represents a particular way of adding bytes to maintain certain alignments
immutable DataAlign
    ttable::Dict
    # default::(Type -> Integer); used for bits types not in ttable
    default::Function
    # aggregate::(Vector{Type} -> Integer); used for composite types not in ttable
    aggregate::Function
end
DataAlign(def::Function, agg::Function) = DataAlign((@compat Dict{Type,Integer}()), def, agg)

immutable Struct
    asize::Dict
    strategy::DataAlign
    endianness::Symbol
end

macro struct(xpr...)
    (typname, typ, asize) = extract_annotations(xpr[1])
    if length(xpr) > 3
        error("too many arguments supplied to @struct")
    end
    if length(xpr) > 2
        if isexpr(xpr[3], :quote) && haskey(endianness_converters, eval(xpr[3]))
            endianness = xpr[3]
        else
            error("$(string(xpr[3])) is not a valid endianness")
        end
        alignment = xpr[2]
    elseif length(xpr) > 1
        if isexpr(xpr[2], :quote) && haskey(endianness_converters, eval(xpr[2]))
            endianness = xpr[2]
            alignment = :(align_default)
        else
            alignment = xpr[2]
            endianness = :(:NativeEndian)
        end
    else
        alignment = :(align_default)
        endianness = :(:NativeEndian)
    end
    new_struct = :(Struct($asize, $alignment, $endianness))
    quote
        $(esc(typ))
        $(esc(:(isdefined(:STRUCT_REGISTRY) || const STRUCT_REGISTRY = ObjectIdDict())))
        $(esc(:(STRUCT_REGISTRY[$typname]))) = $new_struct
        const fisequal = isequal
        fisequal(a::$(esc(typname)), b::$(esc(typname))) = begin
            for name in fieldnames($(esc(typname)))
                if !isequal(getfield(a, name), getfield(b, name))
                    return false
                end
            end
            true
        end
    end
end

headof(xpr, sym) = false
headof(xpr::Expr, sym) = (xpr.head == sym)

# hmm...AST selectors could be useful
function extract_annotations(exprIn)
    fieldnames = Expr[]
    asizes = Expr[]
    typname = nothing
    if isexpr(exprIn, :type)
        if isa(exprIn.args[2],Symbol)
            typname = exprIn.args[2]
        elseif isexpr(exprIn.args[2],:<:)
            typname = exprIn.args[2].args[1]
        else
            error("Unable to extract type name!")
        end
        for field_xpr in exprIn.args[3].args
            if isexpr(field_xpr, :(::)) && isexpr(field_xpr.args[2], :call)
                push!(fieldnames, quot(field_xpr.args[1]))
                push!(asizes, quot(Integer[field_xpr.args[2].args[2:end]...]))
                # turn this back into a legal type expression
                field_xpr.args[2] = field_xpr.args[2].args[1]
            end
        end
    else
        error("only type definitions can be supplied to @struct")
    end
    asize = :(Dict(zip([$(fieldnames...)], Array{Integer,1}[$(asizes...)])))
    (typname, exprIn, asize)
end

endianness_converters = @compat Dict(
    :BigEndian => (hton, ntoh),
    :LittleEndian => (htol, ltoh),
    :NativeEndian => (identity, identity),
    :SwappedEndian => (Base.bswap, Base.bswap))

# A byte of padding
bitstype 8 PadByte
write(s::IO, x::PadByte) = write(s, 0x00)
read(s::IO, ::Type{PadByte}) = read(s, UInt8)

function isbitsequivalent{T}(::Type{T})
    if isbits(T) || T <: AbstractString && !T.abstract
        return true
    elseif isempty(fieldnames(T))
        return false
    end
    # TODO AbstractArray inspect element instead
    for S in T.types
        if S <: AbstractArray
            S = eltype(S)
        end
        if !isbitsequivalent(S)
            return false
        end
    end
    true
end

function chktype{T}(::Type{T})
    if isempty(fieldnames(T))
        error("Type $T is not an aggregate type.")
    end
    if !isbitsequivalent(T)
        error("Type $T is not bits-equivalent.")
    end
end

function unpack{T}(in::IO, ::Type{T}, asize::Dict, strategy::DataAlign, endianness::Symbol)
    chktype(T)
    tgtendianness = endianness_converters[endianness][2]
    offset = 0
    rvar = Any[]
    for (typ, name) in zip(T.types, fieldnames(T))
        dims = get(asize, name, 1)
        intyp = if typ <: AbstractArray
            eltype(typ)
        else
            typ
        end

        # Skip padding before next field
        pad = pad_next(offset,intyp,strategy)
        skip(in,pad)
        offset += pad
        offset += if intyp <: AbstractString
            push!(rvar, rstrip(convert(typ, read(in, UInt8, dims...)), ['\0']))
            prod(dims)
        elseif !isempty(fieldnames(intyp))
            if typ <: AbstractArray
                item = Array(intyp, dims...)
                for i in 1:prod(dims)
                    item[i] = unpack(in, intyp)
                end
                push!(rvar, item)
            else
                push!(rvar, unpack(in, intyp))
            end
            calcsize(intyp)*prod(dims)
        else
            if typ <: AbstractArray
                push!(rvar, map(tgtendianness, read(in, intyp, dims...)))
            else
                push!(rvar, tgtendianness(read(in, intyp)))
            end
            sizeof(intyp)*prod(dims)
        end
    end
    skip(in, pad_next(offset, T, strategy))
    T(rvar...)
end
function unpack{T}(in::IO, ::Type{T}, endianness::Symbol)
    chktype(T)
    reg = T.name.module.STRUCT_REGISTRY[T]::Struct
    unpack(in, T, reg.asize, reg.strategy, endianness)
end
function unpack{T}(in::IO, ::Type{T})
    chktype(T)
    reg = T.name.module.STRUCT_REGISTRY[T]::Struct
    unpack(in, T, reg.asize, reg.strategy, reg.endianness)
end

function pack{T}(out::IO, struct::T, asize::Dict, strategy::DataAlign, endianness::Symbol)
    chktype(T)
    tgtendianness = endianness_converters[endianness][1]
    offset = 0
    for (typ, name) in zip(T.types, fieldnames(T))
        if typ <: AbstractArray
            typ = eltype(typ)
        end
        data = if typ <: AbstractString
            typ = UInt8
            convert(Array{UInt8}, getfield(struct, name))
        else
            getfield(struct, name)
        end

        offset += write(out, zeros(UInt8, pad_next(offset, typ, strategy)))

        numel = prod(get(asize, name, 1))
        idx_end = numel > 1 ? min(numel, length(data)) : 1
        if !isempty(fieldnames(typ))
            if typeof(data) <: AbstractArray
                for i in 1:idx_end
                    offset += pack(out, data[i])
                end
                offset += write(out, zeros(UInt8, calcsize(typ)*(numel-idx_end)))
            else
                offset += pack(out, data)
            end
        else
            offset += if typeof(data) <: AbstractArray
                isempty(data) ? 0 : write(out, map(tgtendianness, data[1:idx_end]))
            else
                write(out, tgtendianness(data))
            end
            offset += write(out, zeros(typ, max(numel-idx_end, 0)))
        end
    end
    offset += write(out, zeros(UInt8, pad_next(offset, T, strategy)))
end
function pack{T}(out::IO, struct::T, endianness::Symbol)
    chktype(T)
    reg = T.name.module.STRUCT_REGISTRY[T]
    pack(out, struct, reg.asize, reg.strategy, endianness)
end
function pack{T}(out::IO, struct::T)
    chktype(T)
    reg = T.name.module.STRUCT_REGISTRY[T]
    pack(out, struct, reg.asize, reg.strategy, reg.endianness)
end

# Convenience methods when you just want to use strings
macro withIOBuffer(iostr, ex)
    quote
        $iostr = IOBuffer()
        $ex
        $iostr
    end
end

pack{T}(struct::T, a::Dict, s::DataAlign, n::Symbol) = @withIOBuffer iostr pack(iostr, a, s, n)
pack{T}(struct::T) = @withIOBuffer iostr pack(iostr, struct)

unpack{T}(str::Union{AbstractString, Array{UInt8,1}}, ::Type{T}) = unpack(IOBuffer(str), T)

## Alignment strategies and utility functions ##

# default alignment for bitstype T is nextpow2(sizeof(::Type{T}))
type_alignment_default{T<:AbstractArray}(::Type{T}) = type_alignment_default(eltype(T))
type_alignment_default{T<:AbstractString}(::Type{T}) = 1
type_alignment_default{T}(::Type{T}) = nextpow2(sizeof(T))

# default strategy
align_default = DataAlign(type_alignment_default, x -> maximum(map(type_alignment_default, x)))

# equivalent to __attribute__ (( __packed__ ))
align_packed = DataAlign(_ -> 1, _ -> 1)

# equivalent to #pragma pack(n)
align_packmax(da::DataAlign, n::Integer) = DataAlign(
    da.ttable,
    _ -> min(type_alignment_default(_), n),
    da.aggregate,
    )

# equivalent to __attribute__ (( align(n) ))
align_structpack(da::DataAlign, n::Integer) = DataAlign(
    da.ttable,
    da.default,
    _ -> n,
    )

# provide an alignment table
align_table(da::DataAlign, ttable::Dict) = DataAlign(
    merge(da.ttable, ttable),
    da.default,
    da.aggregate,
    )

# convenience forms using a default alignment
for fun in (:align_packmax, :align_structpack, :align_table)
    @eval ($fun)(arg) = ($fun)(align_default, arg)
end

# Specific architectures
align_x86_pc_linux_gnu = align_table(align_default,
    @compat Dict(
    Int64 => 4,
    UInt64 => 4,
    Float64 => 4,
    ))

# Get alignment for a given type
function alignment_for(strategy::DataAlign, T::Type)
    if haskey(strategy.ttable, T)
        strategy.ttable[T]
    elseif !isempty(fieldnames(T))
        strategy.aggregate(T.types)
    else
        strategy.default(T)
    end
end

function pad_next(offset, typ, strategy::DataAlign)
    align_to = alignment_for(strategy, typ)
    (align_to - offset % align_to) % align_to
end

function calcsize{T}(::Type{T}, asize::Dict, strategy::DataAlign)
    chktype(T)
    size = 0
    for (typ, name) in zip(T.types, fieldnames(T))
        dims = get(asize, name, 1)
        typ = if typ <: Array
            eltype(typ)
        elseif typ <: AbstractString
            UInt8
        else
            typ
        end
        size += pad_next(size, typ, strategy)
        size += if isbits(typ)
            prod(dims)*sizeof(typ)
        elseif !isempty(fieldnames(typ))
            prod(dims)*sizeof(Struct(typ))
        else
            error("Improper type $typ in struct.")
        end
    end
    size += pad_next(size, T, strategy)
    size
end
calcsize{T}(::Type{T}) = calcsize(T, T.name.module.STRUCT_REGISTRY[T].asize, T.name.module.STRUCT_REGISTRY[T].strategy)

function show_struct_layout{T}(::Type{T}, asize::Dict, strategy::DataAlign, width, bytesize)
    chktype(T)
    offset = 0
    for (typ, name) in zip(T.types, fieldnames(T))
        dims = get(asize, name, 1)
        intyp = if typ <: Array
            eltype(typ)
        elseif typ <: AbstractString
            UInt8
        else
            typ
        end
        padsize = pad_next(offset, intyp, strategy)
        offset = show_layout_format(PadByte, sizeof(PadByte), padsize, width, bytesize, offset)
        offset = show_layout_format(typ, sizeof(intyp), dims, width, bytesize, offset)
    end
    padsize = pad_next(offset, T, strategy)
    offset = show_layout_format(PadByte, sizeof(PadByte), padsize, width, bytesize, offset)
    if offset % width != 0
        println()
    end
end
show_struct_layout{T}(::Type{T}) = show_struct_layout(T, T.name.module.STRUCT_REGISTRY[T].asize, T.name.module.STRUCT_REGISTRY[T].strategy, 8, 10)
show_struct_layout{T}(::Type{T}, width::Integer, bytesize::Integer) = show_struct_layout(T, T.name.module.STRUCT_REGISTRY[T].asize, T.name.module.STRUCT_REGISTRY[T].strategy, width, bytesize)
# show_struct_layout(T::Type, asize::Dict, strategy::DataAlign, width) = show_struct_layout(T, strategy, width, 10)

function show_layout_format(typ, typsize, dims, width, bytesize, offset)
    for i in 1:prod(dims)
        tstr = string(typ)[1:min(typsize*bytesize-2, end)]
        str = "[" * tstr * "-"^(bytesize*typsize-2-length(tstr)) * "]"
        typsize_i = typsize
        while !isempty(str)
            if offset % width == 0
                @printf("0x%04X ", offset)
            end
            len_prn = min(width - (offset % width), typsize_i)
            nprint = bytesize*len_prn
            print(str[1:nprint])
            str = str[nprint+1:end]
            typsize_i -= len_prn
            offset += len_prn
            if offset % width == 0
                println()
            end
        end
    end
    offset
end

## Native layout ##
align_native = align_table(align_default, let
    i8a, i16a, i32a, i64a, f32a, f64a = Array(UInt, 1), Array(UInt, 1), Array(UInt, 1), Array(UInt, 1), Array(UInt, 1), Array(UInt, 1)

    ccall("jl_native_alignment", Void,
          (Ptr{UInt}, Ptr{UInt}, Ptr{UInt}, Ptr{UInt}, Ptr{UInt}, Ptr{UInt}),
          i8a, i16a, i32a, i64a, f32a, f64a)

    @compat Dict(
     Int8 => i8a[1],
     UInt8 => i8a[1],
     Int16 => i16a[1],
     UInt16 => i16a[1],
     Int32 => i32a[1],
     UInt32 => i32a[1],
     Int64 => i64a[1],
     UInt64 => i64a[1],
     Float32 => f32a[1],
     Float64 => f64a[1],
    )
end)

end


text:
using JetPack, Jets, Test

@testset "JopBlend" begin
    A = JopBlend(Float64, 128, [5,10,20])
    lhs,rhs = dot_product_test(A,rand(domain(A)),rand(range(A)))
    @test lhs ≈ rhs
end


text:
import AlphaZero.GI
using AlphaZero: AbstractSchedule, ConstSchedule
using GraphNeuralNetworks
using Random
using Crayons

o2ji(o, M, N) = CartesianIndices((M, N))[o]
ji2o(j, i, M, N) = LinearIndices((M, N))[j, i]

struct GameSpec <: GI.AbstractGameSpec
  M::Pair{AbstractSchedule,AbstractSchedule}
  N::Pair{AbstractSchedule,AbstractSchedule}
  A::Pair{AbstractSchedule,AbstractSchedule}
  P::Pair{UInt8,UInt8}
end

function generate_conjuctive_edges(rng::AbstractRNG, M, N, N_OPP, T, S)
  order = reduce(vcat, [(randperm(rng, M) .+ i * M) for i in 0:N-1])#order of operations
  target = zeros(UInt8, N_OPP + 1 + N)
  target[T] = T
  n = 0
  for i in 1:N*M
    if (i % M == 1) #first in row: s->
      target[S+n] = order[i]
    end
    if (i % M == 0)#last in row: ->t
      target[order[i]] = T
      n += 1
    else #propagate
      target[order[i]] = order[i+1]
    end
  end
  return target
end

function gen_done_time(p_time, conj_tar, start_edges, S, T)
  done_time = zeros(UInt16, S)
  for o in start_edges
    last_done_time = 0
    while (true)#propagate expected done time
      o = conj_tar[o]
      if (o == T)
        done_time[T] = max(done_time[T], last_done_time)
        break
      end
      last_done_time = done_time[o] = last_done_time + minimum(values(p_time[o]))
    end
  end
  return done_time
end

function gen_action_values(p_time, conj_tar, start_edges, S, T)
  src = Matrix{UInt8}(undef, 0, 3)
  tar = Matrix{UInt8}(undef, 0, 3)
  info = Matrix{UInt8}(undef, 0, 2)
  done_time = []
  for (i, e) in enumerate(start_edges)
    o = conj_tar[e]
    for (m, time) in p_time[o]
      self = length(done_time) + S + 1
      src = [src; [S self S]]
      tar = [tar; [self conj_tar[o] self]]
      info = [info; [m o]]
      done_time = [done_time; time]
    end
  end
  return src, tar, info, done_time
end

mutable struct GameEnv <: GI.AbstractGameEnv
  #Problem instance
  process_time::Vector{Dict{UInt8,UInt8}} #Index in Mxn
  M::UInt8
  N::UInt8
  N_OPP::UInt16
  T::UInt16
  S::UInt16
  UB::UInt16
  LB::UInt16
  conj_src::Vector{UInt8}
  conj_tar::Vector{UInt8}
  #State
  action_src::Matrix{UInt8}
  action_tar::Matrix{UInt8}
  action_info::Matrix{UInt8}
  disj_src::Vector{UInt8}
  disj_tar::Vector{UInt8}
  is_done::Vector{Bool}
  done_time::Vector{UInt16}
  action_done_time::Vector{UInt16}
  #Info
  prev_operation::Vector{UInt8}
  prev_machine::Vector{UInt8}
end

function GI.init(spec::GameSpec, itc::Int, rng::AbstractRNG)
  M = rand(rng, spec.M.first[itc]:spec.M.second[itc])
  N = rand(rng, spec.N.first[itc]:spec.N.second[itc])
  N_OPP = M * N #[nodes, T, S]
  T = M * N + 1
  S = M * N + 2
  gen_node() = Dict([UInt8(i) => rand(rng, spec.P.first:spec.P.second) for i in randperm(rng, M)[1:rand(rng, spec.A.first[itc]:spec.A.second[itc])]])
  p_time = [gen_node() for _ in 1:N_OPP]
  conj_tar = generate_conjuctive_edges(rng, M, N, N_OPP, T, S)
  start_edges_n = collect(0:N-1) .+ S
  start_edges_m = collect(0:M-1) .+ S
  node_done = gen_done_time(p_time, conj_tar, start_edges_n, S, T)
  action_src, action_tar, action_info, action_done = gen_action_values(p_time, conj_tar, start_edges_n, S, T)
  return GameEnv(
    #Problem instance
    p_time,
    M,
    N,
    N_OPP,
    T,
    S,
    sum(maximum.(values.(p_time))), #All worst operations in a row
    maximum(sum.([minimum.(values.(p_time[i*M+1:(i+1)*M])) for i in 0:N-1])),#Best operations parrallel
    [collect(1:N_OPP); T; S * ones(N)],
    conj_tar,
    #State
    action_src,
    action_tar,
    action_info,
    [],
    [],
    [falses(N_OPP); false; true], #[Operations, T, S]
    node_done,
    action_done,
    #Info
    start_edges_n,
    start_edges_m
  )
end

GI.init(spec::GameSpec, s) = GameEnv(
  #Static values
  s.process_time,
  s.M,
  s.N,
  s.N_OPP,
  s.T,
  s.S,
  s.UB,
  s.LB,
  s.conj_src,
  s.conj_tar,
  #Mutable values
  copy(s.action_src),
  copy(s.action_tar),
  copy(s.action_info),
  copy(s.disj_src),
  copy(s.disj_tar),
  copy(s.is_done),
  copy(s.done_time),
  copy(s.action_done_time),
  copy(s.prev_operation),
  copy(s.prev_machine),
)

GI.spec(g::GameEnv) = GameSpec(ConstSchedule(g.M) => ConstSchedule(g.M), ConstSchedule(g.N) => ConstSchedule(g.N), ConstSchedule(3) => ConstSchedule(3), 1 => 5)

GI.two_players(::GameSpec) = false

GI.state_dim(spec::GameSpec) = (2, (spec.M.second[1] * spec.N.second[1] + 2))#Opperations + source and sink

function GI.set_state!(g::GameEnv, s)
  g.process_time = s.process_time
  g.M = s.M
  g.N = s.N
  g.N_OPP = s.N_OPP
  g.T = s.T
  g.S = s.S
  g.UB = s.UB
  g.LB = s.LB
  g.conj_src = s.conj_src
  g.conj_tar = s.conj_tar
  g.action_src = copy(s.action_src)
  g.action_tar = copy(s.action_tar)
  g.action_info = copy(s.action_info)
  g.disj_src = copy(s.disj_src)
  g.disj_tar = copy(s.disj_tar)
  g.is_done = copy(s.is_done)
  g.done_time = copy(s.done_time)
  g.action_done_time = copy(s.action_done_time)
  g.prev_operation = copy(s.prev_operation)
  g.prev_machine = copy(s.prev_machine)
end

GI.current_state(g::GameEnv) = (
  process_time=g.process_time,
  M=g.M,
  N=g.N,
  N_OPP=g.N_OPP,
  T=g.T,
  S=g.S,
  UB=g.UB,
  LB=g.LB,
  conj_src=g.conj_src,
  conj_tar=g.conj_tar,
  action_src=copy(g.action_src),
  action_tar=copy(g.action_tar),
  action_info=copy(g.action_info),
  disj_src=copy(g.disj_src),
  disj_tar=copy(g.disj_tar),
  is_done=copy(g.is_done),
  done_time=copy(g.done_time),
  action_done_time=copy(g.action_done_time),
  prev_operation=copy(g.prev_operation),
  prev_machine=copy(g.prev_machine),
)

GI.white_playing(g::GameEnv) = true

GI.game_terminated(g::GameEnv) = g.is_done[g.T]

GI.actions(spec::GameSpec) = collect(1:(spec.M.second[1]*spec.N.second[1]+2))

GI.available_actions(g::GameEnv) = collect(g.S+1:g.S+length(g.action_done_time)+1)

GI.actions_mask(g::GameEnv) = [falses(g.S); trues(length(g.action_done_time))]

function GI.play!(g::GameEnv, action)
  #mark operation scheduled
  m, o = g.action_info[action-g.S, :]
  @assert g.is_done[o] == false
  g.is_done[o] = true
  i = o2ji(o, g.M, g.N)[2]
  #update previous operation and machine
  k = g.prev_machine[m] #previous operation done on machine of todo operation
  l = g.prev_operation[i] #previous operation done in job of todo operation
  #update info vectors
  g.prev_machine[m] = o
  g.prev_operation[i] = o
  #add disjunctive link
  append!(g.disj_src, k)
  append!(g.disj_tar, o)
  #convert from edge to node notation
  l = min(l, g.S)
  k = min(k, g.S)
  g.done_time[o] = max(g.done_time[l], g.done_time[k] * g.is_done[k]) + get(g.process_time[o], m, nothing)
  #remove old actions
  mask = g.action_info[:, 2] .!== o
  g.action_src = g.action_src[mask, :]
  g.action_tar = g.action_tar[mask, :]
  g.action_info = g.action_info[mask, :]
  g.action_done_time = g.action_done_time[mask]
  #fix disjunctive edges
  g.action_src[g.action_info[:, 1] .== m, 3] .= o
  #fix self reference
  g.action_src[:, 2] = collect(g.S+1:g.S+length(g.action_done_time))
  g.action_tar[:, 1] = g.action_src[:, 2]
  g.action_tar[:, 3] = g.action_src[:, 2]
  #add new actions
  next_o = g.conj_tar[o]
  next_next_o = g.conj_tar[next_o]
  for (m, p_time) in get(g.process_time, next_o, [])
    self = g.S + length(g.action_done_time) + 1
    g.action_src = [g.action_src; [o self min(g.prev_machine[m], g.S)]]
    g.action_tar = [g.action_tar; [self next_next_o self]]
    g.action_info = [g.action_info; [m next_o]]
    append!(g.action_done_time, g.done_time[o] + p_time)
  end
  #mark last operation as done
  isempty(g.action_done_time) && (g.is_done[g.T] = true)
  #propagate expected done time
  last_done_time = g.done_time[o]
  while (true)
    o = g.conj_tar[o]
    if (o == g.T)
      g.done_time[g.T] = max(g.done_time[g.T], last_done_time)
      return
    end
    last_done_time = g.done_time[o] = last_done_time + minimum(values(g.process_time[o]))
  end
end

function GI.white_reward(g::GameEnv)
  if (all(g.is_done))
    #return -convert(Float32, g.done_time[g.T])
    return ((g.UB - g.done_time[g.T]) / (g.UB - g.LB))
  end
  return 0
end

function GI.vectorize_state(::GameSpec, state)
  num_actions = length(state.action_done_time)
  return GNNGraph([state.conj_src; state.disj_src; vec(state.action_src)],
    [state.conj_tar; state.disj_tar; vec(state.action_tar)],
    num_nodes = state.S + num_actions,
    ndata=Float32.([[state.done_time; state.action_done_time] [state.is_done; zeros(num_actions)]]'))
end

#####
##### Interaction APIc
#####

function GI.action_string(spec::GameSpec, o)
  mn = o2ji(o, spec.M.second, spec.N.second)
  return string("job: ", mn[2], " machine: ", mn[1])
end

function GI.parse_action(spec::GameSpec, str)
  try
    s = split(str)
    @assert length(s) == 2
    n = parse(Int, s[1])
    m = parse(Int, s[2])
    return ji2o(m, n, spec.M.second, spec.N.second)
  catch e
    return nothing
  end
end

function GI.read_state(::GameSpec)#input problem
  return nothing
end

function GI.render(g::GameEnv)
  for (m, o) in enumerate(g.disj_tar[g.S:end])
    print(crayon"white", "m", m, ":")
    o == g.S && (println(); continue)
    last_op = 0
    while (true)
      mn = o2ji(o, g.M, g.N)
      idle = (crayon"dark_gray", repeat("-", g.done_time[o] - (last_op + g.process_time[o])))
      active = (Crayon(foreground=mn[2]), repeat("=", g.process_time[o]))
      print(idle..., active...)
      last_op = g.done_time[o]
      g.disj_tar[o] == o && break
      o = g.disj_tar[o]
    end
    println(crayon"reset")
  end
  for n in 1:g.N
    print(Crayon(foreground=n), "n", n, " ")
  end
  println(crayon"reset")
end


text:
using Test, Distributions, HypothesisTests, LinearAlgebra
using KernelGoodnessOfFit

xs = randn(1, 100);

@testset "FSSDrand" begin
    xs = randn(100);

    res = KernelGoodnessOfFit.FSSDrand(
        reshape(xs, 1, :),
        MvNormal([1.0], [1.0]),
        GaussianRBF(1.0),
        reshape([0.0], 1, 1),
        nsim = 3000
    )

    println(res)

    # reject
    @test pvalue(res) ≤ 0.05
end

@testset "FSSDopt" begin
    xs = randn(100);

    res = KernelGoodnessOfFit.FSSDopt(
        reshape(xs, 1, :),
        MvNormal([1.0], [1.0]),
        GaussianRBF(1.0),
        reshape([0.0], 1, 1),
        nsim = 3000,
        train_test_ratio = 0.5
    )

    println(res)

    # reject
    @test pvalue(res) ≤ 0.05
end

@testset "FSSDrand multivariate" begin
    d = 5
    xs = randn(d, 100);

    res = KernelGoodnessOfFit.FSSDrand(
        xs,
        MvNormal(ones(d), diagm(0 => ones(d))),
        GaussianRBF(1.0),
        zeros(5, 2);  # two test-locations
        nsim = 3000
    )

    println(res)

    # reject
    @test pvalue(res) ≤ 0.05
end

@testset "FSSDopt multivariate" begin
    d = 5
    xs = randn(d, 100);

    res = pvalue(KernelGoodnessOfFit.FSSDopt(
        xs,
        MvNormal(10 .* ones(d), diagm(0 => ones(d))),
        GaussianRBF(1.0),
        zeros(5, 2);  # two test-locations
        nsim = 3000,
        train_test_ratio = 0.5
    ))

    # reject
    @test res ≤ 0.05
end

@testset "FSSDrand different kernels" begin
    q = MvNormal([10.0], [1.0])

    t_rbf = FSSDrand(xs, q, KernelGoodnessOfFit.GaussianRBF(1.0))
    t_exp = FSSDrand(xs, q, KernelGoodnessOfFit.ExponentialKernel())
    t_matern = FSSDrand(xs, q, KernelGoodnessOfFit.Matern25Kernel(1.0))

    @test pvalue(t_rbf) ≥ 0.0
    @test pvalue(t_exp) ≥ 0.0
    @test pvalue(t_matern) ≥ 0.0
end

@testset "FSSDopt different kernels" begin
    q = MvNormal([10.0], [1.0])

    t_rbf = FSSDopt(xs, q, KernelGoodnessOfFit.GaussianRBF(1.0))
    t_exp = FSSDopt(xs, q, KernelGoodnessOfFit.ExponentialKernel())
    # t_matern = FSSDrand(xs, q, KernelGoodnessOfFit.Matern25Kernel(1.0))  # OPTIMIZATION NOT STABLE

    @test pvalue(t_rbf) ≥ 0.0
    @test pvalue(t_exp) ≥ 0.0
    # @test pvalue(t_matern) ≥ 0.0
end


text:
"""
    swap_our_param_with_jump_param(jm, expr)

This function takes an expression tree and replaces all instances of
`ParameterRef` with the corresponding `JuMP.NLParameter`.
"""
function swap_our_param_with_jump_param(jm, expr)
    return MacroTools.postwalk(expr) do x
        if x isa ParameterRef
            return jm[x.model._parameters[x.index].name]
        else
            return x
        end
    end
end

"""
swap_our_param_with_val(expr)

This function takes an expression tree and replaces all instances of
`ParameterRef` with its value.
"""
function swap_our_param_with_val(expr)
    return MacroTools.postwalk(expr) do x
        if x isa ParameterRef
            return x.model._parameters[x.index].value
        else
            return x
        end
    end
end

function get_jump_variable_for_sector(jm, sector)
    if sector.subindex===nothing
        return jm[get_name(sector)]
    else
        return jm[get_name(sector)][sector.subindex]
    end
end

function get_jump_variable_for_commodity(jm, commodity::CommodityRef)
    if commodity.subindex===nothing
        return jm[get_name(commodity)]
    else
        return jm[get_name(commodity)][commodity.subindex]
    end
end

function get_jump_variable_for_commodity(jm, commodity::ScalarCommodity)
    return jm[get_name(commodity)]
end

function get_jump_variable_for_commodity(jm, commodity::IndexedCommodity)
    return jm[get_name(commodity)][commodity.subindex]
end

function get_jump_variable_for_consumer(jm, consumer)
    return jm[get_name(consumer)]
end

function get_jump_variable_for_intermediate_supply(jm, output)
    return jm[get_comp_supply_name(output)]
end

function get_jump_variable_for_intermediate_demand(jm, input)
    return jm[get_comp_demand_name(input)]
end

function get_jump_variable_for_final_demand(jm, demand)
    return jm[get_final_demand_name(demand)]
end

function get_prod_func_name(x::Production)
    return Symbol("$(get_name(x.sector, true))")
end

function get_demand_func_name(x::DemandFunction)
    return Symbol("$(get_name(x.consumer))")
end

function get_comp_demand_name(i::Input)
    p = i.production_function::Production
    return Symbol("$(get_name(i.commodity, true))†$(get_prod_func_name(p))")
end

function get_comp_supply_name(o::Output)
    p = o.production_function::Production
    return Symbol("$(get_name(o.commodity, true))‡$(get_prod_func_name(p))")
end

function get_final_demand_name(demand::Demand)
    demand_function = demand.demand_function::DemandFunction
    return Symbol("$(get_name(demand.commodity, true))ρ$(get_demand_func_name(demand_function))")
end


text:
"""
  DataLoader(dataset, batchsize::Int; shuffle::Bool=true, droplast::Bool=false)

Interatable minibatch loader.

`dataset` should satisfy `IndexStyle(dataset) == IndexLinear()`. `getindex(dataset, index::Int)` and `size(dataset)` is required too.

If `droplast` is `true`, last batch whoes size is less than `batchsize` is not iterated.
"""
struct DataLoader
  dataset
  batchsize::Int
  ntasks::Int
  shuffle::Bool
  droplast::Bool

  function DataLoader(::IndexLinear, dataset, batchsize::Int; ntasks::Int = 1, shuffle::Bool = true, droplast::Bool = false)
    if shuffle
      new(RandomSampler(dataset), batchsize, ntasks, shuffle, droplast)
    else
      new(dataset, batchsize, ntasks, shuffle, droplast)
    end
  end
end

function DataLoader(dataset, batchsize::Int; ntasks::Int = 1, shuffle::Bool = true, droplast::Bool = false)
  DataLoader(IndexStyle(dataset), dataset, batchsize, ntasks = ntasks, shuffle = shuffle, droplast = droplast)
end

Base.size(dl::DataLoader) = (length(dl),)

Base.length(dl::DataLoader) = round(Int, length(dl.dataset) / dl.batchsize, dl.droplast ? RoundUp : RoundDown)

function getbatch(dl::DataLoader, first::Int, last::Int)
  xdim = length(size(dl.dataset[1][1]))  # dimension of data   (e.g. 2 for image)
  ydim = length(size(dl.dataset[1][2]))  # demension of target (e.g. 1 for lable)

  mapper(x) = @inbounds dl.dataset[x]

  batch = asyncmap(mapper, first:last, ntasks = dl.ntasks)

  xs = cat(getindex.(batch, 1)..., dims=xdim+1)
  ys = cat(getindex.(batch, 2)..., dims=ydim+1)

  return xs, ys
end

function Base.iterate(dl::DataLoader)
  dl.shuffle && Random.shuffle!(dl.dataset)
  return iterate(dl, 1)
end

function Base.iterate(dl::DataLoader, i::Int)
  if length(dl.dataset) < i
    # i is out of range
    return nothing
  elseif length(dl.dataset) < i+dl.batchsize-1
    # the remainings is too few
    if !(dl.droplast)
      return getbatch(dl, i, length(dl)), length(dl.dataset)+1
    else
      return nothing
    end
  end

  return getbatch(dl, i, i+dl.batchsize-1), i+dl.batchsize
end


text:
using Test, DayCounts, Dates
using XLSX
using Tables

xf = XLSX.readxlsx("yearfrac.xlsx")

function range_to_table(rng)
    Tables.table(rng[2:end,:]; header = Symbol.(rng[1,:]))
end

@testset "Excel basis 0" begin
    table = range_to_table(xf["basis0"][:])
    dc = DayCounts.Thirty360Excel()
    for r in Tables.rows(table)
        @test yearfrac(Date(r.startdate), Date(r.enddate), dc) == r.yearfrac
    end
end

@testset "Excel basis 1" begin
    table = range_to_table(xf["basis1"][:])
    dc = DayCounts.ActualActualExcel()
    for r in Tables.rows(table)
        @test yearfrac(Date(r.startdate), Date(r.enddate), dc) == r.yearfrac
    end
end

@testset "Excel basis 2" begin
    table = range_to_table(xf["basis2"][:])
    dc = DayCounts.Actual360()
    for r in Tables.rows(table)
        @test yearfrac(Date(r.startdate), Date(r.enddate), dc) == r.yearfrac
    end
end

@testset "Excel basis 3" begin
    table = range_to_table(xf["basis3"][:])
    dc = DayCounts.Actual365Fixed()
    for r in Tables.rows(table)
        @test yearfrac(Date(r.startdate), Date(r.enddate), dc) == r.yearfrac
    end
end

@testset "Excel basis 4" begin
    table = range_to_table(xf["basis4"][:])
    dc = DayCounts.ThirtyE360()
    for r in Tables.rows(table)
        @test yearfrac(Date(r.startdate), Date(r.enddate), dc) == r.yearfrac
    end
end


text:
function write_xyz(mol::Molecule, filename)
    lines = String[]
    push!(lines, string(length(mol.atoms)))
    push!(lines, "Properties=species:S:1:pos:R:3")
    for atom in mol.atoms
        push!(lines, join([chemical_symbol(atom), position(atom)...], "\t"))
    end
    Base.write(filename, join(lines, "\n"))
end

function read_xyz(filename)
    lines = readlines(filename)
    N = parse(Int, lines[1])
    properties = split(lines[2])
    atoms = Atom[]
    for (idx, line) in enumerate(lines[3:end])
        args = split(line)
        sym = Symbol(args[1])
        pos = map(x -> parse(Float64, x), args[2:4])
        push!(atoms, chemical_symbol_to_atom_type(sym)(pos))
    end
    @assert N == length(atoms)
    return N, atoms
end

function read_xyz!(mol::Molecule, filename)
    N, atoms = read_xyz(filename)
    @assert length(mol) == N
    for (atom1, atom2) in zip(mol.atoms, atoms)
        @assert chemical_symbol(atom1) == chemical_symbol(atom2)
        atom1.position = atom2.position
    end
    return mol
end


text:
# These do most of the work in the package.
# They are all @generated or recusive functions for performance.

const UnionAllTupleOrVector = Union{Vector{UnionAll},Tuple{UnionAll,Vararg}}

"""
Sort dimensions into the order they take in the array.

Missing dimensions are replaced with `nothing`
"""
@inline Base.permutedims(tosort::AbDimTuple, perm::Union{Vector{<:Integer},Tuple{<:Integer,Vararg}}) =
    map(p -> tosort[p], Tuple(perm))
@inline Base.permutedims(tosort::Union{Vector{<:Integer},Tuple{<:Integer,Vararg}}, perm::AbDimTuple) =
    tosort

@inline Base.permutedims(tosort::AbDimTuple, order::UnionAllTupleOrVector ) =
    permutedims(tosort, Tuple(map(u -> u(), order)))
@inline Base.permutedims(tosort::UnionAllTupleOrVector, order::AbDimTuple) =
    permutedims(Tuple(map(u -> u(), tosort)), order)
@inline Base.permutedims(tosort::AbDimTuple, order::AbDimVector) =
    permutedims(tosort, Tuple(order))
@inline Base.permutedims(tosort::AbDimVector, order::AbDimTuple) =
    permutedims(Tuple(tosort), order)

@generated Base.permutedims(tosort::AbDimTuple, order::AbDimTuple) =
    permutedims_inner(tosort, order)
permutedims_inner(tosort::Type, order::Type) = begin
    indexexps = []
    for dim in order.parameters
        index = findfirst(d -> basetype(d) <: basetype(dim), tosort.parameters)
        if index == nothing
            push!(indexexps, :(nothing))
        else
            push!(indexexps, :(tosort[$index]))
        end
    end
    Expr(:tuple, indexexps...)
end

"""
Convert a tuple of AbstractDimension to indices, ranges or Colon.
"""
@inline dims2indices(a, lookup, emptyval=Colon()) =
    dims2indices(dims(a), lookup, emptyval)
@inline dims2indices(dims::Tuple, lookup, emptyval=Colon()) =
    dims2indices(dims, (lookup,), emptyval)
@inline dims2indices(dims::Tuple, lookup::Tuple, emptyval=Colon()) =
    _dims2indices(dims, permutedims(lookup, dims), emptyval)

@inline _dims2indices(dims::Tuple, lookup::Tuple, emptyval) =
    (dims2indices(dims[1], lookup[1], emptyval),
     _dims2indices(tail(dims), tail(lookup), emptyval)...)
@inline _dims2indices(dims::Tuple, lookup::Tuple{}, emptyval) =
    (emptyval, _dims2indices(tail(dims), (), emptyval)...)
@inline _dims2indices(dims::Tuple{}, lookup::Tuple{}, emptyval) = ()

@inline dims2indices(dim::AbDim, lookup::AbDim, emptyval) = val(lookup)
@inline dims2indices(dim::AbDim, lookup::Type{<:AbDim}, emptyval) = Colon()
@inline dims2indices(dim::AbDim, lookup::Nothing, emptyval) = emptyval

"""
Slice the dimensions to match the axis values of the new array

All methods returns a tuple conatining two tuples: the new dimensions,
and the reference dimensions. The ref dimensions are no longer used in
the new struct but are useful to give context to plots.

Called at the array level the returned tuple will also include the
previous reference dims attached to the array.
"""
@inline slicedims(a, dims::AbDimTuple) =
    slicedims(a, dims2indices(a, dims))
@inline slicedims(a, I::Tuple) = slicedims(dims(a), refdims(a), I)
@inline slicedims(dims::Tuple, refdims::Tuple, I::Tuple) = begin
    newdims, newrefdims = slicedims(dims, I)
    newdims, (refdims..., newrefdims...)
end
@inline slicedims(dims::Tuple{}, I::Tuple) = ((), ())
@inline slicedims(dims::AbDimTuple, I::Tuple) = begin
    d = slicedims(dims[1], I[1])
    ds = slicedims(tail(dims), tail(I))
    (d[1]..., ds[1]...), (d[2]..., ds[2]...)
end
@inline slicedims(dims::Tuple{}, I::Tuple{}) = ((), ())
@inline slicedims(d::AbDim, i::Colon) = ((rebuild(d, val(d)),), ())
@inline slicedims(d::AbDim, i::Number) = ((), (rebuild(d, val(d)[i]),))
@inline slicedims(d::AbDim{<:AbstractArray}, i::AbstractArray) =
    ((rebuild(d, val(d)[i]),), ())
@inline slicedims(d::AbDim{<:LinRange}, i::UnitRange) = begin
    range = val(d)
    start, stop, len = range[first(i)], range[last(i)], length(i) ÷ step(i)
    ((rebuild(d, LinRange(start, stop, len)),), ())
end

"""
Get the number of an AbstractDimension as ordered in the array
"""
@inline dimnum(a, lookup) = dimnum(typeof(dims(a)), lookup)
@inline dimnum(dimtypes::Type, lookup::AbstractArray) = dimnum(dimtypes, (lookup...,))
@inline dimnum(dimtypes::Type, lookup::Number) = lookup
@inline dimnum(dimtypes::Type, lookup::Tuple) =
    (dimnum(dimtypes, lookup[1]), dimnum(dimtypes, tail(lookup))...,)
@inline dimnum(dimtypes::Type, lookup::Tuple{}) = ()
@inline dimnum(dimtypes::Type, lookup::AbDim) = dimnum(dimtypes, typeof(lookup))
@generated dimnum(dimtypes::Type{DTS}, lookup::Type{D}) where {DTS,D} = begin
    index = findfirst(dt -> D <: basetype(dt), DTS.parameters)
    if index == nothing
        :(throw(ArgumentError("No $lookup in $dimtypes")))
    else
        :($index)
    end
end

"""
Format the dimension to match internal standards.

Mostily this means converting tuples and UnitRanges to LinRange,
which is easier to handle.

Errors are thrown if dims don't match the array dims or size.
"""
@inline formatdims(a::AbstractArray{T,N}, dims::Tuple) where {T,N} = begin
    dimlen = length(dims)
    dimlen == N || throw(ArgumentError("dims ($dimlen) don't match array dimensions $(N)"))
    formatdims(size(a), dims)
end
@inline formatdims(size::Tuple, dims::Tuple) =
    (formatdims(size[1], dims[1]), formatdims(tail(size), tail(dims))...,)
@inline formatdims(size::Tuple{}, dims::Tuple{}) = ()
@inline formatdims(len::Integer, dim::AbDim{<:AbstractArray}) =
    if length(val(dim)) == len
        dim
    else
        throw(ArgumentError("length of $dim $(length(val(dim))) does not match
                             size of array dimension $len"))
    end
@inline formatdims(len::Integer, dim::AbDim{<:Union{UnitRange,NTuple{2}}}) = linrange(dim, len)

linrange(dim, len) = begin
    range = val(dim)
    start, stop = first(range), last(range)
    rebuild(dim, LinRange(start, stop, len))
end

"""
Replace the specified dimensions with an index of 1 to match
a new array size where the dimension has been reduced to a length
of 1, but the number of dimensions has not changed.

Used in mean, reduce, etc.
"""
@inline reducedims(dim) = reducedims((dim,))
@inline reducedims(dims::Tuple) = map(d -> basetype(d)(1), dims)

"""
Get the dimension(s) matching the type(s) of the lookup dimension.
"""
@inline dims(a::AbstractArray, lookup) = dims(dims(a), lookup)
@inline dims(ds::AbDimTuple, lookup::Integer) = ds[lookup]
@inline dims(ds::AbDimTuple, lookup::Tuple) =
    (dims(ds, lookup[1]), dims(ds, tail(lookup))...)
@inline dims(ds::AbDimTuple, lookup::Tuple{}) = ()
@inline dims(ds::AbDimTuple, lookup) = dims(ds, basetype(lookup))
@generated dims(ds::DT, lookup::Type{L}) where {DT<:AbDimTuple,L} = begin
    index = findfirst(dt -> dt <: L, DT.parameters)
    if index == nothing
        :(throw(ArgumentError("No $lookup in $dims")))
    else
        :(ds[$index])
    end
end


text:
module ProbabilisticExprRules

using ExprRules
using StatsBase
using AbstractTrees
using LinearAlgebra

export
    @grammar,
    RuleNode,

    ProbabilisticGrammar,
    probabilities,

    fit_mle!,
    uniform!

const ProbType = Dict{Symbol,Vector{Float64}}

struct ProbabilisticGrammar
    grammar::Grammar
    probs::ProbType
end
function ProbabilisticGrammar(grammar::Grammar)
    probs = ProbType()
    for nt in nonterminals(grammar)
        n = length(grammar[nt])
        probs[nt] = ones(Float64, n) / n
    end
    ProbabilisticGrammar(grammar, probs)
end

probabilities(pcfg::ProbabilisticGrammar, typ::Symbol) = pcfg.probs[typ]

"""
    rand(::Type{RuleNode}, grammar::ProbabilisticGrammar, typ::Symbol, dmap::AbstractVector{Int},
        max_depth::Int=10)

Generates a random RuleNode of return type typ and maximum depth max_depth.
"""
function Base.rand(::Type{RuleNode}, pcfg::ProbabilisticGrammar, typ::Symbol, dmap::AbstractVector{Int},
                    max_depth::Int=10)
    grammar = pcfg.grammar
    rules = grammar[typ]
    probs = probabilities(pcfg, typ)

    inds = findall(r->dmap[r] ≤ max_depth, rules)
    rules, probs = rules[inds], probs[inds]
    rule_index = StatsBase.sample(rules, weights(probs))

    rulenode = grammar.iseval[rule_index] ?
        RuleNode(rule_index, Core.eval(Main, grammar.rules[rule_index].args[2])) :
        RuleNode(rule_index)

    if !grammar.isterminal[rule_index]
        for ch in child_types(grammar, rule_index)
            push!(rulenode.children, rand(RuleNode, pcfg, ch, dmap, max_depth-1))
        end
    end
    return rulenode
end

"""
    uniform!(pcfg::ProbabilisticGrammar)

Set all probability vectors to uniform distribution
"""
function uniform!(pcfg::ProbabilisticGrammar)
    fill!(pcfg, 0.0)
    normalize!(pcfg)
end

"""
    fit_mle!(pcfg::ProbabilisticGrammar, pop::AbstractVector{RuleNode}; initial_value::Float64=0.0)

Update the probability vectors based on population using MLE.
"""
function fit_mle!(pcfg::ProbabilisticGrammar, pop::AbstractVector{RuleNode}, initial_value::Float64=0.0)
    fill!(pcfg, initial_value)
    for x in pop
        _fit_mle!(pcfg, x)
    end
    normalize!(pcfg)
    pcfg
end
function _fit_mle!(pcfg::ProbabilisticGrammar, x::RuleNode)
    grammar = pcfg.grammar
    typ = return_type(grammar, x)
    i = something(findfirst(isequal(x.ind),grammar[typ]), 0)
    pcfg.probs[typ][i] += 1.0
    for c in x.children
        _fit_mle!(pcfg, c)
    end
    pcfg
end

"""
    Base.fill!(pcfg::ProbabilisticGrammar, x::Float64)

Fill all probability values to x
"""
function Base.fill!(pcfg::ProbabilisticGrammar, x::Float64)
    for v in values(pcfg.probs)
        fill!(v, x)
    end
end

"""
    LinearAlgebra.normalize!(pcfg::ProbabilisticGrammar)

Normalize all probability vectors so that each vector sums to 1.0
"""
function LinearAlgebra.normalize!(pcfg::ProbabilisticGrammar)
    for v in values(pcfg.probs)
        if sum(v) > 0.0
            normalize!(v, 1)
        else
            fill!(v, 1.0/length(v))
        end
    end
end

end # module


text:
using IncrementalInference, KernelDensityEstimate, Distributions
using Gadfly # for draw PDF
using Base.Test

import IncrementalInference: getSample

# switch off y ticks
toggleYTicks()

type MultiModalConditional <: IncrementalInference.FunctorPairwise
  x::Vector{Distribution}
  hypo::Categorical
  MultiModalConditional{D <: Distribution}(x::Vector{D}, p::Categorical) = new(x, p)
end
function getSample(dpl::MultiModalConditional, N::Int=1)
  d = length(dpl.hypo.p)
  p = rand(dpl.hypo, N)
  ret = zeros(1,N)
  for i in 1:N
    ret[i] = rand(dpl.x[p[i]])
  end
  return (ret, p)
end

function (dp::MultiModalConditional)(res::Vector{Float64},
                                    idx::Int,
                                    meas::Tuple{Array{Float64,2},Vector{Int64}},
                                    x1::Array{Float64},
                                    x2::Array{Float64}  )
  #
  res[1] = meas[1][1,idx] - (x2[1,idx]-x1[1,idx])
  nothing
end

fg = emptyFactorGraph()

N=100

doors = [-20.0, 0.0, 20.0]'
pd = kde!(doors,[2.0])
pd = resample(pd,N);
bws = getBW(pd)[:,1]
doors2 = getPoints(pd);
v1 = addNode!(fg,:x1,doors,N=N)
f1  = addFactor!(fg,[v1],Obsv2( doors2, bws', [1.0])) #, samplefnc=getSample

# not initialized
v2 = addNode!(fg,:x2, N=N)

mmc = MultiModalConditional([Normal(-5,0.5),Normal(5,0.5)],Categorical([0.5,0.5]))
f2 = addFactor!(fg, [:x1; :x2], mmc )

# Graphs.plot(fg.g)

pts = approxConv(fg, :x1x2f1, :x2)

q2 = kde!(getSample(mmc,2000)[1])
h1 = plotKDE([getVertKDE(v1), q2],c=["red";"green"],fill=true, xlbl="")
h2 = plotKDE(kde!(pts),fill=true,xlbl="", title="N = 100")

draw(PDF("approxconv.pdf",14cm,10cm),vstack(h1,h2))
# @async run(`evince approxconv.pdf`)

#


text:
@time @testset "[compat] entries" begin
    compat_regex_1 = r"^[\d]*?$"
    compat_regex_2 = r"^[\d]*?\.[\d]*?$"
    compat_regex_3 = r"^[\d]*?\.[\d]*?\.[\d]*?$"

    for (pkgname, pkguuid) in predictmd_direct_deps
        if !Pkg.Types.is_stdlib(Base.UUID(pkguuid))
            @testset "Dependency: $(pkgname)" begin
                @test haskey(predictmd_project_compat, pkgname)
                pkgcompat = predictmd_project_compat[pkgname]
                matches_regex_1 = occursin(compat_regex_1, pkgcompat)
                matches_regex_2 = occursin(compat_regex_2, pkgcompat)
                matches_regex_3 = occursin(compat_regex_3, pkgcompat)
                @test matches_regex_1 || matches_regex_2 || matches_regex_3
                ver = Base.VersionNumber(pkgcompat)
                @test isempty(ver.prerelease)
                @test isempty(ver.build)
                @test ver > v"0"
                @test ver >= v"0.1"
                if ver.major < 1
                    @test ver.minor >= 1
                end
                if matches_regex_1
                    major = parse(Int, pkgcompat)
                    @test major >= 1
                    @test ver >= v"1"
                end
            end
        end
    end
end


text:
using Base.Pkg

try Pkg.clone("https://github.com/camshaft/PatternDispatch.jl") end
Pkg.checkout("PatternDispatch", "new_pattern")


text:
export Dense

"""
    Dense([T::Type=ComplexF32], in::Integer, out::Integer, σ = identity)

Creates a traditional `Dense` layer with parameters `W` and `b`.

    y = σ.(W * x .+ b)

The input `x` must be a vector of length `in`, or a batch of vectors represente
d
as an `in × N` matrix. The out `y` will be a vector or batch of length `out`.

```julia
julia> d = Dense(5, 2)
Dense(5, 2)
julia> d(rand(5))
2-element Array{Float64,1}:
  0.00257447
  -0.00449443
```

"""
struct Dense{Ta,Tb,C}
    W::Ta
    b::Tb
    σ::C
end
functor(d::Dense) = (W=d.W, b=d.b), t -> Dense(t.W,t.b,d.σ)
(l::Dense)(x) = l.σ.(l.W*x .+ l.b)

Dense(in::Integer, args...;kwargs...)= Dense(Complex{STD_REAL_PREC}, in, args...;kwargs...)
function Dense(T::Type, in::Integer, out::Integer, σ = identity;
               initW = glorot_uniform, initb = glorot_uniform)
  return Dense(initW(T, out, in), initb(T, out), σ)
end

function Base.show(io::IO, l::Dense)
  print(io, "Dense(", size(l.W, 2), ", ", size(l.W, 1))
  l.σ == identity || print(io, ", ", l.σ)
  print(io, ")")
end

struct DenseCache{Ta,Tb,Tc,Td,Te}
    σ::Tc
    out::Tb
    δℒℒ::Td
    δℒ::Te

    θ::Ta
    valid::Bool
end

function cache(l::Dense{Ta,Tb}, arr_T, in_T, in_sz) where {Ta,Tb}
    in_T2 = promote_type(in_T, eltype(l.W))
    c = DenseCache(similar(l.W, in_T2, size(l.W,2)),
                   similar(l.b, in_T2),
                   similar(l.W, size(l.W,1)),
                   similar(l.W, 1,   size(l.W,2)),

                   similar(l.b, in_T2),
                   false)
    return c
end

function layer_out_type_size(l::Dense, in_T ,in_sz)
    T1     = promote_type(in_T, eltype(l.W))
    out_T  = promote_type(T1, eltype(l.b))
    out_sz = size(l.b)
    return out_T, out_sz
end

function (l::Dense)(c::DenseCache, x)
    # The preallocated caches
    logℒθ  = c.out
    θ = c.θ

    # Store the input to this layer for the backpropagation
    σ = copyto!(c.σ, x)

    #θ .= net.b .+ net.W * x
    mul!(θ, l.W, σ)
    θ .+= l.b

    # Apply the nonlinear function
    logℒθ .= l.σ.(θ)
    return logℒθ
end

function backprop(∇, l::Dense, c::DenseCache, δℒ)
    # The preallocated caches
    θ = c.θ
    δℒℒ = c.δℒℒ

    # Compute the actual sensitivity
    copyto!(δℒℒ, δℒ)
    δℒℒ .*= fwd_der.(l.σ, θ)

    ∇.W .= δℒℒ.*transpose(c.σ)
    ∇.b .= δℒℒ

    δℒ = mul!(c.δℒ, transpose(δℒℒ), l.W)
    return δℒ
end


text:
"""
    _xi(r)

Returns the value of 'xi' given by the change of variable `r->xi`.
See `Aoki et al. (1979)`.
"""
function _xi(r::Float64)
    r_a = r/a
    return (r_a^2-1)/(r_a^2+1)
end

"""
    radius(xi)

Returns the value of the radius given the value of `xi`.
See `Aoki et al. (1979)`.
"""
function radius(xi::Float64)
    return a*sqrt((1+xi)/(1-xi))
end

"""
    Omega(xi,x,q)

Computes the dimensionless angular velocity of the system.
See `Reddish et al. (2021)`.
"""
function Omega(xi::Float64,x::Float64=0.0, q::Float64=1.0)
    r = radius(xi)
    bulb = (a/c)^3*x/(1-x)*((1-xi)/2)^(-3/2)*(1+(r/c)^2)^(-3/2)
    halo = (a/b)^3*(1/q-1)*((1-xi)/2)^(-3/2)*(1+(r/b)^2)^(-3/2)
    disk = 1-eps/((1-x)^(2/3))
    return sqrt(1-x)*((1-xi)/2)^(3/4)*sqrt(bulb+halo+disk)
end

"""
    dOmegadxi(xi,x,q)

Computes the `xi`-gradient of the dimensionless angular velocity of the system.
See `Aoki et al. (1979)` and `Reddish et al. (2021)`.
"""
function dOmegadxi(xi::Float64,x::Float64=0.0, q::Float64=1.0)
    rc = a/c
    rb = a/b
    num = (-3+3*eps/(1-x)^(2/3)+(12*(-1+q)*rb^5*sqrt((2-2*xi)/(1-xi+rb^2*(1+xi))))/(q*sqrt(1-xi)*(1-xi+rb^2*(1+xi))^2)
           +(12*rc^5*x*sqrt((2-2*xi)/(1-xi+rc^2*(1+xi)))/((-1+x)*sqrt(1-xi)*(1-xi+rc^2*(1+xi))^2)))
    den = (4*2^(3/4)*(1-xi)^(1/4)*sqrt(1-eps/(1-x)^(2/3)+(2*sqrt(2)*(-1+1/q)*rb^3*(1-(rb^2*(1+xi))/(-1+xi))^(-3/2))
            /(1-xi)^(3/2)+(2*sqrt(2)*x*rc^3*(1-(rc^2*(1+xi))/(-1+xi))^(-3/2))
                    /((1-x)*(1-xi)^(3/2))))

    return sqrt(1-x)*(num/den)
end

"""
    kappaSqOverTwoOmega(xi,x,q)

Computes the dimensionless ratio of the squared epicyclic frequency with twice
the angular velocity of the system.
See `Aoki et al. (1979)` and `Reddish et al. (2021)`.
"""
function kappaSqOverTwoOmega(xi::Float64,x::Float64=0.0, q::Float64=1.0)
    omega = Omega(xi,x,q)
    domega = dOmegadxi(xi,x,q)
#    return (1-eps)^(1/2)*((1-xi)/2)^(3/4)*(1/2+(3/2)*(1-xi)/2)
    return 2*omega*(1+(1+xi)*(1-xi)*domega/(2*omega))
end

##############################
# Test
##############################

"""
    dOmegadxi(xi,x,q)

Computes the numerical `xi`-gradient of the dimensionless angular velocity of the system.
"""
function dOmegadxiNum(xi::Float64,x::Float64=0.0, q::Float64=1.0,prec::Float64=0.0001)
    return (Omega(xi+prec,x,q)-Omega(xi-prec,x,q))/(2*prec)
end

max_stars_repo_name: KerwannTEP/LiRA

max_issues_repo_name: KerwannTEP/LiRA

max_forks_repo_name: KerwannTEP/LiRA