text:
"""
# module LaminartEqImfilterGPU_FIR

- Julia version: 1.4
- Author: niallcullinane
- Date: 2020-08-20

LAMINART equations from

R. D. S. Raizada and S. Grossberg, ‘Context-sensitive binding by the laminar circuits of V1 and V2: A unified model of perceptual grouping, attention, and orientation contrast’, Visual Cognition, vol. 8, no. 3–5, p. 431, Jun. 2001, doi: 10/d8phhc.

Uses JuliaImages's imfilter with zeros boarding for convolution on GPU with FIR algorithm.
"""
module LaminartEqImfilterGPU_FIR

using NNlib, ImageFiltering, Images, OffsetArrays, ComputationalResources
addresource(CUDALibs)

filter_resource = CUDALibs(Algorithm.FIR())
# retina

function I_u!(r::AbstractArray, I::AbstractArray, p::NamedTuple)
    imfilter!(filter_resource, r, I, p.k_gauss_1)
    @. r = I - r
    return nothing
end

# lgn feedback

function fun_x_lgn!(x_lgn::AbstractArray, x::AbstractArray, p::NamedTuple)
    @. x_lgn = 0.0f0
    # @inbounds begin
    for k ∈ 1:p.K
        @. x_lgn += @view x[:, :, k]
        # end
    end
    return nothing
end

function fun_f!(f_out::AbstractArray, x::AbstractArray, p::NamedTuple)
    @. f_out = (p.μ * x^p.n) / (p.ν^p.n + x^p.n)
    return nothing
end

function func_filter_W!(
    W_out::AbstractArray,
    img::AbstractArray,
    W::AbstractArray,
    W_temp::AbstractArray,
    p::NamedTuple,
)
    #     temp_k = similar(W_out[:, :, 1])
    # @inbounds begin
    for k ∈ 1:p.K
        #     todo fix W
        img_k = @view img[:, :, k]
        out_k = @view W_out[:, :, k]
        imfilter!(filter_resource, out_k, img_k, centered(W[:, :, k, k]))
        for l ∈ 1:p.K
            if l ≠ k
                img_l = @view img[:, :, l]
                imfilter!(
                    filter_resource,
                    W_temp,
                    img_l,
                    centered(W[:, :, k, l]),
                )
                @. out_k += W_temp
            end
        end
    end
    # end
    return nothing
end

# LGN
function fun_dv!(
    dv::AbstractArray,
    v::AbstractArray,
    u::AbstractArray,
    x_lgn::AbstractArray,
    p::NamedTuple,
)
    imfilter!(filter_resource, dv, x_lgn, centered(p.k_gauss_1))
    @. dv =
        p.δ_v * (
            -v + ((1f0 - v) * max(u, 0f0) * (1f0 + p.C_1 * x_lgn)) -
            ((1f0 + v) * p.C_2 * dv)
        )
    return nothing
end

# lgn to l6 and l4
function fun_v_C!(
    v_C::AbstractArray,
    v_p::AbstractArray,
    v_m::AbstractArray,
    v_C_temp1::AbstractArray,
    v_C_temp2::AbstractArray,
    v_C_tempA::AbstractArray,
    p::NamedTuple,
)
    #     V = similar(v_p)
    #     temp = similar(v_p)

    @. v_C_temp2 = exp(-1.0f0 / 8.0f0) * (max(v_p, 0f0) - max(v_m, 0f0))
    imfilter!(
        filter_resource,
        v_C_temp1[:, :, 1],
        v_C_temp2[:, :, 1],
        centered(p.k_gauss_2),
    )

    #     A = similar(v_C)
    #     allocate B to v_C
    # @inbounds begin
    for k ∈ 1:p.K
        a = @view v_C_tempA[:, :, k]
        b = @view v_C[:, :, k]
        imfilter!(
            filter_resource,
            a,
            v_C_temp1[:, :, 1],
            centered(p.k_C_d[:, :, k]),
        )
        imfilter!(
            filter_resource,
            b,
            v_C_temp1[:, :, 1],
            centered(p.k_C_b[:, :, k]),
        )
    end
    # end
    @. v_C =
        p.γ * (max(v_C_tempA - abs(v_C), 0f0) + max(-v_C_tempA - abs(v_C), 0f0))
    return nothing
end

# L6
function fun_dx_v1!(
    dx::AbstractArray,
    x::AbstractArray,
    C::AbstractArray,
    z::AbstractArray,
    x_v2::AbstractArray,
    p::NamedTuple,
)
    @. dx =
        p.δ_c * (
            -x + (
                (1.0f0 - x) * (
                    (p.α * C) + (p.ϕ * max(z - p.Γ, 0f0)) .+ (p.v_21 * x_v2) +
                    p.att
                )
            )
        )
    return nothing
end

#     L4 excit
function fun_dy!(
    dy::AbstractArray,
    y::AbstractArray,
    C::AbstractArray,
    x::AbstractArray,
    m::AbstractArray,
    W_temp::AbstractArray,
    p::NamedTuple,
)
    func_filter_W!(dy, m, p.k_W_p, W_temp, p)
    @. dy = m * dy
    fun_f!(dy, dy, p)
    @. dy = p.δ_c * (-y + ((1f0 - y) * (C + (p.η_p * x))) - ((1f0 + y) * dy))
    return nothing
end

#     l4 inhib
function fun_dm!(
    dm::AbstractArray,
    m::AbstractArray,
    x::AbstractArray,
    W_temp::AbstractArray,
    p::NamedTuple,
)
    func_filter_W!(dm, m, p.k_W_m, W_temp, p)
    fun_f!(dm, dm, p)
    @. dm = p.δ_m * (-m + (p.η_m * x) - (m * dm))
    return nothing
end

#     L2/3 excit
function fun_dz!(
    dz::AbstractArray,
    z::AbstractArray,
    y::AbstractArray,
    H_z::AbstractArray,
    s::AbstractArray,
    p::NamedTuple,
)
    imfilter!(filter_resource, dz, s, centered(p.k_T_p))
    @. dz =
        p.δ_z * (
            -z +
            ((1f0 - z) * ((p.λ * max(y, 0f0)) + H_z + (p.a_23_ex * p.att))) -
            ((z + p.ψ) * dz)
        )
    return nothing
end

#     L2/3 inhib
function fun_ds!(
    ds::AbstractArray,
    s::AbstractArray,
    H_z::AbstractArray,
    p::NamedTuple,
)
    imfilter!(filter_resource, ds, s, centered(p.k_T_m))
    @. ds = p.δ_s * (-s + H_z + (p.a_23_in * p.att) - (s * ds))
    return nothing
end

function fun_H_z!(
    H_z::AbstractArray,
    z::AbstractArray,
    H_z_temp::AbstractArray,
    p::NamedTuple,
)
    #     temp = similar(z)
    @. H_z_temp = max(z - p.Γ, 0f0)
    # @inbounds begin
    for k ∈ 1:p.K
        H_z_k = @view H_z[:, :, k]
        temp_k = @view H_z_temp[:, :, k]
        imfilter!(filter_resource, H_z_k, temp_k, centered(p.k_H[:, :, k]))
    end
    # end
    return nothing
end

#     V2 L6
function fun_dx_v2!(
    dx::AbstractArray,
    x_v2::AbstractArray,
    z_v2::AbstractArray,
    z::AbstractArray,
    p::NamedTuple,
)
    @. dx =
        p.δ_c * (
            -x_v2 + (
                (1f0 - x_v2) * (
                    (p.v12_6 * max(z - p.Γ, 0f0)) +
                    (p.ϕ * max(z_v2 - p.Γ, 0f0)) +
                    p.att
                )
            )
        )
    return nothing
end

# V2 L4 excit
function fun_dy_v2!(
    dy::AbstractArray,
    y_v2::AbstractArray,
    z::AbstractArray,
    x_v2::AbstractArray,
    m_v2::AbstractArray,
    p::NamedTuple,
)
    func_filter_W!(dy, m_v2, p.k_W_p, p)
    fun_f!(dy, dy, p)
    @. dy =
        p.δ_c * (
            -y_v2 + (
                (1.0f0 - y_v2) *
                ((p.v12_4 * max(z - p.Γ, 0f0)) + (p.η_p * x_v2))
            ) - ((1.0f0 + y_v2) * dy)
        )
    return nothing
end

# # # for equilabrium

# # l6 equilabrium
function fun_x_equ!(
    x::AbstractArray,
    C::AbstractArray,
    z::AbstractArray,
    x_v2::AbstractArray,
    x_temp::AbstractArray,
    p::NamedTuple,
)
    # 	@inbounds begin
    @. x_temp = (p.α * C) + (p.ϕ * max(z - p.Γ, 0f0)) + (p.v_21 * x_v2) + p.att
    @. x = x_temp / (1f0 + x_temp)
    # 	end
    return nothing
end

# # l4 excit equilabrium
function fun_y_equ!(
    y::AbstractArray,
    C::AbstractArray,
    x::AbstractArray,
    m::AbstractArray,
    dy_temp::AbstractArray,
    p::NamedTuple,
)
    # 	@inbounds begin
    imfilter!(filter_resource, dy_temp, m, p.k_W_p)
    @. dy_temp = m * dy_temp
    fun_f!(dy_temp, p)
    @. y = (C + (p.η_p * x) - dy_temp) / (1f0 + C + (p.η_p * x) + dy_temp)
    # 	end
    return nothing
end

end


text:
rename(env, ex) = IRTools.prewalk(
  x -> x isa GlobalRef ? getfield(x.mod, x.name) :
  x isa Variable ? env[x] : x, ex)

returntype(ir) = exprtype(ir, returnvalue(IRTools.blocks(ir)[end]))

function copyblock!(ir::IR, b)
  c = IRTools.block!(ir)
  env = Dict()
  for (arg, T) in zip(arguments(b), argtypes(b))
    env[arg] = argument!(c, nothing, T, insert = false)
  end
  for (v, st) in b
    env[v] = push!(c, rename(env, st))
  end
  for br in IRTools.branches(b)
    push!(branches(c), rename(env, br))
  end
  return c
end

function abstract_group(cfg, i)
  j = IRTools.Inner.idoms(cfg', entry = length(cfg))[i]
  group = Int[]
  group!(b) =
    (b == j || b in group) || (push!(group, b); foreach(group!, cfg[b]))
  foreach(group!, cfg[i])
  return sort!(group), j
end

function extract_group(bl)
  out = IR()
  args = filter(x -> x isa Variable, union(arguments.(branches(block(bl.ir, bl.id)))...))
  args = union(args, filter(x -> x != nothing, map(x -> x.condition, branches(block(bl.ir, bl.id)))))
  innerargs = Dict(a => argument!(out, Union{}) for a in args)
  for br in branches(block(bl.ir, bl.id))
    push!(branches(block(out, 1)), rename(innerargs, br))
  end
  bs, after = abstract_group(CFG(bl.ir), bl.id)
  bmap = Dict()
  for b in bs
    bmap[b] = copyblock!(out, block(bl.ir, b)).id
  end
  ret = IRTools.block!(out)
  bmap[after] = ret.id
  for (arg, T) in zip(arguments(block(bl.ir, after)), argtypes(block(bl.ir, after)))
    argument!(ret, nothing, T, insert = false)
  end
  for b in blocks(out), (i, br) in enumerate(branches(b))
    branches(b)[i] = IRTools.Branch(br, block = bmap[br.block])
  end
  IRTools.domorder!(out)
  return out, args, after
end

function inline!(out, ir, args)
  env = Dict(zip(arguments(ir), args))
  before = blocks(out)[end]
  for br in branches(block(ir, 1))
    push!(branches(before), rename(env, br))
  end
  bmap = Dict()
  for bl in blocks(ir)[2:end]
    bmap[bl.id] = copyblock!(out, bl).id
  end
  for b in before.id:blocks(out)[end].id, (i, br) in enumerate(branches(block(out, b)))
    branches(block(out, b))[i] = IRTools.Branch(br, block = bmap[br.block])
  end
end

function abstract!(out, bl, env)
  ir, args, after = extract_group(bl)
  for i = 1:length(args)
    argtypes(ir)[i] = exprtype(out, rename(env, args[i]))
  end
  infer!(Inference(Frame(ir)))
  inline!(out, ir, rename.((env,), args))
  for (k, v) in zip(arguments(block(bl.ir, after)), arguments(blocks(out)[end]))
    env[k] = v
  end
  return after
end

function openbranches(out, env, bl)
  brs = []
  for br in branches(bl)
    br.condition == nothing && (push!(brs, br); break)
    cond = exprtype(out, env[br.condition])
    cond == Const(true) && continue
    cond == Const(false) && (push!(brs, br); break)
    push!(brs, br)
  end
  return brs
end

_nodetype(x, T::AType) = T
_nodetype(x, T::Type) = Node{T}(x)
nodetype(ir::IR, x) = _nodetype(x, exprtype(ir, x))

function traceblock!(out, env, bl)
  for (k, v) in bl
    ex = v.expr
    if isexpr(ex, :call)
      Ts = map(v -> nodetype(out, rename(env, v)), ex.args)
      if (T = partial(Ts...)) != nothing
        if T isa Node
          env[k] = T.value
        else
          env[k] = push!(out, stmt(rename(env, v.expr), type = T))
        end
      else
        env[k] = tracecall!(out, rename(env, ex).args, Ts)
      end
    elseif isexpr(ex, :meta)
    elseif isexpr(ex)
      error("Can't trace through $(ex.head) expression")
    else
      env[k] = ex
    end
  end
end

function trace!(out, ir, args)
  env = Dict{Any,Any}(zip(arguments(ir), args))
  bl = 1
  while true
    traceblock!(out, env, block(ir, bl))
    brs = openbranches(out, env, block(ir, bl))
    if length(brs) == 1
      isreturn(brs[1]) && return rename(env, returnvalue(brs[1]))
      bl = brs[1].block
      foreach((a, b) -> env[a] = rename(env, b), arguments(block(ir, bl)), arguments(brs[1]))
    else
      bl = abstract!(out, block(ir, bl), env)
    end
  end
end

function tracecall!(tr, args, Ts)
  ir = IR(widen.(Ts)...)
  ir == nothing && error("No IR for $(Tuple{widen.(Ts)...})")
  ir = ir |> merge_returns! |> prepare_ir!
  trace!(tr, ir, args)
end

function trace(Ts...)
  tr = IR()
  args = [argument!(tr, T) for T in Ts]
  return!(tr, tracecall!(tr, args, Ts))
  return cleanup!(tr)
end

atype(T::AType) = T
atype(x) = Const(x)

macro trace(ex)
  @capture(ex, f_(args__)) || error("@trace f(args...)")
  :(trace(atype.(($(esc(f)), $(esc.(args)...)))...))
end

ext: jl


text:
# We create a simple function to test
A = rand(5, 5);
Q = A' * A

f(x) = x' * Q * x
nlp = ADNLPModel(f, zeros(5))
nlp_at_x = NLPAtX(zeros(5))
stop_nlp = NLPStopping(nlp, (x,y) -> unconstrained_check(x,y), nlp_at_x, optimality0 = 0.0)

a = zeros(5)
fill_in!(stop_nlp, a)

# we make sure that the fill_in! function works properly
@test obj(nlp, a) == stop_nlp.current_state.fx
@test grad(nlp, a) == stop_nlp.current_state.gx
@test stop_nlp.meta.optimality0 == 0.0

# we make sure the optimality check works properly
@test stop!(stop_nlp)
@test stop_nlp.current_state.current_score == unconstrained_check(stop_nlp.pb, stop_nlp.current_state)
# we make sure the counter of stop works properly
@test stop_nlp.meta.nb_of_stop == 1

reinit!(stop_nlp, rstate = true, x = ones(5))
@test stop_nlp.current_state.x == ones(5)
@test stop_nlp.current_state.fx == nothing
@test stop_nlp.meta.nb_of_stop == 0

#We know test how to initialize the counter:
test_max_cntrs = Stopping._init_max_counters(obj = 2)
stop_nlp_cntrs = NLPStopping(nlp, max_cntrs = test_max_cntrs)
@test stop_nlp_cntrs.max_cntrs[:neval_obj] == 2
@test stop_nlp_cntrs.max_cntrs[:neval_grad] == 20000
@test stop_nlp_cntrs.max_cntrs[:neval_sum] == 20000*11

reinit!(stop_nlp.current_state)
@test unconstrained_check(stop_nlp.pb, stop_nlp.current_state) >= 0.0
reinit!(stop_nlp.current_state)
@test unconstrained2nd_check(stop_nlp.pb, stop_nlp.current_state) >= 0.0
@test stop_nlp.current_state.Hx != nothing

nlp_bnd = ADNLPModel(f, zeros(5), lvar = zeros(5), uvar = zeros(5))
stop_bnd = NLPStopping(nlp_bnd)
fill_in!(stop_bnd, zeros(5))
@test KKT(stop_bnd.pb, stop_bnd.current_state) == 0.0
reinit!(stop_bnd.current_state)
@test optim_check_bounded(stop_bnd.pb, stop_bnd.current_state) == 0.0

stop_bnd.optimality_check = (x,y) -> NaN
start!(stop_bnd)
@test stop_bnd.meta.domainerror == true
reinit!(stop_bnd)
@test stop_bnd.meta.domainerror == false
stop!(stop_bnd)
@test stop_bnd.meta.domainerror == true

stop_bnd.optimality_check = (x,y) -> 0.0
reinit!(stop_bnd)
fill_in!(stop_bnd, zeros(5), mu = ones(5), lambda = zeros(0))
@test stop_bnd.current_state.mu == ones(5)
@test stop_bnd.current_state.lambda == zeros(0)

update!(stop_bnd.current_state, fx = nothing, current_score = 0.0)
stop!(stop_bnd)
@test stop_bnd.current_state.fx != nothing

try
 NLPStopping(nlp_bnd, KKT, GenericState(x0)) #would not work as there is an entry check
 @test false
catch
 @test true
end


text:
"""
`rescale(vec0, [newlim0, newlim1])` rescale `vec0` on the new scale `[newlim0, newlim1]`.
"""
function rescale(vec0::Vector{<:Real}, vec2::Vector{<:Real})
	if length(vec2) != 2
		error("`vec2` should be a vector of exact two elements.");
	end
	intervals = vec0 |> diff;
	interval1 = intervals |> sum;
	diffv2 = diff(vec2);
	interval2 =  pop!(diffv2);

	sf = interval2/interval1; # scale factor
	int_scaled = sf*intervals;
	vec1 = [first(vec2), int_scaled...] |> cumsum
	return vec1
end


text:
NCdata = Union{NCDataset, NCDatasets.MFDataset}
NCfiles = Union{Vector{<:AbstractString}, AbstractString}

# nc_open = NCDataset
"""
    nc_open is same as NCDataset

@seealso [NCDatasets.NCDataset()]
"""
# nc_open = NCDataset

function nc_open(f::NCfiles, args...; kwargs...)
    NCDataset(path_mnt.(f), args...; kwargs...)
end

function nc_open(f::Function, args...; kwargs...)
    ds = nc_open(args...; kwargs...)
    try
        f(ds)
    finally
        @debug "closing netCDF NCDataset $(ds.ncid) $(NCDatasets.path(ds))"
        close(ds)
    end
end

nc_close(ds::NCdata) = close(ds)

# https://github.com/rafaqz/Rasters.jl/blob/master/src/sources/ncdatasets.jl
function nc_bands(ds::NCdata)
    # v_id = NCDatasets.nc_inq_varids(ds.ncid)
    # vars = NCDatasets.nc_inq_varname.(ds.ncid, v_id)
    vars = keys(ds)
    dims = ["lon", "long", "longitude",
        "lat", "latitude",
        "lev", "level", "mlev",
        "crs",
        # "vertical",
        # "x", "y", "z",
        "time"]
    setdiff(vars, [dims; dims .* "_bnds"; "height"])
end

function nc_bands(file::NCfiles)
    nc_open(file) do ds
        nc_bands(ds)
    end
end

function nc_info(ds::NCdata)
    # vars = nc_bands(ds)[1]
    println(ds)
end

function nc_info(file::NCfiles)
    println(basename(file))
    nc_open(file) do ds
        nc_info(ds)
    end
end

ncinfo = nc_info


text:
_lm(::Type{LabelEnc.TrueFalse},   ::Type{Val{2}}) = LabelEnc.TrueFalse()
_lm(::Type{LabelEnc.ZeroOne},     ::Type{Val{2}}) = LabelEnc.ZeroOne()
_lm(::Type{LabelEnc.MarginBased}, ::Type{Val{2}}) = LabelEnc.MarginBased()
_lm(::Type{LabelEnc.OneOfK},      ::Type{Val{K}}) where {K} = LabelEnc.OneOfK(Val{K})
_lm(::Type{LabelEnc.Indices},     ::Type{Val{K}}) where {K} = LabelEnc.Indices(Val{K})
_lm(::Type{LabelEnc.OneOfK{D}},   ::Type{Val{K}}) where {D,K} = LabelEnc.OneOfK(D,Val{K})
_lm(::Type{LabelEnc.Indices{D}},  ::Type{Val{K}}) where {D,K} = LabelEnc.Indices(D,Val{K})

_lm(::Type{LabelEnc.TrueFalse}, ::Type{T}, ::Type{Val{2}}) where {T} = LabelEnc.TrueFalse()
_lm(::Type{D}, ::Type, ::Type{Val{2}}) where {D<:LabelEnc.ZeroOne} = D()
_lm(::Type{D}, ::Type, ::Type{Val{2}}) where {D<:LabelEnc.MarginBased} = D()
_lm(::Type{LabelEnc.ZeroOne}, ::Type{T}, ::Type{Val{2}}) where {T<:Number} = LabelEnc.ZeroOne(T)
_lm(::Type{LabelEnc.MarginBased}, ::Type{T}, ::Type{Val{2}}) where {T<:Number} = LabelEnc.MarginBased(T)
_lm(::Type{LabelEnc.OneOfK},  ::Type{T}, ::Type{Val{K}}) where {T<:Number,K} = LabelEnc.OneOfK(T,Val{K})
_lm(::Type{LabelEnc.Indices}, ::Type{T}, ::Type{Val{K}}) where {T<:Number,K} = LabelEnc.Indices(T,Val{K})
_lm(::Type{LabelEnc.OneOfK},  ::Type, ::Type{Val{K}}) where {K} = LabelEnc.OneOfK(Val{K})
_lm(::Type{LabelEnc.Indices}, ::Type, ::Type{Val{K}}) where {K} = LabelEnc.Indices(Val{K})
_lm(::Type{LabelEnc.OneOfK{D}},  ::Type, ::Type{Val{K}}) where {D,K} = LabelEnc.OneOfK(D,Val{K})
_lm(::Type{LabelEnc.Indices{D}}, ::Type, ::Type{Val{K}}) where {D,K} = LabelEnc.Indices(D,Val{K})

@inline _array_type(::Type{T}, ::Type{Val{N}}) where {T,N} = Array{T,N}
@inline _array_type(::Type{Bool}, ::Type{Val{N}}) where {N} = BitArray{N}

## Convert views for Vector based using MappedArrays.

convertlabelview(dst::LabelEnc.FuzzyBinary, values, src::VectorLabelEncoding) = throw(MethodError(convertlabelview, (dst,values,src)))
convertlabelview(dst::VectorLabelEncoding,  values, src::LabelEnc.FuzzyBinary) = throw(MethodError(convertlabelview, (dst,values,src)))

function convertlabelview(dst::VectorLabelEncoding, values::AbstractVector)
    convertlabelview(dst, values, labelenc(values))
end

function convertlabelview(dst::VectorLabelEncoding{T,2}, values::AbstractVector{V}, src::LabelEnc.OneVsRest{V}) where {T,V}
    f = x -> convertlabel(dst, x, src)
    ReadonlyMappedArray{T,1,typeof(values),typeof(f)}(f, values)
end

function convertlabelview(dst::VectorLabelEncoding{T,K}, values::AbstractVector{V}, src::VectorLabelEncoding{V,K}) where {T,K,V}
    f    = x -> convertlabel(dst, x, src)
    finv = x -> convertlabel(src, x, dst)
    MappedArray{T,1,typeof(values),typeof(f),typeof(finv)}(f, finv, values)
end

## General Vector based

function convertlabel(dst::VectorLabelEncoding{T,K}, x, src::VectorLabelEncoding{S,K}) where {T,K,S}
    ind2label(label2ind(x, src), dst)::T
end

function convertlabel(dst::VectorLabelEncoding{T,K}, values::AbstractVector, src::VectorLabelEncoding{S,K}) where {T,K,S}
    convertlabel.(dst, values, src)::_array_type(T,Val{1})
end

## Generic types to objects

function convertlabel(::Type{L}, x, src::LabelEncoding{T,K}) where {L<:LabelEncoding,T,K}
    convertlabel(_lm(L,Val{K}), x, src)
end

function convertlabel(::Type{L}, values::AbstractArray{Bool}, src::LabelEncoding{T,K}, args...) where {L<:LabelEncoding,T,K}
    convertlabel(_lm(L,Val{K}), values, src, args...)
end

function convertlabel(::Type{L}, values::AbstractArray{T}, src::LabelEncoding{S,K}, args...) where {L<:LabelEncoding,T,S,K}
    convertlabel(_lm(L,T,Val{K}), values, src, args...)
end

function convertlabel(dst::LabelEncoding{T,K,M}, values::AbstractMatrix) where {T,K,M}
    convertlabel(dst, values, labelenc(values))::_array_type(T,Val{M})
end

function convertlabel(dst::LabelEncoding{T,K,M}, values::AbstractVector) where {T,K,M} # avoid method clash
    convertlabel(dst, values, labelenc(values))::_array_type(T,Val{M})
end

convertlabel(dst, values::AbstractArray) = convertlabel(dst, values, labelenc(values))

## OneVsRest

function convertlabel(dst::LabelEnc.OneVsRest{T}, values::AbstractVector{T}) where {T}
    convertlabel(dst, values, dst)
end

## NativeLabels

function convertlabel(dst, values, src_lbl::AbstractVector, args...)
    convertlabel(dst, values, LabelEnc.NativeLabels(src_lbl), args...)
end

# NativeLabels binary inference helper

function convertlabel(dst_lbl::AbstractVector{T}, values, src::LabelEncoding{S,K}, args...) where {T,S,K}
    convertlabel(LabelEnc.NativeLabels{T,K}(dst_lbl), values, src, args...)
end

function convertlabel(dst::LabelEncoding{S,K}, values, src_lbl::AbstractVector{T}, args...) where {S,K,T}
    convertlabel(dst, values, LabelEnc.NativeLabels{T,K}(src_lbl), args...)
end

function convertlabel(::Type{L}, values, src_lbl::AbstractVector{T}, args...) where {L<:BinaryLabelEncoding,T}
    convertlabel(L, values, LabelEnc.NativeLabels{T,2}(src_lbl), args...)
end

## OneOfK obsdim kw specified

convertlabel(dst, values, src; obsdim = LearnBase.default_obsdim(values)) = convertlabel(dst, values, src, convert(LearnBase.ObsDimension,obsdim))

convertlabel(dst, values; obsdim = LearnBase.default_obsdim(values)) = convertlabel(dst, values, labelenc(values), convert(LearnBase.ObsDimension,obsdim))

function convertlabel(dst, values::AbstractMatrix; obsdim = LearnBase.default_obsdim(values))
    nobsdim = convert(LearnBase.ObsDimension,obsdim)
    convertlabel(dst, values, labelenc(values, nobsdim), nobsdim)
end

# OneOfK default obsdim inserted

function convertlabel(dst::LabelEnc.OneOfK, values::AbstractMatrix, src::LabelEnc.OneOfK)
    convertlabel(dst, values, src, LearnBase.default_obsdim(values))
end

function convertlabel(dst::LabelEncoding, values::AbstractMatrix, src::LabelEnc.OneOfK)
    convertlabel(dst, values, src, LearnBase.default_obsdim(values))
end

function convertlabel(dst::LabelEnc.OneOfK, values::AbstractArray, src::LabelEncoding)
    convertlabel(dst, values, src, LearnBase.default_obsdim(values))
end

## To OneOfK

function convertlabel(dst::LabelEnc.OneOfK{TD,K}, values::AbstractMatrix, src::LabelEnc.OneOfK{TS,K}, ::LearnBase.ObsDimension) where {TD,TS,K}
    _array_type(TD,Val{2})(values)
end

function convertlabel(dst::LabelEnc.OneOfK{T,K}, values::AbstractVector, src::VectorLabelEncoding{S,KS}, ::Union{ObsDim.Last,ObsDim.Constant{2}}) where {T,K,S,KS}
    @assert KS <= K
    n = length(values)
    buffer = _array_type(T,Val{2})(zeros(T, K, n))
    @inbounds for i in 1:n
        buffer[label2ind(values[i], src), i] = one(T)
    end
    buffer
end

function convertlabel(dst::LabelEnc.OneOfK{T,K}, values::AbstractVector, src::VectorLabelEncoding{S,KS}, ::ObsDim.First) where {T,K,S,KS}
    @assert KS <= K
    n = length(values)
    buffer = _array_type(T,Val{2})(zeros(T, n, K))
    @inbounds for i in 1:n
        buffer[i, label2ind(values[i], src)] = one(T)
    end
    buffer
end

## From OneOfK

function convertlabel(dst::VectorLabelEncoding{TD,K}, values::AbstractMatrix{T}, src::LabelEnc.OneOfK{TS,K}, ::Union{ObsDim.Last,ObsDim.Constant{2}}) where {TD,T,TS,K}
    @assert size(values, 1) == K
    n = size(values, 2)
    buffer = _array_type(TD,Val{1})(undef, n)
    @inbounds for i in 1:n
        tind = 1
        tmax = typemin(T)
        for j in 1:K
            tval = values[j,i]
            if tval > tmax
                tind = j
                tmax = tval
            end
        end
        buffer[i] = ind2label(tind, dst)
    end
    buffer
end

function convertlabel(dst::VectorLabelEncoding{TD,K}, values::AbstractMatrix{T}, src::LabelEnc.OneOfK{TS,K}, ::ObsDim.First) where {TD,T,TS,K}
    @assert size(values, 2) == K
    n = size(values, 1)
    buffer = _array_type(TD,Val{1})(undef, n)
    @inbounds for i in 1:n
        tind = 1
        tmax = typemin(T)
        for j in 1:K
            tval = values[i,j]
            if tval > tmax
                tind = j
                tmax = tval
            end
        end
        buffer[i] = ind2label(tind, dst)
    end
    buffer
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 0120e370-f8f5-11ea-1b65-9d29b0b5b647
begin
	using LightGraphs
	using LightGraphs.SimpleGraphs
	using MetaGraphs
	using CitationNetworks
	using GraphPlot
end

# ╔═╡ 0ba77020-f8f5-11ea-025d-61768c938c5c
A = [
 0  0  1  0  0  0  0  0  0  0  0
 0  0  1  1  0  0  0  0  0  0  0
 0  0  0  0  1  0  1  0  0  0  0
 0  0  0  0  0  1  0  0  0  1  1
 0  0  0  0  0  0  0  1  1  0  0
 0  0  0  0  0  0  0  0  1  1  0
 0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0
];

# ╔═╡ 2680e980-f8f5-11ea-10c2-fbf1d8358506
g = SimpleDiGraph(A)

# ╔═╡ 53e56540-f8f5-11ea-1189-99d9421b9146
gplot(g, nodelabel=1:nv(g))

# ╔═╡ 7563a7d0-f8f6-11ea-00ae-bbee824a385c
genetic_knowper(g, normalize=:none)

# ╔═╡ 6b60e68e-f8f5-11ea-25f9-f3075e37b317
mpres = mainpath(g, 3, GBFP(), normalize=:none)

# ╔═╡ 3404b0a0-f8f5-11ea-1cb0-4f35dff35a1f
mpnet, vlabels = induced_subgraph(g, mpres.edges)

# ╔═╡ bd8334f0-f8f5-11ea-0def-dda12adf86e7
gplot(mpnet, nodelabel=vlabels)

# ╔═╡ 5bb78380-f8f9-11ea-2b66-1529d729d0d9
function add_weights(g, w)
    mg = MetaDiGraph(g)
	for (i,e) in enumerate(edges(g))
		set_prop!(mg, e, :weight, w[i])
	end
	return mg
end

# ╔═╡ c9fba690-f8f5-11ea-068f-7393d3245730
begin
	spc = weights_spc(g).edgeweights
	mg = add_weights(g, spc)
	mpfl = mainpath(mg, ForwardLocal())
	mpbl = mainpath(mg, BackwardLocal())
end

# ╔═╡ Cell order:
# ╠═0120e370-f8f5-11ea-1b65-9d29b0b5b647
# ╠═0ba77020-f8f5-11ea-025d-61768c938c5c
# ╠═2680e980-f8f5-11ea-10c2-fbf1d8358506
# ╠═53e56540-f8f5-11ea-1189-99d9421b9146
# ╠═7563a7d0-f8f6-11ea-00ae-bbee824a385c
# ╠═6b60e68e-f8f5-11ea-25f9-f3075e37b317
# ╠═3404b0a0-f8f5-11ea-1cb0-4f35dff35a1f
# ╠═bd8334f0-f8f5-11ea-0def-dda12adf86e7
# ╠═5bb78380-f8f9-11ea-2b66-1529d729d0d9
# ╠═c9fba690-f8f5-11ea-068f-7393d3245730


text:
# This is a ripoff of the sinpi/cospi functions from Base.

function sintau(x::Real)
    if isinf(x)
        return throw(DomainError())
    elseif isnan(x)
        return nan(x)
    end

    rx = copysign(float(rem(x,1)),x)
    arx = abs(rx)

    if arx <= oftype(rx,1/8)
        return sin(tau*rx)
    elseif arx <= oftype(rx,3/8)
        arx = oftype(rx,1/4) - arx
        return copysign(cos(tau*arx),rx)
    elseif arx < oftype(rx,5/8)
        rx = (oftype(rx,1/2) - arx)*sign(rx)
        return sin(tau*rx)
    elseif arx <= oftype(rx,7/8)
        arx = oftype(rx,3/4) - arx
        return -copysign(cos(tau*arx),rx)
    else
        rx = rx - copysign(one(rx),rx)
        return sin(tau*rx)
    end
end

function costau(x::Real)
    if isinf(x)
        return throw(DomainError())
    elseif isnan(x)
        return nan(x)
    end

    rx = abs(float(rem(x,1)))

    if rx <= oftype(rx,1/8)
        return cos(tau*rx)
    elseif rx < oftype(rx,3/8)
        rx = oftype(rx,1/4) - rx
        return sin(tau*rx)
    elseif rx <= oftype(rx,5/8)
        rx = oftype(rx,1/2) - rx
        return -cos(tau*rx)
    elseif rx < oftype(rx,7/8)
        rx = rx - oftype(rx,3/4)
        return sin(tau*rx)
    else
        rx = one(rx) - rx
        return cos(tau*rx)
    end
end

sintau(x::Integer) = zero(x)
costau(x::Integer) = one(x)

function sintau(z::Complex)
    zr, zi = reim(z)
    if !isfinite(zi) && zr == 0 return complex(zr, zi) end
    if isnan(zr) && !isfinite(zi) return complex(zr, zi) end
    if !isfinite(zr) && zi == 0 return complex(oftype(zr, NaN), zi) end
    if !isfinite(zr) && isfinite(zi) return complex(oftype(zr, NaN), oftype(zi, NaN)) end
    if !isfinite(zr) && !isfinite(zi) return complex(zr, oftype(zi, NaN)) end
    tauzi = tau*zi
    complex(sintau(zr)*cosh(tauzi), costau(zr)*sinh(tauzi))
end

function costau(z::Complex)
    zr, zi = reim(z)
    if !isfinite(zi) && zr == 0
        return complex(isnan(zi) ? zi : oftype(zi, Inf),
                       isnan(zi) ? zr : zr*-sign(zi))
    end
    if !isfinite(zr) && isinf(zi)
        return complex(oftype(zr, Inf), oftype(zi, NaN))
    end
    if isinf(zr)
        return complex(oftype(zr, NaN), zi==0 ? -copysign(zi, zr) : oftype(zi, NaN))
    end
    if isnan(zr) && zi==0 return complex(zr, abs(zi)) end
    tauzi = tau*zi
    complex(costau(zr)*cosh(tauzi), -sintau(zr)*sinh(tauzi))
end

Base.@vectorize_1arg Number sintau
Base.@vectorize_1arg Number costau

const sinτ = sintau
const cosτ = costau


text:
module Lattice1D
    using SparseArrays
    export H0, Hnoise, Hamiltonian
    """
        H0(N; ϵ=0., ν=1.)

    Constructs the noiseless Hamiltonian for the 1D lattice with `N` nodes,
    self-energy ``\\epsilon`` and tunneling constant ``\\nu``.
    """
    function H0(n; ϵ=0., ν=1.)
        H = spdiagm(0 => ϵ*ones(n), 1=>ν*ones(n-1), -1=> ν*ones(n-1))
        if n > 2
            H[1,end] = H[end,1] = ν
        end
        return H
    end

    """
        Hnoise(N, noise_id)

    Returns the noisy part of the Hamiltonian for the 1D lattice with `N` nodes
    for a specific state of the noise fluctuators specified by `noise_id`.
    """
    function Hnoise(n, noise_id)
        if noise_id < 1 || noise_id > 2^n
            throw(ArgumentError("noise_id must be between 1 and 2^n"))
        end
        noise = 2 * digits(noise_id .- 1, base=2, pad=n) .- 1
        H = spdiagm(-1 => noise[2:end], 1 => noise[2:end])
        H[1, end] = H[end, 1] = noise[1]
        return H
    end

    """
        Hamiltonian(n, ν, ν0=1.)

    Returns an array of Hamiltonians of size `n` for all possible configurations
    of noise.

    # Arguments
        * `n::Integer`: the dimension of the Hilbert space
        * `ν::Float`: the strength of the noise
        * `ν0::Float`: The rate of the noiseless Hamiltonian (ν0 = 1.)
    """
    function Hamiltonian(n::Integer, ν::Float64, ν0=1.)
        return [H0(n;ν=ν0) + ν*Hnoise(n, id) for id = 1:2^n]
    end

end


text:
const PyJuliaVectorValue_Type__ref = Ref(PyNULL)
PyJuliaVectorValue_Type() = begin
    ptr = PyJuliaVectorValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaArrayValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.VectorValue",
                base = base,
                methods = [
                    (name = "resize", flags = Py_METH_O, meth = pyjlvector_resize),
                    (
                        name = "sort",
                        flags = Py_METH_KEYWORDS | Py_METH_VARARGS,
                        meth = pyjlvector_sort,
                    ),
                    (name = "reverse", flags = Py_METH_NOARGS, meth = pyjlvector_reverse),
                    (name = "clear", flags = Py_METH_NOARGS, meth = pyjlvector_clear),
                    (
                        name = "__reversed__",
                        flags = Py_METH_NOARGS,
                        meth = pyjlvector_reversed,
                    ),
                    (name = "insert", flags = Py_METH_VARARGS, meth = pyjlvector_insert),
                    (name = "append", flags = Py_METH_O, meth = pyjlvector_append),
                    (name = "extend", flags = Py_METH_O, meth = pyjlvector_extend),
                    (name = "pop", flags = Py_METH_VARARGS, meth = pyjlvector_pop),
                    (name = "remove", flags = Py_METH_O, meth = pyjlvector_remove),
                    (name = "index", flags = Py_METH_O, meth = pyjlvector_index),
                    (name = "count", flags = Py_METH_O, meth = pyjlvector_count),
                ],
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        abc = PyMutableSequenceABC_Type()
        isnull(abc) && return PyNULL
        ism1(PyABC_Register(ptr, abc)) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaVectorValue_Type__ref[] = ptr
    end
    ptr
end

PyJuliaVectorValue_New(x::AbstractVector) = PyJuliaValue_New(PyJuliaVectorValue_Type(), x)
PyJuliaValue_From(x::AbstractVector) = PyJuliaVectorValue_New(x)

pyjlvector_resize(xo::PyPtr, arg::PyPtr) =
    try
        x = PyJuliaValue_GetValue(xo)::AbstractVector
        r = PyObject_TryConvert(arg, Int)
        r == -1 && return PyNULL
        r == 0 &&
            (PyErr_SetString(PyExc_TypeError(), "size must be an integer"); return PyNULL)
        resize!(x, takeresult(Int))
        PyNone_New()
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlvector_sort(xo::PyPtr, args::PyPtr, kwargs::PyPtr) =
    try
        x = PyJuliaValue_GetValue(xo)::AbstractVector
        ism1(PyArg_CheckNumArgsEq("sort", args, 0)) && return PyNULL
        ism1(PyArg_GetArg(Bool, "sort", kwargs, "reverse", false)) && return PyNULL
        rev = takeresult(Bool)
        ism1(PyArg_GetArg(Any, "sort", kwargs, "key", nothing)) && return PyNULL
        by = takeresult()
        sort!(x, rev = rev, by = by === nothing ? identity : by)
        PyNone_New()
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlvector_reverse(xo::PyPtr, ::PyPtr) =
    try
        x = PyJuliaValue_GetValue(xo)::AbstractVector
        reverse!(x)
        PyNone_New()
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlvector_clear(xo::PyPtr, ::PyPtr) =
    try
        x = PyJuliaValue_GetValue(xo)::AbstractVector
        empty!(x)
        PyNone_New()
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlvector_reversed(xo::PyPtr, ::PyPtr) =
    try
        x = PyJuliaValue_GetValue(xo)::AbstractVector
        PyObject_From(reverse(x))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlvector_insert(xo::PyPtr, args::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractVector
    a = axes(x, 1)
    ism1(PyArg_CheckNumArgsEq("insert", args, 2)) && return PyNULL
    ism1(PyArg_GetArg(Int, "insert", args, 0)) && return PyNULL
    k = takeresult(Int)
    k′ = k < 0 ? (last(a) + 1 + k) : (first(a) + k)
    checkbounds(Bool, x, k′) || (
        PyErr_SetString(PyExc_IndexError(), "array index out of bounds"); return PyNULL
    )
    ism1(PyArg_GetArg(eltype(x), "insert", args, 1)) && return PyNULL
    v = takeresult(eltype(x))
    try
        insert!(x, k′, v)
        PyNone_New()
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjlvector_append(xo::PyPtr, vo::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractVector
    r = PyObject_TryConvert(vo, eltype(x))
    r == -1 && return PyNULL
    r == 0 && (
        PyErr_SetString(
            PyExc_TypeError(),
            "array value must be a Julia '$(eltype(x))', not a '$(PyType_Name(Py_Type(vo)))'",
        );
        return PyNULL
    )
    v = takeresult(eltype(x))
    try
        push!(x, v)
        PyNone_New()
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjlvector_extend(xo::PyPtr, vso::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractVector
    r = PyIterable_Map(vso) do vo
        r = PyObject_Convert(vo, eltype(x))
        r == -1 && return -1
        v = takeresult(eltype(x))
        push!(x, v)
        return 1
    end
    r == -1 ? PyNULL : PyNone_New()
end

pyjlvector_pop(xo::PyPtr, args::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractVector
    a = axes(x, 1)
    ism1(PyArg_CheckNumArgsLe("pop", args, 1)) && return PyNULL
    ism1(PyArg_GetArg(Int, "pop", args, 0, -1)) && return PyNULL
    k = takeresult(Int)
    k′ = k < 0 ? (last(a) + 1 + k) : (first(a) + k)
    checkbounds(Bool, x, k′) || (
        PyErr_SetString(PyExc_IndexError(), "array index out of bounds"); return PyNULL
    )
    try
        if k′ == last(a)
            v = pop!(x)
        elseif k′ == first(a)
            v = popfirst!(x)
        else
            v = x[k′]
            deleteat!(x, k′)
        end
        PyObject_From(v)
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjlvector_remove(xo::PyPtr, vo::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractVector
    r = PyObject_TryConvert(vo, eltype(x))
    r == -1 && return PyNULL
    r == 0 &&
        (PyErr_SetString(PyExc_ValueError(), "value not in vector"); return PyNULL)
    v = takeresult(eltype(x))
    try
        k = findfirst(x -> x == v, x)
        if k === nothing
            PyErr_SetString(PyExc_ValueError(), "value not in vector")
            PyNULL
        else
            deleteat!(x, k)
            PyNone_New()
        end
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjlvector_index(xo::PyPtr, vo::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractVector
    r = PyObject_TryConvert(vo, eltype(x))
    r == -1 && return PyNULL
    r == 0 &&
        (PyErr_SetString(PyExc_ValueError(), "value not in vector"); return PyNULL)
    v = takeresult(eltype(x))
    try
        k = findfirst(x -> x == v, x)
        if k === nothing
            PyErr_SetString(PyExc_ValueError(), "value not in vector")
            PyNULL
        else
            PyObject_From(k - first(axes(x, 1)))
        end
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjlvector_count(xo::PyPtr, vo::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractVector
    r = PyObject_TryConvert(vo, eltype(x))
    r == -1 && return PyNULL
    r == 0 && return PyObject_From(0)
    v = takeresult(eltype(x))
    try
        n = count(x -> x == v, x)
        PyObject_From(n)
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end


text:
"""Extension of QuantumPropagators.propagate for control objectives."""

import QuantumPropagators
import .ConditionalThreads: @threadsif

"""Construct a `genfunc` suitable for propagating an objective.

```julia
genfunc = objective_genfunc(obj, tlist; controls_map=IdDict())
```

can be passed to [`QuantumPropagators.propagate`](@ref) to propagate under the
dynamical generator in `obj`.

If given, `control_map` defines replacements for the controls in
`obj.generator`. This allows, e.g., to replace the controls with those
resulting from an optimization.
"""
function objective_genfunc(obj::AbstractControlObjective, tlist; controls_map=IdDict())
    controls = getcontrols(obj.generator)
    pulses = [
        discretize_on_midpoints(get(controls_map, control, control), tlist) for
        control in controls
    ]

    zero_vals = IdDict(control => 0 for control in controls)
    G = evalcontrols(obj.generator, zero_vals)
    vals_dict = IdDict(control => pulses[j][1] for (j, control) in enumerate(controls))

    function genfunc(tlist, i; kwargs...)
        for (j, control) in enumerate(controls)
            vals_dict[control] = pulses[j][i]
        end
        evalcontrols!(G, obj.generator, vals_dict)
        return G
    end

    return genfunc
end

# Internal method for getting a propagation method from one or more keyword
# arguments, with a default in a corresponding field/property of `obj`. For use
# in setting up the workspace for optimization methods.
#
# For example,
#
#     get_objective_prop_method(obj, :bw_prop_method, :prop_method; kwargs...)
#
# uses `kwargs[:bw_prop_method]`, `obj.bw_prop_method`, `kwargs[:prop_method]`,
# `obj.prop_method`, `:auto`, in that order.
#
function get_objective_prop_method(obj, symbols...; kwargs...)
    for symbol in symbols
        if symbol ∈ keys(kwargs)
            return kwargs[symbol]
        elseif symbol ∈ propertynames(obj)
            return Symbol(getproperty(obj, symbol))
        end
    end
    return :auto
end

"""Propagate with the dynamical generator of a control objective.

```julia
propagate_objective(obj, tlist; method=:auto, initial_state=obj.initial_state,
                    controls_map=IdDict(), kwargs...)
```

propagates `initial_state` under the dynamics described by `obj.generator`.

The optional dict `control_map` may be given to replace the controls in
`obj.generator` (as obtained by [`getcontrols`](@ref)) with custom functions or
vectors, e.g. with the controls resulting from optimization.

If `obj` has a property/field `prop_method` or `fw_prop_method`, its value will
be used as the default for `method` instead of :auto. An explicit keyword
argument for `method` always overrides the default.

All other `kwargs` are forwarded to the underlying
[`QuantumPropagators.propagate`](@ref) method for `obj.initial_state`.
"""
function propagate_objective(
    obj,
    tlist;
    initial_state=obj.initial_state,
    controls_map=IdDict(),
    kwargs...
)
    method = :auto
    if :method in keys(kwargs)
        method = kwargs[:method]
    else
        for symbol ∈ (:prop_method, :fw_prop_method)
            if symbol ∈ propertynames(obj)
                method = Symbol(getproperty(obj, symbol))
            end
        end
    end
    return propagate_objective(
        obj,
        tlist,
        Val(method);
        initial_state=initial_state,
        controls_map=controls_map,
        kwargs...
    )
end

#! format: off
function propagate_objective(obj, tlist, method::Symbol; initial_state,
                             controls_map=IdDict(), kwargs...)
    return propagate_objective(obj, tlist, Val(method);
                               initial_state=initial_state,
                               controls_map=controls_map, kwargs...)
end

function propagate_objective(obj, tlist, method::Val; initial_state,
                             controls_map=IdDict(), kwargs...)
    wrk = initobjpropwrk(obj, tlist, method; initial_state=initial_state,
                         kwargs...)
    return _propagate_objective(
        obj, tlist, wrk; initial_state=initial_state,
        controls_map=controls_map, kwargs...
    )

end
#! format: on

# `propagate_objective` backend (note `wrk` argument instead of `method`)
function _propagate_objective(
    obj,
    tlist,
    wrk;
    initial_state=obj.initial_state,
    controls_map=IdDict(),
    kwargs...
)
    genfunc = objective_genfunc(obj, tlist; controls_map=controls_map)
    return QuantumPropagators._propagate(initial_state, genfunc, tlist, wrk; kwargs...)
end

"""Propagate multiple objectives in parallel.

```julia
result = propagate_objectives(objectives, tlist; use_threads=true, kwargs...)
```

runs [`propagate_objective`](@ref) for every objective in `objectives`,
collects and returns a vector of results. The propagation happens in parallel
if `use_threads=true` (default). All keyword parameters are passed to
[`propagate_objective`](@ref), except that if `initial_state` is given, it must
be a vector of initial states, one for each objective. Likewise, to pass
pre-allocated storage arrays to `storage`, a vector of storage arrays must be
passed. A simple `storage=true` will still work to return a vector of storage
results.
"""
function propagate_objectives(
    objectives,
    tlist;
    use_threads=true,
    storage=nothing,
    initial_state=[obj.initial_state for obj in objectives],
    kwargs...
)
    result = Vector{Any}(undef, length(objectives))
    @threadsif use_threads for (k, obj) in collect(enumerate(objectives))
        if isnothing(storage) || (storage isa Bool)
            result[k] = propagate_objective(
                obj,
                tlist;
                storage=storage,
                initial_state=initial_state[k],
                kwargs...
            )
        else
            result[k] = propagate_objective(
                obj,
                tlist;
                storage=storage[k],
                initial_state=initial_state[k],
                kwargs...
            )
        end
    end
    return [result...]  # choose an atomatic eltype
end

"""
```julia
wrk = initobjpropwrk(obj, tlist, method; kwargs...)
```

initializes a workspace for the propagation of
an [`AbstractControlObjective`](@ref).
"""
function initobjpropwrk(
    obj::AbstractControlObjective,
    tlist,
    method::Val;
    initial_state=obj.initial_state,
    info_msg=nothing,
    verbose=false,
    kwargs...
)
    if verbose && !isnothing(info_msg)
        @info info_msg
    end
    controls = getcontrols(obj.generator)
    controlvals = [discretize(control, tlist) for control in controls]
    # We'll construct some examplary dynamical generators. This is mainly for
    # method=:cheby, but it's the "safe" thing to do for any method that uses
    # generators for initpropwrk, so we're using this as the default
    # implementation
    zero_vals =
        IdDict(control => zero(controlvals[i][1]) for (i, control) in enumerate(controls))
    G_zero = evalcontrols(obj.generator, zero_vals)
    max_vals =
        IdDict(control => maximum(controlvals[i]) for (i, control) in enumerate(controls))
    G_max = evalcontrols(obj.generator, max_vals)
    min_vals =
        IdDict(control => minimum(controlvals[i]) for (i, control) in enumerate(controls))
    G_min = evalcontrols(obj.generator, min_vals)
    return QuantumPropagators.initpropwrk(
        initial_state,
        tlist,
        method,
        G_zero,
        G_max,
        G_min;
        verbose=verbose,
        kwargs...
    )
end

function initobjpropwrk(
    obj::AbstractControlObjective,
    tlist,
    method::Symbol;
    initial_state=obj.initial_state,
    kwargs...
)
    return initobjpropwrk(obj, tlist, Val(method); initial_state=initial_state, kwargs...)
end

function initobjpropwrk(
    obj::AbstractControlObjective,
    tlist,
    method::Val{:newton};
    initial_state=obj.initial_state,
    kwargs...
)
    return QuantumPropagators.initpropwrk(initial_state, tlist, method; kwargs...)
end

function initobjpropwrk(
    obj::AbstractControlObjective,
    tlist,
    method::Val{:expprop};
    initial_state=obj.initial_state,
    kwargs...
)
    return QuantumPropagators.initpropwrk(initial_state, tlist, method; kwargs...)
end


text:
const ARDRegressor_ = sklm(:ARDRegression)
@sk_reg mutable struct ARDRegressor <: MMI.Deterministic
    n_iter::Int               = 300::(_ > 0)
    tol::Float64              = 1e-3::(_ > 0)
    alpha_1::Float64          = 1e-6::(_ > 0)
    alpha_2::Float64          = 1e-6::(_ > 0)
    lambda_1::Float64         = 1e-6::(_ > 0)
    lambda_2::Float64         = 1e-6::(_ > 0)
    compute_score::Bool       = false
    threshold_lambda::Float64 = 1e4::(_ > 0)
    fit_intercept::Bool       = true
    normalize::Bool           = false
    copy_X::Bool              = true
    verbose::Bool             = false
end
MMI.fitted_params(model::ARDRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    alpha     = fitresult.alpha_,
    lambda    = fitresult.lambda_,
    sigma     = fitresult.sigma_,
    scores    = fitresult.scores_
    )
descr(ARDRegressor, "Bayesian ARD regression.")

# =============================================================================
const BayesianRidgeRegressor_ = sklm(:BayesianRidge)
@sk_reg mutable struct BayesianRidgeRegressor <: MMI.Deterministic
    n_iter::Int         = 300::(_ ≥ 1)
    tol::Float64        = 1e-3::(_ > 0)
    alpha_1::Float64    = 1e-6::(_ > 0)
    alpha_2::Float64    = 1e-6::(_ > 0)
    lambda_1::Float64   = 1e-6::(_ > 0)
    lambda_2::Float64   = 1e-6::(_ > 0)
    compute_score::Bool = false
    fit_intercept::Bool = true
    normalize::Bool     = false
    copy_X::Bool        = true
    verbose::Bool       = false
end
MMI.fitted_params(model::BayesianRidgeRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    alpha     = fitresult.alpha_,
    lambda    = fitresult.lambda_,
    sigma     = fitresult.sigma_,
    scores    = fitresult.scores_
    )
descr(BayesianRidgeRegressor, "Bayesian ridge regression.")

# =============================================================================
const ElasticNetRegressor_ = sklm(:ElasticNet)
@sk_reg mutable struct ElasticNetRegressor <: MMI.Deterministic
    alpha::Float64      = 1.0::(_ ≥ 0)   # 0 is OLS
    l1_ratio::Float64   = 0.5::(0 ≤ _ ≤ 1)
    fit_intercept::Bool = true
    normalize::Bool     = false
    precompute::Union{Bool,AbstractMatrix} = false
    max_iter::Int       = 1_000::(_ ≥ 1)
    copy_X::Bool        = true
    tol::Float64        = 1e-4::(_ > 0)
    warm_start::Bool    = false
    positive::Bool      = false
    random_state::Any   = nothing  # Int, random state, or nothing
    selection::String   = "cyclic"::(_ in ("cyclic","random"))
end
MMI.fitted_params(model::ElasticNetRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    )
descr(ElasticNetRegressor, "Elastic net regression.")

# =============================================================================
const ElasticNetCVRegressor_ = sklm(:ElasticNetCV)
@sk_reg mutable struct ElasticNetCVRegressor <: MMI.Deterministic
    l1_ratio::Union{Float64,Vector{Float64}} = 0.5::(all(0 .≤ _ .≤ 1))
    eps::Float64        = 1e-3::(_ > 0)
    n_alphas::Int       = 100::(_ > 0)
    alphas::Any         = nothing::(_ === nothing || all(0 .≤ _))
    fit_intercept::Bool = true
    normalize::Bool     = false
    precompute::Union{Bool,String,AbstractMatrix} = "auto"
    max_iter::Int       = 1_000::(_ > 0)
    tol::Float64        = 1e-4::(_ > 0)
    cv::Any             = 5 # can be Int, Nothing or an iterable / cv splitter
    copy_X::Bool        = true
    verbose::Union{Bool, Int}  = 0
    n_jobs::Option{Int} = nothing
    positive::Bool      = false
    random_state::Any   = nothing
    selection::String   = "cyclic"::(_ in ("cyclic","random"))
end
MMI.fitted_params(model::ElasticNetCVRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    l1_ratio  = fitresult.l1_ratio_,
    mse_path  = fitresult.mse_path_,
    alphas    = fitresult.alphas_
    )
descr(ElasticNetCVRegressor, "Elastic net regression $CV.")

# =============================================================================
const HuberRegressor_ = sklm(:HuberRegressor)
@sk_reg mutable struct HuberRegressor <: MMI.Deterministic
    epsilon::Float64    = 1.35::(_ > 1.0)
    max_iter::Int       = 100::(_ > 0)
    alpha::Float64      = 1e-4::(_ > 0)
    warm_start::Bool    = false
    fit_intercept::Bool = true
    tol::Float64        = 1e-5::(_ > 0)
end
MMI.fitted_params(model::HuberRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    scale     = fitresult.scale_,
    outliers  = fitresult.outliers_
    )
descr(HuberRegressor, "Huber regression.")

# =============================================================================
const LarsRegressor_ = sklm(:Lars)
@sk_reg mutable struct LarsRegressor <: MMI.Deterministic
    fit_intercept::Bool      = true
    verbose::Union{Bool,Int} = false
    normalize::Bool          = true
    precompute::Union{Bool,String,AbstractMatrix} = "auto"
    n_nonzero_coefs::Int     = 500::(_ > 0)
    eps::Float64   = eps(Float64)::(_ > 0)
    copy_X::Bool   = true
    fit_path::Bool = true
end
MMI.fitted_params(model::LarsRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    alphas    = fitresult.alphas_,
    active    = fitresult.active_,
    coef_path = fitresult.coef_path_
    )
descr(LarsRegressor, "Lars regression.")

# =============================================================================
const LarsCVRegressor_ = sklm(:LarsCV)
@sk_reg mutable struct LarsCVRegressor <: MMI.Deterministic
    fit_intercept::Bool      = true
    verbose::Union{Bool,Int} = false
    max_iter::Int     = 500::(_ > 0)
    normalize::Bool   = true
    precompute::Union{Bool,String,AbstractMatrix} = "auto"
    cv::Any           = 5
    max_n_alphas::Int = 1_000::(_ > 0)
    n_jobs::Option{Int} = nothing
    eps::Float64      = eps(Float64)::(_ > 0)
    copy_X::Bool      = true
end
MMI.fitted_params(model::LarsCVRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    alpha     = fitresult.alpha_,
    alphas    = fitresult.alphas_,
    cv_alphas = fitresult.cv_alphas_,
    mse_path  = fitresult.mse_path_,
    coef_path = fitresult.coef_path_
    )
descr(LarsCVRegressor, "Lars regression $CV.")

# =============================================================================
const LassoRegressor_ = sklm(:Lasso)
@sk_reg mutable struct LassoRegressor <: MMI.Deterministic
    alpha::Float64      = 1.0::(_ ≥ 0) # should use alpha > 0 (or OLS)
    fit_intercept::Bool = true
    normalize::Bool     = false
    precompute::Union{Bool,AbstractMatrix} = false
    copy_X::Bool        = true
    max_iter::Int       = 1_000::(_ > 0)
    tol::Float64        = 1e-4::(_ > 0)
    warm_start::Bool    = false
    positive::Bool      = false
    random_state::Any   = nothing
    selection::String   = "cyclic"::(_ in ("cyclic","random"))
end
MMI.fitted_params(model::LassoRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    )
descr(LassoRegressor, "Lasso regression.")

# =============================================================================
const LassoCVRegressor_ = sklm(:LassoCV)
@sk_reg mutable struct LassoCVRegressor <: MMI.Deterministic
    eps::Float64        = 1e-3::(_ > 0)
    n_alphas::Int       = 100::(_ > 0)
    alphas::Any         = nothing::(_ === nothing || all(0 .≤ _))
    fit_intercept::Bool = true
    normalize::Bool     = false
    precompute::Union{Bool,String,AbstractMatrix} = "auto"
    max_iter::Int       = 1_000::(_ > 0)
    tol::Float64        = 1e-4::(_ > 0)
    copy_X::Bool        = true
    cv::Any             = 5
    verbose::Union{Bool,Int} = false
    n_jobs::Option{Int} = nothing
    positive::Bool      = false
    random_state::Any   = nothing
    selection::String   = "cyclic"::(_ in ("cyclic","random"))
end
MMI.fitted_params(model::LassoCVRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    alpha     = fitresult.alpha_,
    alphas    = fitresult.alphas_,
    mse_path  = fitresult.mse_path_,
    dual_gap  = fitresult.dual_gap_
    )
descr(LassoCVRegressor, "Lasso regression $CV.")

# =============================================================================
const LassoLarsRegressor_ = sklm(:LassoLars)
@sk_reg mutable struct LassoLarsRegressor <: MMI.Deterministic
    alpha::Float64      = 1.0::(_ ≥ 0) # 0 should be OLS
    fit_intercept::Bool = true
    verbose::Union{Bool, Int} = false
    normalize::Bool     = true
    precompute::Union{Bool,String,AbstractMatrix} = "auto"
    max_iter::Int       = 500::(_ > 0)
    eps::Float64        = eps(Float64)::(_ > 0)
    copy_X::Bool        = true
    fit_path::Bool      = true
    positive::Any       = false
end
MMI.fitted_params(model::LassoLarsRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    alphas    = fitresult.alphas_,
    active    = fitresult.active_,
    coef_path = fitresult.coef_path_
    )
descr(LassoLarsRegressor, "Lasso model fit with least angle regression (LARS).")

# =============================================================================
const LassoLarsCVRegressor_ = sklm(:LassoLarsCV)
@sk_reg mutable struct LassoLarsCVRegressor <: MMI.Deterministic
    fit_intercept::Bool = true
    verbose::Union{Bool, Int} = false
    max_iter::Int       = 500::(_ > 0)
    normalize::Bool     = true
    precompute::Union{Bool,String,AbstractMatrix} = "auto"
    cv::Any             = 5
    max_n_alphas::Int   = 1_000::(_ > 0)
    n_jobs::Option{Int} = nothing
    eps::Float64        = eps(Float64)::(_ > 0.0)
    copy_X::Bool        = true
    positive::Any       = false
end
MMI.fitted_params(model::LassoLarsCVRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    coef_path = fitresult.coef_path_,
    alpha     = fitresult.alpha_,
    alphas    = fitresult.alphas_,
    cv_alphas = fitresult.cv_alphas_,
    mse_path  = fitresult.mse_path_
    )
descr(LassoLarsCVRegressor, "Lasso model fit with least angle regression (LARS) $CV.")

# =============================================================================
const LassoLarsICRegressor_ = sklm(:LassoLarsIC)
@sk_reg mutable struct LassoLarsICRegressor <: MMI.Deterministic
    criterion::String   = "aic"::(_ in ("aic","bic"))
    fit_intercept::Bool = true
    verbose::Union{Bool, Int} = false
    normalize::Bool     = true
    precompute::Union{Bool,String,AbstractMatrix} = "auto"
    max_iter::Int       = 500::(_ > 0)
    eps::Float64        = eps(Float64)::(_ > 0.0)
    copy_X::Bool        = true
    positive::Any       = false
end
MMI.fitted_params(model::LassoLarsICRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    alpha     = fitresult.alpha_
    )
descr(LassoLarsICRegressor, "Lasso model with Lars using BIC or AIC for model selection.")

# =============================================================================
const LinearRegressor_ = sklm(:LinearRegression)
@sk_reg mutable struct LinearRegressor <: MMI.Deterministic
    fit_intercept::Bool = true
    normalize::Bool     = false
    copy_X::Bool        = true
    n_jobs::Option{Int} = nothing
end
MMI.fitted_params(model::LinearRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing)
    )
descr(LinearRegressor, "Ordinary least-square regression (OLS).")

# =============================================================================
const OrthogonalMatchingPursuitRegressor_ = sklm(:OrthogonalMatchingPursuit)
@sk_reg mutable struct OrthogonalMatchingPursuitRegressor <: MMI.Deterministic
    n_nonzero_coefs::Option{Int} = nothing
    tol::Option{Float64} = nothing
    fit_intercept::Bool  = true
    normalize::Bool      = true
    precompute::Union{Bool,String,AbstractMatrix} = "auto"
end
MMI.fitted_params(model::OrthogonalMatchingPursuitRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing)
    )
descr(OrthogonalMatchingPursuitRegressor, "Orthogonal Matching Pursuit (OMP) model.")

# =============================================================================
const OrthogonalMatchingPursuitCVRegressor_ = sklm(:OrthogonalMatchingPursuitCV)
@sk_reg mutable struct OrthogonalMatchingPursuitCVRegressor <: MMI.Deterministic
    copy::Bool            = true
    fit_intercept::Bool   = true
    normalize::Bool       = false
    max_iter::Option{Int} = nothing::(_ === nothing||_ > 0)
    cv::Any               = 5
    n_jobs::Option{Int}   = 1
    verbose::Union{Bool,Int} = false
end
MMI.fitted_params(model::OrthogonalMatchingPursuitCVRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    n_nonzero_coefs = fitresult.n_nonzero_coefs_
    )
descr(OrthogonalMatchingPursuitCVRegressor, "Orthogonal Matching Pursuit (OMP) model $CV.")

# =============================================================================
const PassiveAggressiveRegressor_ = sklm(:PassiveAggressiveRegressor)
@sk_reg mutable struct PassiveAggressiveRegressor <: MMI.Deterministic
    C::Float64                   = 1.0::(_ > 0)
    fit_intercept::Bool          = true
    max_iter::Int                = 1_000::(_ > 0)
    tol::Float64                 = 1e-4::(_ > 0)
    early_stopping::Bool         = false
    validation_fraction::Float64 = 0.1::(0 < _ < 1)
    n_iter_no_change::Int        = 5::(_ > 0)
    shuffle::Bool                = true
    verbose::Union{Bool,Int}     = 0
    loss::String                 = "epsilon_insensitive"::(_ in ("epsilon_insensitive","squared_epsilon_insensitive"))
    epsilon::Float64             = 0.1::(_ > 0)
    random_state::Any            = nothing
    warm_start::Bool             = false
    average::Union{Bool,Int}     = false
end
MMI.fitted_params(model::PassiveAggressiveRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing)
    )
descr(PassiveAggressiveRegressor, "Passive Aggressive Regressor.")

# =============================================================================
const RANSACRegressor_ = sklm(:RANSACRegressor)
@sk_reg mutable struct RANSACRegressor <: MMI.Deterministic
    base_estimator::Any         = nothing
    min_samples::Union{Int,Float64}     = 5::(_ isa Int ? _ ≥ 1 : (0 ≤ _ ≤ 1))
    residual_threshold::Option{Float64} = nothing
    is_data_valid::Any          = nothing
    is_model_valid::Any         = nothing
    max_trials::Int             = 100::(_ > 0)
    max_skips::Int              = typemax(Int)::(_ > 0)
    stop_n_inliers::Int         = typemax(Int)::(_ > 0)
    stop_score::Float64         = Inf::(_ > 0)
    stop_probability::Float64   = 0.99::(0 ≤ _ ≤ 1.0)
    loss::Union{Function,String}= "absolute_loss"::((_ isa Function) || _ in ("absolute_loss","squared_loss"))
    random_state::Any           = nothing
end
MMI.fitted_params(m::RANSACRegressor, (f, _, _)) = (
    estimator             = f.estimator_,
    n_trials              = f.n_trials_,
    inlier_mask           = f.inlier_mask_,
    n_skips_no_inliers    = f.n_skips_no_inliers_,
    n_skips_invalid_data  = f.n_skips_invalid_data_,
    n_skips_invalid_model = f.n_skips_invalid_model_
    )
descr(RANSACRegressor, "RANSAC regressor.")

# =============================================================================
const RidgeRegressor_ = sklm(:Ridge)
@sk_reg mutable struct RidgeRegressor <: MMI.Deterministic
    alpha::Union{Float64,Vector{Float64}} = 1.0::(all(_ .> 0))
    fit_intercept::Bool = true
    normalize::Bool     = false
    copy_X::Bool        = true
    max_iter::Int       = 1_000::(_ > 0)
    tol::Float64        = 1e-4::(_ > 0)
    solver::String      = "auto"::(_ in ("auto","svd","cholesky","lsqr","sparse_cg","sag","saga"))
    random_state::Any   = nothing
end
MMI.fitted_params(model::RidgeRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing)
    )
descr(RidgeRegressor, "Ridge regression.")

# =============================================================================
const RidgeCVRegressor_ = sklm(:RidgeCV)
@sk_reg mutable struct RidgeCVRegressor <: MMI.Deterministic
    alphas::Any              = (0.1, 1.0, 10.0)::(all(_ .> 0))
    fit_intercept::Bool      = true
    normalize::Bool          = false
    scoring::Any             = nothing
    cv::Any                  = 5
    gcv_mode::Option{String} = nothing::(_ === nothing || _ in ("auto","svd","eigen"))
    store_cv_values::Bool    = false
end
MMI.fitted_params(model::RidgeCVRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    alpha     = fitresult.alpha_,
    cv_values = model.store_cv_values ? fitresult.cv_values_ : nothing
    )
descr(RidgeCVRegressor, "Ridge regression $CV.")

# =============================================================================
const SGDRegressor_ = sklm(:SGDRegressor)
@sk_reg mutable struct SGDRegressor <: MMI.Deterministic
    loss::String             = "squared_loss"::(_ in ("squared_loss","huber","epsilon_insensitive","squared_epsilon_insensitive"))
    penalty::String          = "l2"::(_ in ("none","l2","l1","elasticnet"))
    alpha::Float64           = 1e-4::(_ > 0)
    l1_ratio::Float64        = 0.15::(_ > 0)
    fit_intercept::Bool      = true
    max_iter::Int            = 1_000::(_ > 0)
    tol::Float64             = 1e-3::(_ > 0)
    shuffle::Bool            = true
    verbose::Union{Int,Bool} = 0
    epsilon::Float64         = 0.1
    random_state::Any        = nothing
    learning_rate::String    = "invscaling"::(_ in ("constant","optimal","invscaling","adaptive"))
    eta0::Float64            = 0.01::(_ > 0)
    power_t::Float64         = 0.25::(_ > 0)
    early_stopping::Bool     = false
    validation_fraction::Float64 = 0.1::(0 < _ < 1)
    n_iter_no_change::Int    = 5::(_ > 0)
    warm_start::Bool         = false
    average::Union{Int,Bool} = false
end
MMI.fitted_params(model::SGDRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    average_coef      = model.average ? fitresult.average_coef_ : nothing,
    average_intercept = model.average ? ifelse(model.fit_intercept, fitresult.average_intercept_, nothing) : nothing
    )
descr(SGDRegressor, "Stochastic gradient descent-based regressor.")

# =============================================================================
const TheilSenRegressor_ = sklm(:TheilSenRegressor)
@sk_reg mutable struct TheilSenRegressor <: MMI.Deterministic
    fit_intercept::Bool = true
    copy_X::Bool        = true
    max_subpopulation::Int    = 10_000::(_ > 0)
    n_subsamples::Option{Int} = nothing::(_ === nothing||_ > 0)
    max_iter::Int       = 300::(_ > 0)
    tol::Float64        = 1e-3::(_ > 0)
    random_state::Any   = nothing
    n_jobs::Option{Int} = nothing
    verbose::Bool       = false
end
MMI.fitted_params(model::TheilSenRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    breakdown       = fitresult.breakdown_,
    n_subpopulation = fitresult.n_subpopulation_
    )
descr(TheilSenRegressor, "Theil-Sen regressor.")

# Metadata for Continuous -> Vector{Continuous}
const SKL_REGS_SINGLE = Union{
        Type{<:ARDRegressor},
        Type{<:BayesianRidgeRegressor},
        Type{<:ElasticNetRegressor},
        Type{<:ElasticNetCVRegressor},
        Type{<:HuberRegressor},
        Type{<:LarsRegressor},
        Type{<:LarsCVRegressor},
        Type{<:LassoRegressor},
        Type{<:LassoCVRegressor},
        Type{<:LassoLarsRegressor},
        Type{<:LassoLarsCVRegressor},
        Type{<:LassoLarsICRegressor},
        Type{<:LinearRegressor},
        Type{<:OrthogonalMatchingPursuitRegressor},
        Type{<:OrthogonalMatchingPursuitCVRegressor},
        Type{<:PassiveAggressiveRegressor},
        Type{<:RANSACRegressor},
        Type{<:RidgeRegressor},
        Type{<:RidgeCVRegressor},
        Type{<:SGDRegressor},
        Type{<:TheilSenRegressor}
    }

MMI.input_scitype(::SKL_REGS_SINGLE)  = Table(Continuous)
MMI.target_scitype(::SKL_REGS_SINGLE) = AbstractVector{Continuous}


text:
# AUTO GENERATED FILE - DO NOT EDIT

export slider

"""
    slider(;kwargs...)

A Slider component.
Slider
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `name` (String; optional): Form item name
- `inputType` (String; optional): The form element type
- `ariaLabelInput` (String; optional): The ARIA label for the <input>
- `labelText` (String; required): Label text
- `disabled` (Bool; optional): Disabled
- `light` (Bool; optional): Light variant
- `hideTextInput` (Bool; optional): Without text input
- `value` (Real; required): The value
- `min` (Real; required): The minimum value
- `max` (Real; required): The maximum value
- `step` (Real; required): The step
- `stepMultiplier` (Real; optional): The step factor for Shift+arrow keys
"""
function slider(; kwargs...)
        available_props = Symbol[:id, :name, :inputType, :ariaLabelInput, :labelText, :disabled, :light, :hideTextInput, :value, :min, :max, :step, :stepMultiplier]
        wild_props = Symbol[]
        return Component("slider", "Slider", "dash_carbon_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: Matheus-Rangel/dash-carbon-components

max_issues_repo_name: Matheus-Rangel/dash-carbon-components

max_forks_repo_name: Matheus-Rangel/dash-carbon-components


text: Dict{String,Any}("bestdayever" => Dict{String,Any}("value" => "2017-06-06T12:34:56-05:00","type" => "datetime"))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
function circle_len1(θ; n_pts=1000000)
    xmin, xmax = cos(θ), 1
    ts = sort(rand(n_pts))
    xs = @. (xmax-xmin)/(1. - 0.) * ts + xmin
    ys = @. sqrt(1 - xs * xs)
    d = zero(eltype(xs))
    for i in 1:(length(xs)-1)
        d += sqrt((xs[i+1]-xs[i])^2+(ys[i+1]-ys[i])^2)
    end
    println(d)
    println(θ)
end

function circle_len2(θ; n_pts=1000000)
    xmin, xmax = cos(θ), 1
    t = sort(rand(n_pts))
    x = @. (xmax-xmin)/(1. - 0.) * t + xmin
    y = @. sqrt(1 - x * x)
    ds = @. sqrt((x[2:end] - x[1:end-1])^2+(y[2:end] - y[1:end-1])^2)
    d=sum(ds)
    println(d)
    println(θ)
end

function circle_len3(θ; n_pts=1000000)
    ts = θ * sort(rand(n_pts))
    xs = sin.(ts)
    ys = cos.(ts)
    d = zero(eltype(xs))
    for i in 1:(length(xs)-1)
        d += sqrt((xs[i+1]-xs[i])^2+(ys[i+1]-ys[i])^2)
    end
    println(d)
    println(θ)
end

circle_len3(2π)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
@cdef :PyLong_FromLongLong PyPtr (Clonglong,)
@cdef :PyLong_FromUnsignedLongLong PyPtr (Culonglong,)
@cdef :PyLong_FromString PyPtr (Cstring, Ptr{Cvoid}, Cint)
@cdef :PyLong_AsLongLong Clonglong (PyPtr,)
@cdef :PyLong_AsUnsignedLongLong Culonglong (PyPtr,)

const PyLong_Type__ref = Ref(PyNULL)
PyLong_Type() = pyglobal(PyLong_Type__ref, :PyLong_Type)

PyLong_Check(o) = Py_TypeCheckFast(o, Py_TPFLAGS_LONG_SUBCLASS)

PyLong_CheckExact(o) = Py_TypeCheckExact(o, PyLong_Type())

PyLong_From(
    x::Union{Bool,Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128,BigInt},
) =
    if x isa Signed && typemin(Clonglong) ≤ x ≤ typemax(Clonglong)
        PyLong_FromLongLong(x)
    elseif typemin(Culonglong) ≤ x ≤ typemax(Culonglong)
        PyLong_FromUnsignedLongLong(x)
    else
        PyLong_FromString(string(x), C_NULL, 10)
    end

PyLong_From(x::Integer) = begin
    y = tryconvert(BigInt, x)
    y === PYERR() && return PyNULL
    y === NOTIMPLEMENTED() && (
        PyErr_SetString(
            PyExc_NotImplementedError(),
            "Cannot convert this Julia '$(typeof(x))' to a Python 'int'",
        );
        return PyNULL
    )
    PyLong_From(y::BigInt)
end

# "Longable" means an 'int' or anything with an '__int__' method.
PyLongable_TryConvertRule_integer(o, ::Type{S}) where {S<:Integer} = begin
    # first try to convert to Clonglong (or Culonglong if unsigned)
    x = S <: Unsigned ? PyLong_AsUnsignedLongLong(o) : PyLong_AsLongLong(o)
    if !ism1(x) || !PyErr_IsSet()
        # success
        return putresult(tryconvert(S, x))
    elseif PyErr_IsSet(PyExc_OverflowError())
        # overflows Clonglong or Culonglong
        PyErr_Clear()
        if S in (
               Bool,
               Int8,
               Int16,
               Int32,
               Int64,
               Int128,
               UInt8,
               UInt16,
               UInt32,
               UInt64,
               UInt128,
           ) &&
           typemin(typeof(x)) ≤ typemin(S) &&
           typemax(S) ≤ typemax(typeof(x))
            # definitely overflows S, give up now
            return 0
        else
            # try converting to String then BigInt then S
            so = PyObject_Str(o)
            isnull(so) && return -1
            s = PyUnicode_AsString(so)
            Py_DecRef(so)
            isempty(s) && PyErr_IsSet() && return -1
            y = tryparse(BigInt, s)
            y === nothing && (
                PyErr_SetString(
                    PyExc_ValueError(),
                    "Cannot convert this '$(PyType_Name(Py_Type(o)))' to a Julia 'BigInt' because its string representation cannot be parsed as an integer",
                );
                return -1
            )
            return putresult(tryconvert(S, y))
        end
    else
        # other error
        return -1
    end
end

PyLongable_TryConvertRule_tryconvert(o, ::Type{S}) where {S} = begin
    r = PyLongable_TryConvertRule_integer(o, Integer)
    r == 1 ? putresult(tryconvert(S, takeresult(Integer))) : r
end


text:
module Alpha

include("main.jl")

if abspath(PROGRAM_FILE) == @__FILE__
    alphamain()
end

end


text:
# ## Linear Regression

# ### Added snippet used as a reference for all models

# This model is based on the TuringTutorial example [LinearRegression](https://github.com/TuringLang/TuringTutorials/blob/csp/linear/LinearRegression.ipynb) by Cameron Pfiffer.

# Turing is powerful when applied to complex hierarchical models, but it can also be put to task at common statistical procedures, like linear regression. This tutorial covers how to implement a linear regression model in Turing.

# We begin by importing all the necessary libraries.

using StatisticalRethinking, CmdStan, GLM
#gr(size=(600,600))

ProjDir = rel_path("..", "scripts", "00")
cd(ProjDir)

# Import the dataset.

howell1 = CSV.read(rel_path("..", "data", "Howell1.csv"), delim=';')
df = convert(DataFrame, howell1);

# Use only adults

data = filter(row -> row[:age] >= 18, df)

# Show the first six rows of the dataset.

first(data, 6)

# The next step is to get our data ready for testing. We'll split the mtcars dataset into two subsets, one for training our model and one for evaluating our model. Then, we separate the labels we want to learn (MPG, in this case) and standardize the datasets by subtracting each column's means and dividing by the standard deviation of that column.

# The resulting data is not very familiar looking, but this standardization process helps the sampler converge far easier. We also create a function called unstandardize, which returns the standardized values to their original form. We will use this function later on when we make predictions.

# Split our dataset 70%/30% into training/test sets.

n = size(data, 1)
test_ind = sample(1:n, Int(floor(0.3*n)), replace=false);
train_ind = [(i) for i=1:n if !(i in test_ind)];
test = data[test_ind, :];
train = data[train_ind, :];

# Save dataframe versions of our dataset.

train_cut = DataFrame(train)
test_cut = DataFrame(test)

# Create our labels. These are the values we are trying to predict.

train_label = train[:, :height]
test_label = test[:, :height]

# Get the list of columns to keep.

remove_names = filter(x->!in(x, [:height, :age, :male]), names(data))

# Filter the test and train sets.

train = Matrix(train[:, remove_names]);
test = Matrix(test[:, remove_names]);

# A handy helper function to rescale our dataset.

function standardize(x)
    return (x .- mean(x, dims=1)) ./ std(x, dims=1), x
end

# Another helper function to unstandardize our datasets.

function unstandardize(x, orig)
    return x .* std(orig, dims=1) .+ mean(orig, dims=1)
end

# Standardize our dataset.

(train, train_orig) = standardize(train)
(test, test_orig) = standardize(test)
(train_label, train_l_orig) = standardize(train_label)
(test_label, test_l_orig) = standardize(test_label);

# Design matrix

dmat = [ones(size(train, 1)) train]

# Bayesian linear regression.

lrmodel = "
data {
  int N; //the number of observations
  int K; //the number of columns in the model matrix
  real y[N]; //the response
  matrix[N,K] X; //the model matrix
}
parameters {
  vector[K] beta; //the regression parameters
  real sigma; //the standard deviation
}
transformed parameters {
  vector[N] linpred;
  linpred <- X*beta;
}
model {
  beta[1] ~ cauchy(0,10); // prior for the intercept following Gelman 2008

  for(i in 2:K)
   beta[i] ~ cauchy(0,2.5); // prior for the slopes following Gelman 2008

  y ~ normal(linpred,sigma);
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

stanmodel = Stanmodel(name="linear_regression",
  model=lrmodel, output_format=:mcmcchains);

# Input data for cmdstan

lrdata = Dict("N" => size(train, 1), "K" => size(dmat, 2), "y" => train_label, "X" => dmat);

# Sample using cmdstan

rc, chain, cnames = stan(stanmodel, lrdata, ProjDir, diagnostics=false,
  summary=false, CmdStanDir=CMDSTAN_HOME);

# Convert to a  Chain object

chns = set_section(chain, Dict(
    :parameters => ["beta.1", "beta.2", "sigma"],
    :linpred => ["linpred.$i" for i in 1:247],
    :internals => ["lp__", "accept_stat__", "stepsize__", "treedepth__",
      "n_leapfrog__", "divergent__", "energy__"]
  )
)

# Describe the chains.

MCMCChains.describe(chns)

# Perform multivariate OLS.

ols = lm(@formula(height ~ weight), train_cut)

# Store our predictions in the original dataframe.
train_cut.OLSPrediction = predict(ols);
test_cut.OLSPrediction = predict(ols, test_cut);

# Make a prediction given an input vector.

function prediction(chn, x)
    α = Array(chn[Symbol("beta.1")]);
    β = Array(chn[Symbol("beta.2")]);
    return  mean(α) .+ x .* mean(β)
end

# Calculate the predictions for the training and testing sets.

train_cut.BayesPredictions = unstandardize(prediction(chns, train), train_l_orig)[:,1];
test_cut.BayesPredictions = unstandardize(prediction(chns, test), test_l_orig)[:,1];

# Show the first side rows of the modified dataframe.

remove_names = filter(x->!in(x, [:age, :male]), names(test_cut));
test_cut = test_cut[remove_names];
first(test_cut, 6)

bayes_loss1 = sum((train_cut.BayesPredictions - train_cut.height).^2);
ols_loss1 = sum((train_cut.OLSPrediction - train_cut.height).^2);

bayes_loss2 = sum((test_cut.BayesPredictions - test_cut.height).^2);
ols_loss2 = sum((test_cut.OLSPrediction - test_cut.height).^2);

println("\nTraining set:")
println("  Bayes loss: $bayes_loss1")
println("  OLS loss: $ols_loss1")

println("Test set:")
println("  Bayes loss: $bayes_loss2")
println("  OLS loss: $ols_loss2")

# Plot the chains.

#plot(chain)


text:
using Poirot.Abstract, IRTools, Test
using Poirot.Abstract: Const, Partial, Inference, return_type, @trace, exprtype, returntype
using IRTools: var, returnvalue

ir = @code_ir identity(1)
@test return_type(ir, Nothing, Int) == Int

add(a, b) = a+b

ir = @code_ir add(1, 2)

@test return_type(ir, Nothing, Const(1.0), Const(2)) == Const(3.0)

@test return_type(ir, Nothing, Int, Const(2.0)) == Float64

g = 2
addg(x) = x+g
ir = @code_ir addg(2)
@test return_type(ir, Nothing, Const(2)) == Const(4)
@test return_type(ir, Nothing, Float64) == Float64

f() = rand() < 0.5 ? "foo" : "bar"

ir = @code_ir f()
@test return_type(ir, Nothing) == String

function pow(x, n)
  r = 1
  while n > 0
    n -= 1
    r *= x
  end
  return r
end

ir = @code_ir pow(2, 3)
@test return_type(ir, Nothing, Int64, Const(5)) == Int
@test return_type(ir, Nothing, Float64, Int) == Union{Float64, Int}
@test return_type(ir, Nothing, Float64, Const(5)) == Float64
@test return_type(ir, Nothing, Const(2), Const(3)) == Const(8)

foo(a, b) = a + b
bar(a, b) = foo(a, b)

ir = @code_ir bar(1, 1)
@test return_type(ir, Nothing, Int, Int) == Int

fact(n) = n == 0 ? 1 : n*fact(n-1)

ir = @code_ir fact(1)
@test return_type(ir, Nothing, Int) == Int

# Tracing

tr = @trace f()
@test returntype(tr) == String

tr = @trace pow(Int, 3)
@test length(tr.blocks) == 1
@test returntype(tr) == Int

tr = @trace pow(2, 3)
@test length(tr.blocks) == 1
@test returntype(tr) == Const(8)

tr = @trace pow(2, Int)
@test returntype(tr) == Int

tr = @trace pow(2.0, Int)
@test returntype(tr) == Union{Float64,Int}

tr = @trace pow(1, Int)
@test returntype(tr) == Const(1)

tr = @trace bar(Int, Int)
@test returntype(tr) == Int

function foo(x)
  r = Ref{Any}()
  r[] = x
  return r[]
end

tr = @trace foo(1)
@test returntype(tr) == Const(1)

tr = @trace foo(Int)
@test returntype(tr) == Int

function foo(x)
  r = Ref{Any}()
  r[] = x
  if rand(Bool)
    r[] = 1
  end
  return r[]
end

tr = @trace foo(1)
@test returntype(tr) == Const(1)

tr = @trace foo(2)
@test returntype(tr) == Int

struct Foo
  x
end

function foo(x)
  r = Foo(x)
  return r.x
end

tr = @trace foo(2)
@test returntype(tr) == Const(2)

foo(x) = Foo(x)

tr = @trace foo(1)
@test returntype(tr) == Const(Foo(1))

function pow(x, n)
  env = Dict()
  env[:r] = 1
  env[:n] = n
  while env[:n] > 0
    env[:r] *= x
    env[:n] -= 1
  end
  return env[:r]
end

tr = @trace pow(2, 3)
@test returntype(tr) == Const(8)

tr = @trace pow(Int, Int)
@test returntype(tr) == Int


text:
using ThreadPinning
using Test
using Base.Threads: @threads, nthreads
using Random: shuffle

nthreads() ≥ 2 ||
    error("Can't run tests with single Julia thread! Forgot to set `JULIA_NUM_THREADS`?")

function check_compact_within_socket(cpuids)
    socket_cpuids = cpuids_per_socket()
    for s in 1:nsockets()
        cpuids_filtered = filter(x -> x in socket_cpuids[s], cpuids)
        if cpuids_filtered != socket_cpuids[s][1:length(cpuids_filtered)]
            return false
        end
    end
    return true
end

function check_compact_within_numa(cpuids)
    numa_cpuids = cpuids_per_numa()
    for s in 1:nnuma()
        cpuids_filtered = filter(x -> x in numa_cpuids[s], cpuids)
        if cpuids_filtered != numa_cpuids[s][1:length(cpuids_filtered)]
            return false
        end
    end
    return true
end

@testset "ThreadPinning.jl" begin
    @testset "Helper" begin
        @test ThreadPinning.interweave([1, 2, 3, 4], [5, 6, 7, 8]) ==
            [1, 5, 2, 6, 3, 7, 4, 8]
        @test ThreadPinning.interweave(1:4, 5:8) == [1, 5, 2, 6, 3, 7, 4, 8]
        @test ThreadPinning.interweave(1:4, 5:8, 9:12) ==
            [1, 5, 9, 2, 6, 10, 3, 7, 11, 4, 8, 12]
        # different size inputs
        @test_throws ArgumentError ThreadPinning.interweave([1, 2, 3, 4], [5, 6, 7, 8, 9])
    end

    @testset "gather_sysinfo_lscpu (NOCTUA2LOGIN)" begin
        sinfo = ThreadPinning.gather_sysinfo_lscpu(ThreadPinning.lscpu_NOCTUA2LOGIN)
        @test typeof(sinfo) == ThreadPinning.SysInfo
        @test sinfo.nsockets == 1
        @test sinfo.nnuma == 4
        @test sinfo.hyperthreading == true
        @test sinfo.cpuids == 0:127
        @test length(sinfo.cpuids_sockets) == 1
        @test sinfo.cpuids_sockets[1] == 0:127
        @test length(sinfo.cpuids_numa) == 4
        @test sinfo.cpuids_numa[1] == vcat(0:15, 64:79)
        @test sinfo.cpuids_numa[2] == vcat(16:31, 80:95)
        @test sinfo.cpuids_numa[3] == vcat(32:47, 96:111)
        @test sinfo.cpuids_numa[4] == vcat(48:63, 112:127)
        @test sinfo.ishyperthread == vcat(falses(64), trues(64))
    end

    @testset "gather_sysinfo_lscpu (NOCTUA2)" begin
        sinfo = ThreadPinning.gather_sysinfo_lscpu(ThreadPinning.lscpu_NOCTUA2)
        @test typeof(sinfo) == ThreadPinning.SysInfo
        @test sinfo.nsockets == 2
        @test sinfo.nnuma == 8
        @test sinfo.hyperthreading == false
        @test sinfo.cpuids == 0:127
        @test length(sinfo.cpuids_sockets) == 2
        @test sinfo.cpuids_sockets[1] == 0:63
        @test sinfo.cpuids_sockets[2] == 64:127
        @test length(sinfo.cpuids_numa) == 8
        @test sinfo.cpuids_numa[1] == 0:15
        @test sinfo.cpuids_numa[2] == 16:31
        @test sinfo.cpuids_numa[3] == 32:47
        @test sinfo.cpuids_numa[4] == 48:63
        @test sinfo.cpuids_numa[5] == 64:79
        @test sinfo.cpuids_numa[6] == 80:95
        @test sinfo.cpuids_numa[7] == 96:111
        @test sinfo.cpuids_numa[8] == 112:127
        @test sinfo.ishyperthread == falses(128)
    end

    @testset "gather_sysinfo_lscpu (NOCTUA1)" begin
        sinfo = ThreadPinning.gather_sysinfo_lscpu(ThreadPinning.lscpu_NOCTUA1)
        @test typeof(sinfo) == ThreadPinning.SysInfo
        @test sinfo.nsockets == 2
        @test sinfo.nnuma == 2
        @test sinfo.hyperthreading == false
        @test sinfo.cpuids == 0:39
        @test length(sinfo.cpuids_sockets) == 2
        @test sinfo.cpuids_sockets[1] == 0:19
        @test sinfo.cpuids_sockets[2] == 20:39
        @test length(sinfo.cpuids_numa) == 2
        @test sinfo.cpuids_numa[1] == 0:19
        @test sinfo.cpuids_numa[2] == 20:39
        @test sinfo.ishyperthread == falses(40)
    end

    @testset "gather_sysinfo_lscpu (FUGAKU)" begin
        sinfo = ThreadPinning.gather_sysinfo_lscpu(ThreadPinning.lscpu_FUGAKU)
        @test typeof(sinfo) == ThreadPinning.SysInfo
        @test sinfo.nsockets == 1
        @test sinfo.nnuma == 6
        @test sinfo.hyperthreading == false
        @test sinfo.cpuids == vcat([0, 1], 12:59)
        @test length(sinfo.cpuids_sockets) == 1
        @test sinfo.cpuids_sockets[1] == sinfo.cpuids
        @test length(sinfo.cpuids_numa) == 6
        @test sinfo.cpuids_numa[1] == [0]
        @test sinfo.cpuids_numa[2] == [1]
        @test sinfo.cpuids_numa[3] == 12:23
        @test sinfo.cpuids_numa[4] == 24:35
        @test sinfo.cpuids_numa[5] == 36:47
        @test sinfo.cpuids_numa[6] == 48:59
        @test sinfo.ishyperthread == falses(50)
    end

    @testset "Querying" begin
        pinthreads(:random)
        @test typeof(getcpuid()) == Int
        @test typeof(getcpuids()) == Vector{Int}
        @test getcpuids() == getcpuid.(1:Threads.nthreads())
        @test typeof(nsockets()) == Int
        @test nsockets() >= 1
        @test typeof(hyperthreading_is_enabled()) == Bool
        @test typeof(cpuids_per_socket()) == Vector{Vector{Int}}
        @test ishyperthread(0) == false
    end

    @testset "threadinfo()" begin
        @test isnothing(threadinfo())
        @test isnothing(threadinfo(; groupby=:numa))
    end

    @testset "Thread Pining (explicit)" begin
        cpuid_before = getcpuid()
        cpuid_new = cpuid_before != 1 ? 1 : 0
        @test pinthread(cpuid_new)
        @test getcpuid() == cpuid_new
        cpuids_new = shuffle(0:(nthreads() - 1))
        @test isnothing(pinthreads(cpuids_new))
        @test getcpuids() == cpuids_new
        cpuids_new = reverse(0:(nthreads() - 1))
        @test isnothing(pinthreads(cpuids_new))
        @test getcpuids() == cpuids_new

        rand_thread = rand(1:Threads.nthreads())
        for cpuid in rand(cpuids_all(), 5)
            @test isnothing(pinthread(rand_thread, cpuid))
            @test getcpuid(rand_thread) == cpuid
        end
    end

    @testset "Thread Pinning (compact)" begin
        @assert getcpuids() != 0:(nthreads() - 1)
        @test isnothing(pinthreads(:compact; nthreads=2))
        @test getcpuids()[1:2] == 0:1
        @test isnothing(pinthreads(:compact))
        @test getcpuids() == 0:(nthreads() - 1)
    end

    @testset "Thread Pinning (scatter)" begin
        @test isnothing(pinthreads(:scatter))
        cpuids_after = getcpuids()
        @test check_compact_within_socket(cpuids_after)
    end

    @testset "Thread Pinning (numa)" begin
        @test isnothing(pinthreads(:numa))
        cpuids_after = getcpuids()
        @test check_compact_within_numa(cpuids_after)
    end

    @testset "Core2CoreLatency" begin
        latencies = ThreadPinning.bench_core2core_latency()
        @test typeof(latencies) == Matrix{Float64}
        @test size(latencies) == (Sys.CPU_THREADS, Sys.CPU_THREADS)
    end
end


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module SeqUtils

using Nemo

export ModuleSeqUtils
export SeqShow, SeqPrint, SeqName, SeqNum, Println, ZArray, Nemofmpz

"""
Nemo is a library designed, developed and maintained by William Hart with the help of others. Many functions in our project are based on Nemo.
"""
const ModuleSeqUtils = ""

"""
The basic integer type of our project is 'Nemo.fmpz'. Think of a 'Nemo.fmpz' as a 'BigInt'. 'fmpz' stands for 'fast multiple precision zahl' (zahl=integer). The mathematical symbol for the ring of integers is the blackboard (double-struck) Z, also written ZZ. In reference to this ZZ(n) defines the integer n as a fmpz.
"""
function Nemofmpz(n::Int) return ZZ(n) end

"""
Return an array of type `fmpz` of length ``n`` preset with ``0``.
"""
ZArray(len::Int) = zeros(ZZ, len)

"""
Return an array of type fmpz and of size ``n`` filled by the values of the function ``f`` in the range `offset:n`.
"""
function ZArray(n::Int, f::Function, offset=1)
    n ≤ 0 && return fmpz[]
    fmpz[f(k) for k in offset:offset+n-1] # !
end

"""
Return the size of the array ``A``.
"""
SeqSize(A) = Base.length(A)

"""
Return the range of a SeqArray.
"""
SeqRange(A) = 1:SeqSize(A)

"""
Print the array ``A`` in the format ``n ↦ A[n]`` for n in the given range.
"""
function SeqShow(A::Array, R::AbstractRange, offset=0)
    for k in R
        if isassigned(A, k)
            println(k + offset, " ↦ ", A[k])
        else
            println(k + offset, " ↦ ", "undef")
        end
    end
end

"""
Print the array ``A`` in the format ``n ↦ A[n]`` for n in the range first:last.
"""
SeqShow(A, first::Int, last::Int) = SeqShow(A, first:last)

"""
Print the array ``A`` in the format ``n ↦ A[n]``.
"""
SeqShow(A::Array, offset=0) = SeqShow(A, 1:length(A), offset - 1)

function print_without_type(io, v::AbstractVector)
    print(io, "[")
    for (i, el) in enumerate(v)
        i > 1 && print(io, ", ")
        print(io, el)
    end
    println(io, "]")
end

"""
Print the array without typeinfo.
"""
Println(v::AbstractVector) = print_without_type(IOContext(stdout), v)

"""
Print the array with or without typeinfo.
"""
function SeqPrint(v::AbstractVector, typeinfo = false)
    typeinfo ? println(v) : Println(v)
end

"""
Valid prefixes to the numerical part of the OEIS A-numbers.

* C => Channel
* F => Filter (all below n)
* G => Generating function
* I => Iteration
* L => List (array based)
* M => Matrix
* R => RealFunction
* S => Staircase (iteration)
* T => Triangle (iteration)
* TA => Triangle (triangular array)
* TL => Triangle (flat-list array)
* V => Value
* is => is a (predicate), boolean
"""
const ValidPrefixes = ['C', 'F', 'G', 'I', 'L', 'M', 'R', 'S', 'T', 'V']

"""
Return the name of a OEIS sequence given a similar named function as a string.
"""
function SeqName(fun)
    aname = string(fun)
    for X in ValidPrefixes
        aname = replace(aname, X => 'A')
    end

    if !occursin(r"^A[0-9]{6}$", aname)
        fullname = split(aname, ".")
        aname = String(fullname[end])
        if !occursin(r"^A[0-9]{6}$", aname)
            @error(aname, " is not a valid A-name!")
            return
        end
    end
    aname
end

"""
Return the A-number of a OEIS sequence given a similar named function as an integer.
"""
function SeqNum(seq)
    aname = SeqName(seq)
    parse(Int, aname[2:end])
end

#START-TEST-########################################################

using Test

function test()
    #@testset "SeqUtils" begin
    #end
end

function demo()
    A = fmpz[1, 2, 3, 4, 5]

    println()
    SeqShow(A)
    println()

    SeqShow(A, -4)
    println()

    A |> Println
end

function perf()
end

function main()
    test()
    demo()
    perf()
end

main()

end # module


text:
using Sigma
using Lens

srand(1234)

# Paramter Estimation
λ = uniform(0,2)
x = mvexponential(λ, 3)
observations = x == [0.083, 0.55, 2.37]
rand(λ, observations)
rand(x, observations)

m = 3
n = 3
A = Sigma.mvuniform(-1,1,m,n);

# Sample matrix whose sum is 3
As = rand(A, sum(A) == 3)
sum(As)

# Sample Orthogonal Matrix
orthogonal_A = rand(A, (A'*A == eye(m,n)))
orthogonal_A * orthogonal_A'
map(round, orthogonal_A * orthogonal_A')

# trace
B = mvexponential(0.5, m, n);
C = mvlogistic(1, 1, m, n);

Bs, Cs = rand((B,C), (trace(B) == trace(C)), 1)
trace(Bs[1])
trace(Cs[1])

# Simple
A = uniform(0,1)
B = exponential(.5)
rand(A, A^2/sin(A) > 0.45 * atan(B))
rand(A, A < -1)

## Non Negative Matrix Factorisation
## ==================================
dist(A,B; epsilon = 1.0) = sum([(A[i] - B[i])^2 for i = 1:length(A)]) < epsilon

function nmf(V, W, H)
  dist(V, W*H)
end

W = mvuniform(0.0,10.0,4,2)
H = mvuniform(0.0,10.0,2,4)
V = rand(4,4)

nmf(V, W, H)


text:
"""
    add_to_visicon!(actr, vo; stuff=false)

Adds a visual object to the visicon. If `stuff` is set to true, the visual object
is added to the visual location buffer.

# Arguments

- `actr`: an ACT-R model object
- `vo`: visual object
- `stuff`: buffer stuffing if true
"""
function add_to_visicon!(actr::AbstractACTR, vo; stuff=false)
    push!(actr.visicon, deepcopy(vo))
    if stuff
       chunk = vo_to_chunk(actr, vo)
       add_to_buffer!(actr.visual_location, chunk)
    end
    return nothing
end

function add_to_visicon!(models, vo; stuff=false)
    for model in models
        add_to_visicon!(model, vo; stuff)
    end
    return nothing
end

"""
    clear_visicon!(actr)

Clear all visual objects in visicon
"""
function clear_visicon!(visicon)
    empty!(visicon)
end

clear_visicon!(actr::AbstractACTR) = clear_visicon!(actr.visicon)

"""
    remove_visual_object!(actr::AbstractACTR, vo)

Removes object from visicon.

# Arguments

- `actr`: an ACT-R model object
- `vo`: a visual object
"""
remove_visual_object!(actr::AbstractACTR, vo) = remove_visual_object!(actr.visual_location.visicon, vo)

"""
    remove_visual_object!(visicon, vo)

Removes object from visicon.

# Arguments

- `visicon`: a vector of visual objects
- `vo`: a visual object
"""
function remove_visual_object!(visicon, vo)
    filter!(x->x != vo, visicon)
end

"""
    vo_to_chunk(vo=VisualObject())

# Arguments

Converts visible object to a chunk with color and text slots.
"""
function vo_to_chunk(vo=VisualObject())
     return Chunk(;color=vo.color, text=vo.text, x=vo.x, y=vo.y)
end

"""
    vo_to_chunk(actr, vo)

Converts visible object to a chunk with color and text slots, and sets time created to current time.

# Arguments

- `actr`: an ACT-R model object
- `vo`: visual object
"""
function vo_to_chunk(actr, vo)
    time_created = get_time(actr)
    return Chunk(;time_created, color=vo.color, text=vo.text, x=vo.x, y=vo.y)
end

function move_vo!(actr, x, y)
    vo = actr.visicon[1]
    vo.x = x
    vo.y = y
end

"""
    attending!(actr, chunk, args...; kwargs...)

Sets visual module as busy and registers a new event to attend to a `chunk`
created by a visual object

# Arguments

- `actr`: an ACT-R model object
- `chunk`: a memory chunk
"""
function attending!(actr, chunk)
    actr.visual.state.busy = true
    description = "Attend"
    type = "model"
    id = get_name(actr)
    tΔ = rnd_time(.085)
    register!(actr, attend!, after, tΔ , actr, chunk; id, description, type)
end

function attending!(actr, chunk, task)
    actr.visual.state.busy = true
    description = "Attend"
    type = "model"
    id = get_name(actr)
    tΔ = rnd_time(.085)
    register!(actr, attend!, after, tΔ , actr, chunk, task; id, description, type)
end

"""
    attend!(actr, chunk, args...; kwargs...)

Completes an attention shift by adding a `chunk` to the visual buffer and setting
states to busy = false and empty = false.

# Arguments

- `actr`: an ACT-R model object
- `chunk`: a memory chunk
"""
function attend!(actr, chunk)
    slots = chunk.slots
    actr.visual.focus = [slots.x,slots.y]
    actr.visual.state.busy = false
    actr.visual.state.empty = false
    add_to_buffer!(actr.visual, chunk)
    return nothing
end

function attend!(actr, chunk, task)
    slots = chunk.slots
    actr.visual.focus = [slots.x,slots.y]
    actr.visual.state.busy = false
    actr.visual.state.empty = false
    add_to_buffer!(actr.visual, chunk)
    task.visible ? repaint!(task, actr) : nothing
    return nothing
end


text:
using Pkg.Artifacts
using Random
using Test

using DataFrames
using CSV

using LearningTrees

RNGSEED = 53
Random.seed!(RNGSEED)

loadiris() = begin
    artifact_toml = joinpath(@__DIR__, "Artifacts.toml")
    iris_hash = artifact_hash("iris", artifact_toml)
    if iris_hash === nothing || !artifact_exists(iris_hash)
        iris_hash = create_artifact() do artifact_dir
            iris_url_base = "https://archive.ics.uci.edu/ml/machine-learning-databases/iris"
            download("$(iris_url_base)/iris.data", joinpath(artifact_dir, "iris.data"))
            download("$(iris_url_base)/bezdekIris.data", joinpath(artifact_dir, "bezdekIris.csv"))
            download("$(iris_url_base)/iris.names", joinpath(artifact_dir, "iris.names"))
        end
        bind_artifact!(artifact_toml, "iris", iris_hash)
    end
    artifact_path(iris_hash)
end

testiris() = begin
    datasetpath = loadiris()
    names = [:sepal_length, :sepal_width, :petal_length, :petal_width, :class]
    df = CSV.File(joinpath(datasetpath, "iris.data"); header = names) |> DataFrame
    topredict = :class

    dftrain, dftest = traintestsplit(df, 0.33)

    c = true
    for i in 2:(length(propertynames(df))-1)
        clf = DecisionTree(;maxdepth = i, criterion = gini)
        fit!(clf, dftrain, topredict)
        score = accuracyscore(dftest[:, topredict], predict(clf, dftest))
        c = c && score >= 0.9
    end

    for i in [2,4,8,16,32,64]
        clf = RandomForest(;ntrees = i, criterion = gini)
        fit!(clf, dftrain, topredict)
        score = accuracyscore(dftest[:, topredict], predict(clf, dftest))
        c = c && score >= 0.9
    end

    c
end

@test testiris()