text:
struct FeatureDataset{T, D} <: ImageDataset{T, D}
    X
    Xt
end

n_display(::FeatureDataset) = 64

function FeatureDataset(
    n_data::Int,
    features::Matrix{T1};
    seed::Int=1,
    test_ratio=1/6,
    n_test::Int=ratio2num(n_data, test_ratio),
    alpha::T2=0f0,
    is_link::Bool=false,
    rng=MersenneTwister(seed),
) where {T1, T2}
    D, n_features = size(features)
    X = features * rand(rng, Bool, n_features, n_data)
    X = preprocess(rng, X, false, alpha, is_link)
    Xt = features * rand(rng, Bool, n_features, n_test)
    Xt = preprocess(rng, Xt, false, alpha, is_link)
    @info "Oh you just get the feature dataset" n_data=n_data n_test=n_test alpha=alpha is_link=is_link
    return FeatureDataset{Val{is_link}, D}(X, Xt)
end

"""
    get_features_griffiths2011()

Generate the same features for the synthesised dataset used in (Griffiths and Ghahramani, 2011).

Ref: http://www.jmlr.org/papers/v12/griffiths11a.html
"""
function get_features_griffiths2011()
    features = []

    push!(features, vec(Int[0 1 0 0 0 0;
                            1 1 1 0 0 0;
                            0 1 0 0 0 0;
                            0 0 0 0 0 0;
                            0 0 0 0 0 0;
                            0 0 0 0 0 0]))

    push!(features, vec(Int[0 0 0 1 1 1;
                            0 0 0 1 0 1;
                            0 0 0 1 1 1;
                            0 0 0 0 0 0;
                            0 0 0 0 0 0;
                            0 0 0 0 0 0]))

    push!(features, vec(Int[0 0 0 0 0 0;
                            0 0 0 0 0 0;
                            0 0 0 0 0 0;
                            1 0 0 0 0 0;
                            1 1 0 0 0 0;
                            1 1 1 0 0 0]))

    push!(features, vec(Int[0 0 0 0 0 0;
                            0 0 0 0 0 0;
                            0 0 0 0 0 0;
                            0 0 0 1 1 1;
                            0 0 0 0 1 0;
                            0 0 0 0 1 0]))

    features = hcat(features...)

    return features
end

"""
    get_features_xu2019()

Generate 35 12x12 features, the same features for the synthesised dataset used in (Xu et al., 2019).

Ref: http://proceedings.mlr.press/v97/xu19e.html
"""
function get_features_xu2019()

    features = []

    push!(features, vec(Int[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 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 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]))

    push!(features, vec(Int[1 0 0 0 0 0  0 0 0 0 0 0;
                            1 0 0 0 0 0  0 0 0 0 0 0;
                            1 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 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 0 0 0 0 0]))

    push!(features, vec(Int[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 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 0 0 0 0 0;
                            0 0 0 0 0 0  0 0 0 0 0 0]))

    push!(features, vec(Int[0 0 0 0 1 0  0 0 0 0 0 0;
                            0 0 0 1 1 1  0 0 0 0 0 0;
                            0 0 0 0 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 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 0 0 0]))

    push!(features, vec(Int[0 0 0 1 0 1  0 0 0 0 0 0;
                            0 0 0 0 0 0  0 0 0 0 0 0;
                            0 0 0 1 0 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 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 0 0]))

    push!(features, vec(Int[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;
                            1 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 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]))

    push!(features, vec(Int[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;
                            1 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 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]))

    push!(features, vec(Int[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 0 0 0 0;
                            1 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 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]))

    push!(features, vec(Int[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 1 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 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 0 0 0 0  0 0 0 0 0 0;
                            0 0 0 0 0 0  0 0 0 0 0 0]))

    push!(features, vec(Int[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 1 0 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 0 0 0 0 0  0 0 0 0 0 0;
                            0 0 0 0 0 0  0 0 0 0 0 0]))

    push!(features, vec(Int[0 0 0 0 0 0  1 1 1 0 0 0;
                            0 0 0 0 0 0  1 1 0 0 0 0;
                            0 0 0 0 0 0  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 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 0]))

    push!(features, vec(Int[0 0 0 0 0 0  0 0 0 0 0 0;
                            0 0 0 0 0 0  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 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 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]))

    push!(features, vec(Int[0 0 0 0 0 0  0 0 0 1 1 1;
                            0 0 0 0 0 0  0 0 0 1 0 1;
                            0 0 0 0 0 0  0 0 0 1 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 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]))

    push!(features, vec(Int[0 0 0 0 0 0  0 0 0 0 0 0;
                            0 0 0 0 0 0  0 0 0 0 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 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 0 0 0  0 0 0 0 0 0]))

    push!(features, vec(Int[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 1 0 0;
                            0 0 0 0 0 0  0 0 0 1 1 1;
                            0 0 0 0 0 0  0 0 0 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  0 0 0 0 0 0;
                            0 0 0 0 0 0  0 0 0 0 0 0]))

    push!(features, vec(Int[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 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 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]))

    push!(features, vec(Int[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 0 0 0 0;
                            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 0 0 0 0;
                            0 0 0 0 0 0  0 0 0 0 0 0]))

    push!(features, vec(Int[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 1 0 0 0;
                            0 0 0 0 0 0  1 1 1 0 0 0;
                            0 0 0 0 0 0  0 0 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 0  0 0 0 0 0 0;
                            0 0 0 0 0 0  0 0 0 0 0 0]))

    push!(features, vec(Int[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  1 1 0 0 0 0;
                            0 0 0 0 0 0  0 0 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 0 0  0 0 0 0 0 0;
                            0 0 0 0 0 0  0 0 0 0 0 0]))

    push!(features, vec(Int[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 0 0 0 0;
                            0 0 0 0 0 0  0 0 0 0 0 0;

                            0 1 0 0 0 0  0 0 0 0 0 0;
                            0 1 0 0 0 0  0 0 0 0 0 0;
                            1 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]))

    push!(features, vec(Int[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 0 0 0 0;
                            0 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 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 0]))

    push!(features, vec(Int[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 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;
                            1 1 1 0 0 0  0 0 0 0 0 0;
                            0 1 0 0 0 0  0 0 0 0 0 0;
                            0 1 0 0 0 0  0 0 0 0 0 0]))

    push!(features, vec(Int[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 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 0 0  0 0 0 0 0 0;
                            1 0 1 0 0 0  0 0 0 0 0 0;
                            1 0 1 0 0 0  0 0 0 0 0 0]))

    push!(features, vec(Int[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 0 0 0 0;
                            0 0 0 0 0 0  0 0 0 0 0 0;

                            0 0 0 1 0 0  0 0 0 0 0 0;
                            0 0 0 1 1 0  0 0 0 0 0 0;
                            0 0 0 1 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]))

    push!(features, vec(Int[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 0 0 0 0;
                            0 0 0 0 0 0  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 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]))

    push!(features, vec(Int[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 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 0 1  0 0 0 0 0 0;
                            0 0 0 0 1 1  0 0 0 0 0 0;
                            0 0 0 1 1 1  0 0 0 0 0 0]))

    push!(features, vec(Int[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 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 1 1 0  0 0 0 0 0 0;
                            0 0 0 1 0 0  0 0 0 0 0 0;
                            0 0 0 0 0 0  0 0 0 0 0 0]))

    push!(features, vec(Int[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 0 0 0 0;
                            0 0 0 0 0 0  0 0 0 0 0 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  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]))

    push!(features, vec(Int[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 0 0 0 0;
                            0 0 0 0 0 0  0 0 0 0 0 0;

                            0 0 0 0 0 0  0 0 1 0 0 0;
                            0 0 0 0 0 0  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 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]))

    push!(features, vec(Int[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 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 1 1 1;
                            0 0 0 0 0 0  0 0 0 0 0 1;
                            0 0 0 0 0 0  0 0 0 1 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]))

    push!(features, vec(Int[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 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 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]))

    push!(features, vec(Int[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 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 0 0  1 0 1 1 1 1;
                            0 0 0 0 0 0  1 0 1 0 0 1;
                            0 0 0 0 0 0  1 1 1 0 0 1]))

    push!(features, vec(Int[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 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 0 0  0 0 0 0 0 0;
                            0 0 0 0 0 0  0 0 0 1 1 0;
                            0 0 0 0 0 0  0 0 0 1 1 0]))

    push!(features, vec(Int[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 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 0 0  0 1 0 0 0 0;
                            0 0 0 0 0 0  0 1 0 0 0 0;
                            0 0 0 0 0 0  0 0 0 0 0 0]))

    features = hcat(features...)
    return features[:,[1,17,16,6,22,19,11,2,28,33,4,29,31,13,27,20,18,5,9,34,23,24,25,15,21,7,32,30,12,10,8,26,14,3]]
end


text:
using HTTP: request as rq

function url2xml(url::String)
    r = rq(:GET, url)
    xml_string = parse_string(String(r.body))
    return root(xml_string)
end

ext: jl


text:
"""
    ldu_update!(ldu::LDUPivoted, u::Vector, f::Real)

Given the LDU factorization `ldu: P' * A * P = L * D * L'` of a symmetric real or hermitian matrix `A`, update
`ldu` to obtain a factorization of `A + u * u' * f` without changing `L`.

This attempt may fail with a `SingularException` if the new pivot elements becomes zero.
"""
function ldu_update!(ldu::LDUPivoted, u::AbstractVector, f::Real)
    A = ldu.factors
    n = checksquare(A)
    length(u) == n || throw(ArgumentError("vector size $(length(u)) does not match matrix size $n"))
    u = ldu.P' * u

    for k = 1:n
        d1 = real(A[k,k])
        d2 = f
        u1 = u[k]
        dis = d1 + d2 * abs2(u1)
        A[k,k] = dis
        iszero(dis) && throw(SingularException(k))
        c = d1 / dis
        s = d2 * adjoint(u1) / dis
        for i = k+1:n
            ui = u[i] - A[k,i] * u1
            A[k,i] = A[k,i] * c + u[i] * s
            u[i] = ui
        end
        u[k] = 0
        f = d1 * d2 / dis
    end
    ldu
end

"""
    ldu_update!(ldu::LDUPivoted, u::Vector, v::Vector, f::Real)

Given the LDU factorization `ldu: P' * A * P = L * D * L'` of a symmetric real or hermitian matrix `A`, update
`ldu` to obtain a factorization of `A + (u * v' + v * u') * f` without changing `L`.

This attempt may fail with a `SingularException` if the new pivot elements becomes zero.
"""
function ldu_update!(ldu::LDUPivoted, u::AbstractVector, v::AbstractVector, f::Real)
    un = maximum(abs, u)
    vn = maximum(abs, v)
    (iszero(un) || iszero(vn) ) && return ldu
    g = sqrtapp(un / vn)
    u = u / g
    v = v * g
    f = f / 2
    upv = u + v
    umv = u - v
    ldu_update!(ldu, upv, f)
    ldu_update!(ldu, umv, -f)
end

splitexp(a::AbstractFloat) = frexp(a)
function splitexp(a::Union{Integer,Rational})
    q, x = frexp(float(a))
    x >= 0 ? a / 2^x : a * 2^(-x), x
end

"""
    sqrtapp(x)

Rational approximation of `sqrt(x)`.
"""
function sqrtapp(a::Real)
    q, x = splitexp(a)
    if isodd(x)
        q += q
        x -= 1
    end
    x ÷= 2
    c = (12 - (3 - q)^2) / 8
    # c = (q / c + c) / 2
    x >= 0 ? c * 2^x : c / 2^(-x)
end


text:
using Harmonie_namelists
using Test

# These are very basic test. Once we have JSONSchema files for namelists
# we could validate the namelists much better

# Check that the mechanism where we take variables from the ENVironment works
@testset "ENV" begin
    dicts = read_namelist.(["global", "canari"])
    totdict = merge_namelists(dicts)
    # ENV["ENSMBR"] is not defined so replace_env! throws an error
    @test_throws ErrorException replace_env!(totdict)

    # After setting ENV["ENSMBR"] the same call should succeed (returns nothing)
    Harmonie_namelists.MYENV["ENSMBR"] = 123
    @test replace_env!(totdict) === nothing

    # Check that substitution worked
    @test totdict["NAMGRIB"]["NENSFNB"] == 123

    # Check that call to dict2namelist succeeds
    @test dict2namelist(stdout,totdict) === nothing

end

# Check that the mechanism where we optionally take variables from the ENVironment works
@testset "ENV?default" begin
    dicts = read_namelist.(["arome_screening"])
    totdict = merge_namelists(dicts)
    @test replace_env!(totdict)=== nothing
    @test totdict["NAMMKODB"]["NTBMAR"]  == -90
    @test totdict["NAMMKODB"]["NTFMAR"]  == 90

    # After setting ENV
    Harmonie_namelists.MYENV["NTBMAR"] = 1
    Harmonie_namelists.MYENV["NTFMAR"] = 1
    dicts = read_namelist.(["arome_screening"])
    totdict = merge_namelists(dicts)
    @test replace_env!(totdict)=== nothing
    @test totdict["NAMMKODB"]["NTBMAR"]  == 1
    @test totdict["NAMMKODB"]["NTFMAR"]  == 1
end

BINDIR="../bin"
@testset "Namelists " begin
    # This should be improved we only test that gen_namelist and
    # Get_namelist produce output
    @test !isempty(read(`$BINDIR/gen_namelist.jl args`))
    @test !isempty(read(`$BINDIR/Get_namelist.jl`))
end


text:
## --------------------------------------------------------
function find_imports(pkgdir; subs = ["src", "scripts"])

    SPACE = hex_escape(" ")
    TAB = hex_escape("\t")
    ENTER = hex_escape("\n")
    COMMA = hex_escape(",")
    SEMICOLON = hex_escape(";")
    DOUBLEPOINT = hex_escape(":")
    BLANK = "(?:$SPACE|$TAB)"
    PKG_NAME = "[^$SPACE$TAB$COMMA$SEMICOLON$ENTER]+"

    imports_pool = Dict{String, Set{String}}()
    for dirname in subs
        srcdir = joinpath(pkgdir, dirname)
        !isdir(srcdir) && continue
        imports = get!(imports_pool, dirname, Set{String}())

        jl_files = filterdown(srcdir) do path
            endswith(path, ".jl")
        end

        for file in jl_files
            file_str = read(file, String)

            # @info("Doing", file)

            # predigest
            _repeat(file_str) do
                file_str = replace(file_str, r"(?:\"\"\")[^\"]*(?:\"\"\")"=>"")
                file_str = replace(file_str, r"(?:\")[^\"]*(?:\")"=>"")
                file_str = replace(file_str, r"\#\=(.|\n)*\=\#"=>"")
                file_str = replace(file_str, r"\#[^\n\"]*\n"=>"\n")
                file_str = replace(file_str, r"\#[^\n\"]*"=>"")
                file_str = replace(file_str, r"\n\n"=>"\n")
            end

            # import/using Pkg [\n;:]
            reg_str = string(
                "(?:import|using)", "(?:$BLANK|$ENTER)+", "(?<pkg>$PKG_NAME)", "$BLANK*", "[\n;:]"
            )
            reg = Regex(reg_str)

            _repeat(file_str) do
                m = match(reg, file_str)
                isnothing(m) && return file_str
                file_str = replace(file_str, m.match => "")

                # @show m
                # clear Package name
                pkg_str = m[:pkg]
                startswith(pkg_str, ".") && return file_str
                pkg_str = first(split(pkg_str, "."))
                push!(imports, pkg_str)
                return file_str
            end

            # import/using Pkg1 [\n] , [\n] Pkg1
            reg_str = string(
                "(?:import|using)", "(?:$BLANK|$ENTER)+", "(?:$PKG_NAME)", "$BLANK*",
                "(?:", "$BLANK*", "$COMMA", "(?:$BLANK|$ENTER)*",  "(?:$PKG_NAME)", ")*",
                "(?:$BLANK|$ENTER|$SEMICOLON)*"
            )
            reg = Regex(reg_str)
            _repeat(file_str) do
                m = match(reg, file_str)
                isnothing(m) && return file_str
                file_str = replace(file_str, m.match => "")

                # @show m
                # clear Package names
                imp_str = m.match
                clear_reg = Regex("import|using|$ENTER|$TAB|$SPACE|$SEMICOLON")
                _repeat(file_str) do
                    imp_str = replace(imp_str, clear_reg => "")
                end

                pkg_strs = split(imp_str, ",")
                for pkg_str in pkg_strs
                    startswith(pkg_str, ".") && continue
                    pkg_str = first(split(pkg_str, "."))
                    push!(imports, pkg_str)
                end
                return file_str
            end

        end # for file
    end # for dirname

    imports_pool
end


text:
module SumProductNetworks

# loading dependencies into workspaces
using Reexport
@reexport using Distributions
using StatsFuns
using SpecialFunctions
using AxisArrays
using DataStructures
using Distances
using Clustering
using StatsBase: countmap

import Base: getindex, map, parent, length, size, show, isequal, getindex, keys, eltype, rand
import Distributions.logpdf
import StatsBase.nobs

# add Base modules
@reexport using Statistics
using LinearAlgebra
using SparseArrays
using Random
using Printf

# include general implementations
include("nodes.jl")
include("nodeFunctions.jl")
include("networkFunctions.jl")

# include approach specific implementations
include("bmiTest.jl")
include("structureUtilities.jl")

# include utilities
include("utilityFunctions.jl")
include("io.jl")

end # module


text:
mutable struct DataStats
    num_observations::Int
    num_dimensions::Int
    train_fraction::Float64
    test_fraction::Float64
    train_indices::Vector{Int}
    test_indices::Vector{Int}
    DataStats() = new(0, 0, 0.0, 0.0, Int[], Int[])
    DataStats(num_observations::Int, num_dimensions::Int, train_fraction::Float64, test_fraction::Float64, train_indices::Vector{Any}, test_indices::Vector{Any}) = new(num_observations, num_dimensions, train_fraction, test_fraction, convert(Vector{Int}, train_indices), convert(Vector{Int}, test_indices))
    DataStats(num_observations::Int, num_dimensions::Int, train_fraction::Float64, test_fraction::Float64, train_indices::Vector{Int}, test_indices::Vector{Int}) = new(num_observations, num_dimensions, train_fraction, test_fraction, train_indices, test_indices)
end

struct Result
    id
    experiment::Dict
    worker_info::Dict{Symbol, String}
    data_stats::DataStats
    al_history::MVHistory
    al_summary::Dict
    status::Dict{Symbol, Symbol}
    Result(experiment) = new(experiment[:hash],
                     experiment,
                     Dict{Symbol, String}(),
                     DataStats(),
                     MVHistory(),
                     Dict{Symbol, Dict{Symbol, Any}}(),
                     Dict(:exit_code => :initialized))
    Result(id, experiment, worker_info, al_history, al_summary) = new(id, experiment, worker_info, DataStats(), al_history, al_summary, Dict(:exit_code => :initialized))
    Result(id, experiment, worker_info, data_stats, al_history, al_summary) = new(id, experiment, worker_info, data_stats, al_history, al_summary, Dict(:exit_code => :initialized))
    Result(id, experiment, worker_info, data_stats, al_history, al_summary, status) = new(id, experiment, worker_info, data_stats, al_history, al_summary, status)
end

function set_data_stats!(res::Result, x::Array{T, 2}, ds::DataSplits) where T <: Real
    res.data_stats.num_observations = size(x, 2)
    res.data_stats.num_dimensions = size(x, 1)
    res.data_stats.train_fraction = sum(ds.train) / res.data_stats.num_observations
    res.data_stats.test_fraction = sum(ds.test) / res.data_stats.num_observations
    res.data_stats.train_indices = findall(ds.train)
    res.data_stats.test_indices = findall(ds.test)
    return nothing
end

function set_model_fitted!(res, model)
    res.experiment[:model][:fitted] = Dict(:kernel => get_kernel(model),
                              :model_params => get_model_params(model))
    return nothing
end

function get_worker_info(;debug = false)
    worker_info = Dict{Symbol, String}()
    worker_info[:julia_version] = string(VERSION)
    debug && (worker_info[:installed_packages] = join( ["$k [$v]" for (k,v) in Pkg.installed()], ','))
    cmd = Sys.iswindows() ? `cmd /c git -C $(@__DIR__) rev-parse HEAD` : `git -C $(@__DIR__) rev-parse HEAD`
    worker_info[:git_commit] = strip(string(read(cmd, String)))
    worker_info[:utc_time] = string(Dates.now(Dates.UTC))
    worker_info[:hostname] = gethostname()
    return worker_info
end

function set_worker_info!(res::Result)
    merge!(res.worker_info, OneClassActiveLearning.get_worker_info())
    return nothing
end

function al_summarize!(res::Result)
    for e in keys(res.al_history)
        # Skip calculation of summary statistics for some metrics
        if e ∈ [:cm, :time_fit, :mem_fit, :time_qs, :mem_qs, :query_label,
            :query_history, :runtime, :time_set_data]
            continue
        end
        res.al_summary[e] = Dict()
        scores = values(res.al_history, e)
        score_changes = scores[2:end] - scores[1:end-1]
        res.al_summary[e][:start_quality] = scores[1]
        res.al_summary[e][:end_quality] = scores[end]
        res.al_summary[e][:maximum] = maximum(scores)
        res.al_summary[e][:ramp_up] = scores[2:end] .- scores[1]
        res.al_summary[e][:quality_range] = scores[end] .- scores[1]
        res.al_summary[e][:total_quality_range] = maximum(scores) - minimum(scores)
        res.al_summary[e][:average_end_quality] = [mean(scores[k:end]) for k in 1:length(scores)]
        res.al_summary[e][:average_quality_change] = mean(score_changes)
        if res.al_summary[e][:quality_range] > 0
            res.al_summary[e][:learning_stability] = [((scores[end] - scores[end-k]) / k) / (res.al_summary[e][:quality_range] / (length(scores) - 1)) for k in 1:(length(scores) - 1)]
        else
            res.al_summary[e][:learning_stability] = [0.0 for k in 1:(length(scores) - 1)]
        end
        res.al_summary[e][:average_gain] = any(score_changes .> 0) ? mean(score_changes[score_changes .> 0]) : 0
        res.al_summary[e][:average_loss] = any(score_changes .< 0) ? mean(score_changes[score_changes .< 0]) : 0
        query_labels = values(res.al_history, :query_label)
        if :outlier in query_labels
            res.al_summary[e][:ratio_of_outlier_queries] = sum(query_labels .== :outlier) / length(query_labels)
        else
            res.al_summary[e][:ratio_of_outlier_queries] = 0.0
        end
    end
    return nothing
end


text:
module TestExpressionsSafeTupleName

import Test: @testset, @test
import Volcanito: safe_tuple_name

@testset "safe_tuple_name" begin
    @test safe_tuple_name(:(x + y)) == :t
    @test safe_tuple_name(:(x + t)) != :t
end

end


text:
using Pkg
Pkg.activate(@__DIR__)
pkg"up"
pkg"instantiate"

using LibPQ, JSON3

using LibPQ: load!

using HTTP: request

ENV["SAM_KEY"] = "MwvGHSGiyuISuxArFq5wMQ2K2W4G1FhcMElsvCJ2"

conn = LibPQ.Connection("""
                        host = $(get(ENV, "PGHOST", ""))
                        dbname = sdad
                        user = $(get(ENV, "DB_USR", ""))
                        password = $(get(ENV, "DB_PWD", ""))
                        """);

execute(conn,
        """
        CREATE SCHEMA IF NOT EXISTS us_sam
        AUTHORIZATION ncses_oss;
        COMMENT ON SCHEMA us_sam
        IS 'Beta.SAM.Gov Federal Hierarchy Public API
        https://open.gsa.gov/api/fh-public-api/';
        DROP TABLE IF EXISTS us_sam.fh_public_api CASCADE;
        CREATE TABLE IF NOT EXISTS us_sam.fh_public_api (
            fhorgid text,
            fhorgtype text NOT NULL,
            fhorgname text NOT NULL,
            fhdeptindagencyorgid text NOT NULL,
            PRIMARY KEY (fhorgid)
        );
        ALTER TABLE us_sam.fh_public_api OWNER to ncses_oss;
        """
    )
f
helper(elem) =
        (fhorgid = elem.fhorgid,
         fhorgname = :fhorgname ∈ propertynames(elem) ? elem.fhorgname : missing,
         lastupdateddate = :lastupdateddate ∈ propertynames(elem) ? elem.lastupdateddate : missing)
for file in readdir(joinpath(@__DIR__, "data", "oss", "original", "us_sam_fh_public_api"))
    # file = readdir(joinpath(@__DIR__, "data", "oss", "original", "us_sam_fh_public_api"))[1]
    # file = joinpath(@__DIR__, "data", "oss", "original", "us_sam_fh_public_api", "100003657.json")
    println(file)
    obj = JSON3.read(open(joinpath(@__DIR__, "data", "oss", "original", "us_sam_fh_public_api", file)))
    fhdeptindagencyorgid = obj[1].fhdeptindagencyorgid
    fhorgtype = "Sub-Tier"
    data = DataFrame(helper(elem) for elem in obj)
    data = data[.!ismissing.(data.fhorgname),:]
    sort!(data, (:fhorgname, order(:lastupdateddate, rev=true)))
    data = by(data[!,[:fhorgid, :fhorgname]], :fhorgname, x -> first(x, 1))
    data[!,:fhdeptindagencyorgid] .= obj[1].fhdeptindagencyorgid
    data[!,:fhorgtype] .= "Sub-Tier"
    push!(data, (fhorgid = obj[1].fhdeptindagencyorgid,
                 fhorgname = obj[1].fhorgname,
                 fhdeptindagencyorgid = obj[1].fhdeptindagencyorgid,
                 fhorgtype = obj[1].fhorgtype))
    data = sort(data[!,[:fhorgid, :fhorgtype, :fhorgname, :fhdeptindagencyorgid]])
    load!(data, conn, "INSERT INTO us_sam.fh_public_api VALUES ($(join(("\$$i" for i in 1:4), ','))) ON CONFLICT DO NOTHING;")
end

#
CREATE MATERIALIZED VIEW IF NOT EXISTS us_sam.fh_public_clean AS (
	WITH A AS (
		SELECT *
		FROM us_sam.fh_public_api
		WHERE fhdeptindagencyorgid <> '100000000'
	),
	B AS (
		SELECT DISTINCT fhdeptindagencyorgid, fhorgname
		FROM A
		WHERE fhorgtype = 'Department/Ind. Agency'
	)
	SELECT B.fhorgname AS parent, A.fhorgtype, A.fhorgname
	FROM B
	INNER JOIN A
	ON B.fhdeptindagencyorgid = A.fhdeptindagencyorgid
	ORDER BY parent ASC, fhorgtype ASC, fhorgname ASC
)
;
#

    for elem in obj
        lastupdateddate = :lastupdateddate ∈ propertynames(elem) ? elem.lastupdateddate : missing
        push!(output,
              (string(elem.fhorgid),
               elem.fhorgname,
               elem.fhorgtype,
               string(elem.fhdeptindagencyorgid),
               lastupdateddate
               )
              )
    end
    sort!(output, by = x -> (x[2], x[end]))
    for fhorgname in (elem.fhorgname for elem in output)
        idx = findall(isequal(fhorgname), [])

    tmp = [ (elem.fhorgid, elem.fhorgname, elem.fhorgtype, elem.fhdeptindagencyorgid, elem.lastupdateddate)
            for elem in obj ]
    findall(elem -> elem[2] == "VIRGINIA CONTRACTING AGENCY", tmp)
    obj[596]
    obj[597]

end
d
"fhorgid":100000000,
"fhorgname":"DEPT OF DEFENSE",
"fhorgtype":"Department/Ind. Agency",
"status":"ACTIVE",
"createddate":"2006-04-14 00:00",
"updatedby":"FPDS_START_DATE_FIX",
"lastupdateddate":"2019-11-21 23:04",
"fhdeptindagencyorgid":100000000,
"fhagencyorgname":"DEPT OF DEFENSE",
"agencycode":"9700",
"oldfpdsofficecode":"9700",
"cgaclist":[{"cgac":"057"},{"cgac":"021"},{"cgac":"017"},{"cgac":"096"},{"cgac":"097"}],
"fhorgnamehistory":[{"fhorgname":"DEPT OF DEFENSE","effectivedate":"1957-07-01 00:00"}],
"fhorgparenthistory":[{"fhfullparentpathid":"100000000","fhfullparentpathname":"DEPT OF DEFENSE","effectivedate":"1957-07-01 00:00"}],
"links":[{"rel":"self","href":"https://api.sam.gov/prod/federalorganizations/v1/orgs?fhorgid=100000000","hreflang":null,"media":null,"title":null,"type":null,"deprecation":null},
         {"rel":"nextlevelchildren","href":"https://api.sam.gov/prod/federalorganizations/v1/org/hierarchy?fhorgid=100000000","hreflang":null,"media":null,"title":null,"type":null,"deprecation":null}
        ]
}

response = request("GET",
                   "https://api.sam.gov/prod/federalorganizations/v1/orgs?api_key=$(ENV["SAM_KEY"])&status=active&limit=200")
response.status = 200
json = JSON3.read(response.body)
length(json.orglist)

function find_all_records(output, url)
    limit = parse(Int, match(r"(?<=limit=)\d+", url).match)
    current_offset = parse(Int, match(r"(?<=offset=)\d+", url).match)
    url = replace(url, r"(?<=offset=)\d+" => current_offset + limit)
    response = request("GET", url)
    # output, url = vcat(output, 1), url
    json = JSON3.read(response.body)
    if length(json.orglist) < 200
        println(response)
    end
    println("A: $(length(output))")
    output = vcat(output, json.orglist)
    println("B: $(length(output))")
    output, url
end
function find_all_records(key = ENV["SAM_KEY"])
    url = "https://api.sam.gov/prod/federalorganizations/v1/orgs?api_key=$key&status=active&limit=200"
    response = request("GET", url)
    json = JSON3.read(response.body)
    totalrecords = json.totalrecords
    output = json.orglist
    url = string(url, "&offset=200")

    # response_200 = request("GET", url)
    # json_200 = JSON3.read(response_200.body)
    # totalrecords_200 = json_200.totalrecords
    # output_200 = json_200.orglist
    # url = replace(url, r"(?<=&offset=).*" => 400)

    # response_400 = request("GET", url)
    # json_400 = JSON3.read(response_400.body)
    # totalrecords_400 = json_400.totalrecords
    # output_400 = json_400.orglist
    # url = replace(url, r"(?<=&offset=).*" => 600)
    url == "https://api.sam.gov/prod/federalorganizations/v1/orgs?api_key=MwvGHSGiyuISuxArFq5wMQ2K2W4G1FhcMElsvCJ2&status=active&limit=200&offset=600"
    # response_600 = request("GET", url)
    # json_600 = JSON3.read(response_600.body)
    # totalrecords_600 = json_600.totalrecords
    # output_600 = json_600.orglist
    # url = replace(url, r"(?<=&offset=).*" => 800)

    # response_800 = request("GET", url)
    # json_800 = JSON3.read(response_800.body)
    # totalrecords_800 = json_800.totalrecords
    # output_800 = json_800.orglist
    # url = replace(url, r"(?<=&offset=).*" => 1000)

    # awesome = vcat(output, output_200, output_400, output_600, output_800)
    # length(awesome)

    # totalrecords = 3
    # output = [0]
    sleep(7)
    while length(output) < totalrecords
        println(url)
        output, url = find_all_records(output, url)
        println(length(output))
        sleep(7)
    end
    output
end
orgs = find_all_records()
io = open(joinpath(@__DIR__, "data", "oss", "original", "gsa_federalorganizations", "orglist.json"), write = true)
propertynames(io)
JSON3.write(io, awesome)

orgs = JSON3.read(read(joinpath(@__DIR__, "data", "oss", "original", "gsa_federalorganizations", "orglist.json"), String))
org = orgs[1]
length(orgs)
org
setdiff(orgs, orgs_parents)[1]
findfirst(elem -> elem.fhorgid == 500154577, orgs)
findfirst(elem -> elem.fhorgid == 500155065, orgs)

already = replace.(filter!(filename -> occursin(r"\d+(?=\.json$)", filename), readdir(joinpath(@__DIR__, "data", "oss", "original", "gsa_federalorganizations"))),
                   ".json" => "")
orgs_parents = filter(elem -> length(elem.links) == 2, orgs)
length(orgs_parents)
second_tier = sort!([elem.links[2].href[72:end] for elem in orgs_parents])
to_process = setdiff(second_tier, already)
# This one is weird
filter!(!isequal("300000201"), to_process)
intersect(second_tier, already)[1]
arg = filter(filename -> occursin(r"\d+(?=.json)", filename), readdir(joinpath(@__DIR__, "data", "oss", "original", "gsa_federalorganizations")))

maximum(x -> findfirst(elem -> elem.fhorgid == parse(Int, replace(x, ".json" => "")), orgs), arg)
length("https://api.sam.gov/prod/federalorganizations/v1/org/hierarchy?fhorgid=")
org = to_process[1]
for org in orgs[304:end]
    idx = findfirst(elem -> isequal("nextlevelchildren", elem.rel), org.links)
    if isa(idx, Integer)
        url = string(org.links[idx].href, "&api_key=$key&limit=200")
        fhorgid = org.fhorgid
        response = request("GET", url)
        json = JSON3.read(response.body)
        totalrecords = json.totalrecords
        output = json.orglist
        url = string(url, "&offset=200")
        while length(output) < totalrecords
            response = request("GET", url)
            json = JSON3.read(response.body)
            output = vcat(output, json.orglist)
            limit = parse(Int, match(r"(?<=limit=)\d+", url).match)
            current_offset = parse(Int, match(r"(?<=offset=)\d+", url).match)
            url = replace(url, r"(?<=offset=)\d+" => current_offset + limit)
        end
        io = open(joinpath(@__DIR__, "data", "oss", "original", "gsa_federalorganizations", "$(org.fhorgid).json"), write = true)
        JSON3.write(io, output)
        close(io)
        println(org.fhorgid)
    end
end
org = "300000201"
org.fhorgid
org.links[2].href

for org in to_process
    url = "https://api.sam.gov/prod/federalorganizations/v1/org/hierarchy?fhorgid=$org&status=Active&api_key=$key&limit=200"
    fhorgid = org
    response = request("GET", url)
    json = JSON3.read(response.body)
    totalrecords = json.totalrecords
    println(totalrecords)
    output = json.orglist[1]
    url = string(url, "&offset=200")
    while length(output) < totalrecords
        response = request("GET", url)
        json = JSON3.read(response.body)
        output = vcat(output, json.orglist)
        limit = parse(Int, match(r"(?<=limit=)\d+", url).match)
        current_offset = parse(Int, match(r"(?<=offset=)\d+", url).match)
        url = replace(url, r"(?<=offset=)\d+" => current_offset + limit)
    end
    io = open(joinpath(@__DIR__, "data", "oss", "original", "gsa_federalorganizations", "$org.json"), write = true)
    JSON3.write(io, output)
    close(io)
    println(org)
end

    if length(orgs[1].links) == 2
        @assert orgs[1].links[2].rel == "nextlevelchildren"

json.orglist
response_hier = request("GET",
                        "https://api.sam.gov/prod/federalorganizations/v1/org/hierarchy?api_key=$(ENV["SAM_KEY"])&fhorgid=300000015&status=active&limit=100")
json_hier = JSON3.read(response_hier.body)

response_chk = request("GET",
                        "https://api.sam.gov/prod/federalorganizations/v1/org/hierarchy?api_key=$(ENV["SAM_KEY"])&fhorgid=300000117&status=active&limit=100")
json_chk = JSON3.read(response_chk.body)
json_hier.orglist
json_chk.orglist |> length
close(conn)


text:
using Test
using DefaultArrays
using LinearAlgebra
@testset "Basic Tests" begin
    A=DefaultArray(Inf,(2,3))
    A[1,2]=2
    @test maximum(A)==Inf
    @test minimum(A)==2
    @test A[3]==A[1,2]
    A[4]=3.0
    @test A[4]==A[2,2]
    @test all(collect(eachnondefault(A)).==[4, 3])
    @test collect(eachnondefault(A,IndexCartesian())) == [CartesianIndex(2, 2), CartesianIndex(1, 2)]
    l1=length(eachnondefault(A)|>collect)
    A[1,2]=Inf
    l2=length(eachnondefault(A)|>collect)
    @test l1>l2
    B=DefaultArray(Inf,(3,))
    B[3]=2
    @test B[3]==2
end

@testset "similar(..) Tests + broadcasting" begin
    M=similar(DefaultArray(1.0),Float64,5,5)

    @test typeof(M)==DefaultArray{Float64,2}
    @test all(@inferred  M.==1.0)
    M[1,2]=0.5
    @test typeof(@inferred M.+2)==DefaultArray{Float64,2}
    M.=M.+M'

    @test all(M.==[2.0 1.5 2.0 2.0 2.0; 1.5 2.0 2.0 2.0 2.0; 2.0 2.0 2.0 2.0 2.0; 2.0 2.0 2.0 2.0 2.0; 2.0 2.0 2.0 2.0 2.0])
    M[1:2,1:2].=0
    @test all(M.==[0.0 0.0 2.0 2.0 2.0; 0.0 0.0 2.0 2.0 2.0; 2.0 2.0 2.0 2.0 2.0; 2.0 2.0 2.0 2.0 2.0; 2.0 2.0 2.0 2.0 2.0])
    M[4:5,4:5].=0
    @test all(M.==[0.0 0.0 2.0 2.0 2.0; 0.0 0.0 2.0 2.0 2.0; 2.0 2.0 2.0 2.0 2.0; 2.0 2.0 2.0 0.0 0.0; 2.0 2.0 2.0 0.0 0.0])
    Q=M[2:3,2:3].+I(2)
    @test all(Q.==[1.0 2.0; 2.0 3.0])
    @test typeof(M)==DefaultArray{Float64,2}
    @test all((DefaultArray(0.0,2,2).+true).==1)
end

@testset "Various Constructors" begin
    A=[1 1;0 1]
    M=DefaultArray(A)
    @test all(A.==M)
    @test M==DefaultArray(M)
    ugly=DefaultArray(0,(2,2)) .+ Any[1 1;2 2]
    @test repr((typeof(ugly),ugly))=="(DefaultArray{Int64,2}, [1 1; 2 2])"
end

max_stars_repo_name: francescoalemanno/DefaultArrays

max_issues_repo_name: francescoalemanno/DefaultArrays

max_forks_repo_name: francescoalemanno/DefaultArrays


text:
module CompressedCellValuesTests

using Test
using Gridap
using Gridap.CellValuesGallery
using ..MapsMocks

values = [10,20,31]
ptrs = [1,2,3,3,2,2]
r = values[ptrs]
cv = IterCompressedCellValue(values,ptrs)
test_iter_cell_value(cv,r)

cv = IndexCompressedCellValue(values,ptrs)
test_index_cell_value(cv,r)

l = 10
v = 3
ccv = ConstantCellValue(v,l)
cv = IndexCompressedCellValue(ccv)
test_index_cell_value(cv,fill(v,l))

@test cv == cv
@test cv ≈ cv

values = [10,20,31]
ptrs = [1,2,3,3,2,2]
cv = CompressedCellValue(values,ptrs)

cv2 = apply(-,cv)
@test isa(cv2,IndexCompressedCellValue)
r = -values[ptrs]
test_index_cell_value(cv2,r)

values1 = [10,20,31]
values2 = [11,21,51]
ptrs = [1,2,3,3,2,2]
cv1 = CompressedCellValue(values1,ptrs)
cv2 = CompressedCellValue(values2,ptrs)
cv3 = apply(-,cv1,cv2)
@test isa(cv3,IndexCompressedCellValue)
r = values1[ptrs] - values2[ptrs]
test_index_cell_value(cv3,r)

values1 = [10,20,31]
values2 = [11,21,51]
ptrs1 = [1,2,3,3,2,2]
ptrs2 = [1,1,3,2,3,2]
cv1 = CompressedCellValue(values1,ptrs1)
cv2 = CompressedCellValue(values2,ptrs2)
cv3 = apply(-,cv1,cv2)
@test ! isa(cv3,IndexCompressedCellValue)
r = values1[ptrs1] - values2[ptrs2]
test_index_cell_value(cv3,r)

values1 = [10,20,31]
values2 = [11,21,51]
ptrs1 = [1,2,3,3,2,2]
ptrs2 = [1,1,3,2,3,2]
cv1 = IterCompressedCellValue(values1,ptrs1)
cv2 = IterCompressedCellValue(values2,ptrs2)
cv3 = apply(-,cv1,cv2)
@test ! isa(cv3,IndexCompressedCellValue)
r = values1[ptrs1] - values2[ptrs2]
test_iter_cell_value(cv3,r)

values = [[10,30],[20,10,40],[31,]]
ptrs = [1,2,3,3,2,2]
cv = CompressedCellValue(values,ptrs)

cv2 = apply(-,cv,broadcast=true)
@test isa(cv2,IndexCompressedCellValue)
r = -values[ptrs]
test_index_cell_array(cv2,r)

values1 = [[10,33],[20,12,40],[71,]]
values2 = [[10,30],[20,10,40],[31,]]
ptrs = [1,2,3,3,2,2]
cv1 = CompressedCellValue(values1,ptrs)
cv2 = CompressedCellValue(values2,ptrs)
cv3 = apply(-,cv1,cv2,broadcast=true)
@test isa(cv3,IndexCompressedCellValue)
r = values1[ptrs] - values2[ptrs]
test_index_cell_array(cv3,r)

values1 = [[10,33,2],[20,12,40],[71,1,4]]
values2 = [[10,30,3],[20,10,40],[31,3,5]]
ptrs1 = [1,3,2,3,2,2]
ptrs2 = [1,2,3,3,2,2]
cv1 = CompressedCellValue(values1,ptrs1)
cv2 = CompressedCellValue(values2,ptrs2)
cv3 = apply(-,cv1,cv2,broadcast=true)
r = values1[ptrs1] - values2[ptrs2]
@test ! isa(cv3,IndexCompressedCellValue)
test_index_cell_array(cv3,r)

a = VectorValue(10,10)
b = VectorValue(15,20)
p1 = VectorValue(1,1)
p2 = VectorValue(2,2)
p3 = VectorValue(3,3)
p = [p1,p2,p3]
m = MockMap(a)

l = 5
cm = ConstantCellMap(m,l)
vals = [[p1,p3], [p1,p3,p2], ]
ptrs = [1,2,1,1,2]
ca = CompressedCellValue(vals,ptrs)
cv = evaluate(cm,ca)
@test isa(cv,CompressedCellValue)
rvals = [ evaluate(m,pj) for pj in vals  ]
cr = rvals[ptrs]
test_index_cell_map_with_index_arg(cm,ca,cr)

m1 = MockMap(a)
m2 = MockMap(b)
mvals = [m1,m2]
pvals = [[p1,p3], [p1,p3,p2], ]
ptrs = [1,2,1,1,2]
cm = CompressedCellMap(mvals,ptrs)
ca = CompressedCellArray(pvals,ptrs)
cv = evaluate(cm,ca)
@test isa(cv,CompressedCellValue)
rvals = [ evaluate(mi,ai) for (mi,ai) in zip(mvals,pvals) ]
cr = CompressedCellValue(rvals,ptrs)
test_index_cell_map_with_index_arg(cm,ca,cr)

values = [10,20,31]
ptrs = [1,2,3,3,2,2]
inds = [2,1,4]
cv = CompressedCellValue(values,ptrs)
indices = CellValueFromArray(inds)
cv2 = reindex(cv,indices)
@test isa(cv2,CompressedCellValue)
r = values[ptrs[inds]]
test_index_cell_value(cv2,r)

f = ( x -> 3*x )
fca = broadcast(f,cv)
@test all(fca .== cv*3)

end # module


text:
using Gadfly
using Statistics
using Clustering
using DataFrames

import Distances: pairwise, Euclidean
import StatsBase: cov2cor
import RLEVectors: rep

#include("helpers.jl")

export plot_corr
function plot_corr(chain, m::Int64, labs::Array{String,1}; reorder = true, key = false, min_value = -1, max_value = 1, linkage = :ward, reds = false)
"""
    Plot various correlation matrices
    1. Obtain some covariance chain and average over the estimates
    2. Convert to a correlation matrix
    3. Cluster using hierarchical clustering, and reorder the dimensions based on the output
    4. Convert to a molten data frame and plot
"""
    Schain = get_Sigma_chain(chain, m)
    covmat = mapslices(x -> mean(x), Schain; dims = 3)[:,:,1]
    cormat = cov2cor(covmat, sqrt.(diag(covmat)))

    if reorder
            D = pairwise(Euclidean(), cormat)
            try
                    cl = Clustering.hclust(D, linkage=linkage)
                    ordering = getproperty(cl, :order)
            catch
                    ordering = 1:size(cormat,1)
            end
    else
            ordering = 1:size(cormat,1)
    end

    df = DataFrame(cormat[ordering, ordering], Symbol.(labs[ordering]))
    molten = melt(df)
    molten.v2 = rep(labs[ordering], times = size(df,1))
    if reds
            Gadfly.plot(molten, x="variable", y="v2", color="value",
                Geom.rectbin,
                Coord.cartesian(fixed = true),
                Guide.xlabel(nothing), Guide.ylabel(nothing),
                Scale.color_continuous(minvalue=min_value, maxvalue=max_value, colormap = Scale.lab_gradient("white", "red")),
                Guide.colorkey(title=""),
                key ? Theme(key_position = :right) : Theme(key_position = :none))
        else
                Gadfly.plot(molten, x="variable", y="v2", color="value",
                    Geom.rectbin,
                    Coord.cartesian(fixed = true),
                    Guide.xlabel(nothing), Guide.ylabel(nothing),
                    Scale.color_continuous(minvalue=min_value, maxvalue=max_value),
                    Guide.colorkey(title=""),
                    key ? Theme(key_position = :right) : Theme(key_position = :none))
        end
end

export plot_eigvec
function plot_eigvec(chain, m::Int64, labs::Array{String,1}, eig_num::Int64; pal=missing, reorder=missing)
        Schain = get_Sigma_chain(chain, m)
        covmat = mapslices(x -> mean(x), Schain; dims = 3)[:,:,1]

        if !ismissing(reorder)
                covmat = covmat[reorder,reorder]
                labs = labs[reorder]
                pal = pal[reorder]
        end

        evals, evecs = (eigvals(covmat), eigvecs(covmat))

        # Make names dataframe so plot has names
        # The kth ROW is the kth eigenvector here
        edf = DataFrame(evecs', Symbol.(labs))
        edf.vec = 1:size(covmat, 1)
        molten = melt(edf)

        # Plot first eigenvector
        if ismissing(pal)
                p = Gadfly.plot(molten[molten.vec .== eig_num, :],
                        x=:variable,
                        y=:value,
                        Guide.xlabel(nothing), Guide.ylabel(nothing),
                        Geom.bar)#,
                        #yintercept = 0,
                        #Geom.hline(color="black", style = :dash))
        else
                p = Gadfly.plot(molten[molten.vec .== eig_num, :],
                        x=:variable,
                        y=:value,
                        color=pal,
                        Guide.xlabel(nothing), Guide.ylabel(nothing),
                        Geom.bar)#,
                        #yintercept = 0,
                        #Geom.hline(color="black", style = :dash))
        end
        (p, sort(evals; rev=true))
end

export plot_effects
function plot_effects(chain, obsidx::Int64, labs::Array{String,1}, M::Int64, mu_ests; pal = missing)
        z_chain = get_z_chain(chain)
        MAPs = mapslices(x -> [ mean(x .== y) for y in 1:M ], z_chain; dims = 2)
        emap = MAPs[obsidx,:]

        # Effect size estimates in each cell type,
        #       weighted by the posterior probability the this site is in any given cluster
        #       plus the standard errors
        emeans = mapslices(x -> sum(x), hcat([ mu_ests[m,:] .* emap[m] for m in 1:M ]...)'; dims = 1)[1,:]
        eses = mapslices(y -> sum(y), hcat([ mapslices(x -> std(x), get_mu_chain(chain, m);
                                dims = 2) .* emap[m] for m in 1:M ]...)'; dims = 1)[1,:]

        # Plot first eigenvector
        if ismissing(pal)
                Gadfly.plot(y = labs, x = emeans,
                                xmin = emeans .- eses,
                                xmax = emeans .+ eses,
                                Geom.point,
                                Guide.xlabel(nothing), Guide.ylabel(nothing),
                                Geom.xerrorbar)
        else
                Gadfly.plot(y = labs, x = emeans,
                                xmin = emeans .- eses,
                                xmax = emeans .+ eses,
                                Geom.point,
                                Geom.xerrorbar,
                                Guide.xlabel(nothing), Guide.ylabel(nothing),
                                color = pal)
        end
end


text:
export  DeleteFrom,
        setreturning!,
        deletefrom

"""
```
DeleteFrom( name::Union{Var, AbstractString}, cond::SQLCondition=nocond )
```
This mutable struct represents a `DELETE FROM` SQL statement where `name` is the name of the affected table, and `cond` is the SQL condition that describes which entries should be deleted from the table.
"""
mutable struct DeleteFrom <: SQLStatement
    name::String
    wherec::SQLCondition
    returning::Vector{String}

    DeleteFrom( name::SVar, cond::SQLCondition=nocond ) = new( name, cond, String[] )
end

function Base.show( io::IO, df::DeleteFrom )
    print( io, "DELETE FROM \"", df.name, "\"" )
    isempty(df.wherec) || print( io, "\n    WHERE ", df.wherec )
    isempty(df.returning) || print( io, "\n    RETURNING ", join( df.returning, ", " ) )
end

"""
```
setreturning!( names::Vector{T} ) where T <: AbstractString
setreturning!( names::Vector{Var} )
setreturning!( names::Union{Var, AbstractString}... )
```
This function permits the user to add a command that sets the `RETURNING` clause of a `InsertInto`, `Update`, or `DeleteFrom` struct to `names`.

The output of this function is a function that takes an argument of the aforementioned types and returns the updated argument, such that Julia's function chaining syntax can be used.

Example:
```
sq |> setreturning!( "var1", "var2 + var3 AS varcomb" )
```
sets the commands to return `var1` and `var2 + var3`, aliased by `varcomb`.
"""
setreturning!( names::Vector{T} ) where T <: AbstractString =
    dm::Union{InsertInto, Update, DeleteFrom} -> begin
        dm.returning = names |> Vector{String}
        dm
    end
setreturning!( names::Vector{Var} ) = string.(names) |> setreturning!
setreturning!( names::SVar... ) = string.( names |> collect ) |> setreturning!

"""
```
deletefrom( conn::Conn, df::Bool=true )
```
This function permits the user to execute the SQL equivalent of a `DeleteFrom` struct on the database that `conn` connects to. The flag `df` determines the format of the output (`DataFrame` or `NamedTuple`) if the statement will return an output.

The output of this function is a function that takes a `DeleteFrom` as argument and returns the result of the executed SQL statement.

Example:
```
sq |> deletefrom(conn)
```
"""
deletefrom( conn::Conn, df::Bool=true ) = def::DeleteFrom -> begin
    res = LibPQ.execute( conn, string( def, ";" ) )
    isempty(def.returning) && return res
    res |> makeunique! |> tableresult(df)
end


text:
"""
Copied from `cleanup_yacn` branch in neuroglancer:
https://github.com/seung-lab/neuroglancer/blob/cleanup_yacn/python/ext/third_party/yacn/chunked_regiongraph/tasks.jl
"""

module Precomputed

using PyCall
@pyimport neuroglancer.pipeline as pl
const pyslice=pybuiltin(:slice)

function cached(f)
	cache=Dict()
	function my_f(args...)
		if !haskey(cache, args)
			cache[args] = f(args...)
		else
			println("restoring from cache")
		end
		return cache[args]
	end
end

CachedStorage = cached(pl.Storage)

immutable PrecomputedWrapper
	val
	function PrecomputedWrapper(storage_string, scale_idx=0)
		return new(pl.Precomputed(CachedStorage(storage_string), scale_idx))
	end
end

function Base.getindex(x::PrecomputedWrapper, slicex::UnitRange,
                                        slicey::UnitRange, slicez::UnitRange)
    return squeeze(get(x.val,
            (pyslice(slicex.start,slicex.stop+1),
            pyslice(slicey.start,slicey.stop+1),
            pyslice(slicez.start,slicez.stop+1))),4)
end

end # module Precomputed


text:
# Create a special type for permutations. The real point here is to be able to
# unambiguously identify an RRPermArray (see below) so that we may "unwrap" in
# expressions like `channelview(colorview(C, A))`.

"""
ColorChanPerm(perm)

Construct a reordering permutation for the color channel.
This handles swaps between memory layout and constructor argument order for `AbstractRGB` and
various `AlphaChannel` and `ChannelAlpha` color types.
"""
struct ColorChanPerm{N} <: AbstractVector{Int}
    perm::NTuple{N,Int}
end
Base.IndexStyle(::Type{<:ColorChanPerm}) = IndexLinear()
Base.size(v::ColorChanPerm{N}) where N = (N,)
Base.getindex(v::ColorChanPerm, i::Int) = v.perm[i]

dimorder(::Type{<:RGB}) = ColorChanPerm((1, 2, 3))
dimorder(::Type{<:BGR}) = ColorChanPerm((3, 2, 1))
dimorder(::Type{<:XRGB}) = ColorChanPerm((1, 1, 2, 3))
dimorder(::Type{<:RGBX}) = ColorChanPerm((1, 2, 3, 3)) # this causes problems for setindex!, fixed below
dimorder(::Type{<:BGRA}) = ColorChanPerm((3, 2, 1, 4))
dimorder(::Type{<:ABGR}) = ColorChanPerm((4, 3, 2, 1))
dimorder(::Type{<:AlphaColor{<:Color1,T,N}}) where {T,N} = ColorChanPerm((2, 1))
dimorder(::Type{<:AlphaColor{<:Color3,T,N}}) where {T,N} = ColorChanPerm((4, 1, 2, 3))

const ColorChanPermIndexType{NC} = Tuple{<:ColorChanPerm,Vararg{<:Base.Slice,NC}}
const ColorChanPermSubArray{T,N,P,I<:ColorChanPermIndexType,L} =
    SubArray{T,N,P,I,L}
const RRPermArray{To,From,N,M,P<:ColorChanPermSubArray} =
    RRArray{To,From,N,M,P}

# This type exists solely to set multiple values in the color channel axis
struct NVector{T,N} <: AbstractVector{T}
    v::NTuple{N,T}
end
Base.IndexStyle(::Type{<:NVector}) = IndexLinear()
Base.size(v::NVector{T,N}) where {T,N} = (N,)
Base.getindex(v::NVector, i::Int) = v.v[i]
NVector(x::Vararg{T,N}) where {T,N} = NVector{T,N}(x)

@inline Base.setindex!(A::RRPermArray{<:RGBX,<:Number,N}, val::AbstractRGB, i::Vararg{Int,N}) where N =
    setindex!(parent(parent(parent(A))), NVector(red(val), green(val), blue(val)), :, i...)

"""
    channelview(A)

returns a view of `A`, splitting out (if necessary) the color channels
of `A` into a new first dimension.

Of relevance for types like RGB and BGR, the channels of the returned
array will be in constructor-argument order, not memory order (see
`reinterpretc` if you want to use memory order).

# Example
```julia
img = rand(RGB{N0f8}, 10, 10)
A = channelview(img)   # a 3×10×10 array

See also: [`colorview`](@ref)
"""
channelview(A::AbstractArray{T}) where {T<:Number} = A
channelview(A::RRPermArray{<:Colorant,<:Number}) = parent(parent(parent(A)))
channelview(A::RRArray{<:Colorant,<:Number}) = parent(parent(A))
channelview(A::Base.ReinterpretArray{<:AbstractGray,M,<:Number}) where M = parent(A)
channelview(A::AbstractArray{RGB{T}}) where {T} = reinterpretc(T, A)
function channelview(A::AbstractArray{C}) where {C<:AbstractRGB}
    # BGR, XRGB, etc don't satisfy conditions for reinterpret
    CRGB = RGB{eltype(C)}
    channelview(of_eltype(CRGB, A))
end
channelview(A::AbstractArray{C}) where {C<:Color} = reinterpretc(eltype(C), A)
channelview(A::AbstractArray{C}) where {C<:ColorAlpha} = _channelview(color_type(C), A)
_channelview(::Type{<:RGB}, A) = reinterpretc(eltype(eltype(A)), A)
function _channelview(::Type{C}, A) where {C<:AbstractRGB}
    CRGBA = RGBA{eltype(C)}
    channelview(of_eltype(CRGBA, A))
end
_channelview(::Type{C}, A) where {C<:Color} = reinterpretc(eltype(C), A)
function channelview(A::AbstractArray{AC}) where {AC<:AlphaColor}
    CA = coloralpha(base_color_type(AC)){eltype(AC)}
    channelview(of_eltype(CA, A))
end

"""
    colorview(C, A)

returns a view of the numeric array `A`, interpreting successive
elements of `A` as if they were channels of Colorant `C`.

Of relevance for types like RGB and BGR, the elements of `A` are
interpreted in constructor-argument order, not memory order (see
`reinterpretc` if you want to use memory order).

# Example
```jl
A = rand(3, 10, 10)
img = colorview(RGB, A)
```

See also: [`channelview`](@ref)
"""
colorview(::Type{C}, A::AbstractArray{T}) where {C<:Colorant,T<:Number} =
    _ccolorview(ccolor_number(C, T), A)
_ccolorview(::Type{C}, A::RRPermArray{T,C}) where {C<:Colorant,T<:Number} =
    parent(parent(parent(A)))
_ccolorview(::Type{C}, A::RRArray{T,C}) where {C<:Colorant,T<:Number} =
    parent(parent(A))
_ccolorview(::Type{C}, A::Base.ReinterpretArray{T,M,C}) where {C<:AbstractGray,T<:Number,M} =
    parent(A)
_ccolorview(::Type{C}, A::Base.ReinterpretArray{T,M,C}) where {C<:RGB,T<:Number,M} =
    reshape(parent(A), Base.tail(axes(parent(A))))
_ccolorview(::Type{C}, A::Base.ReinterpretArray{T,M,C}) where {C<:AbstractRGB,T<:Number,M} =
    _colorview_reorder(C, A)
_ccolorview(::Type{C}, A::Base.ReinterpretArray{T,M,C}) where {C<:Color,T<:Number,M} =
    reshape(parent(A), Base.tail(axes(parent(A))))
_ccolorview(::Type{C}, A::AbstractArray{T}) where {C<:Colorant,T<:Number} =
    __ccolorview(C, A)  # necessary to avoid ambiguities from dispatch on eltype
__ccolorview(::Type{C}, A::AbstractArray{T}) where {T<:Number,C<:RGB{T}} = reinterpretc(C, A)
__ccolorview(::Type{C}, A::AbstractArray{T}) where {T<:Number,C<:AbstractRGB} =
    _colorview_reorder(C, A)
__ccolorview(::Type{C}, A::AbstractArray{T}) where {T<:Number,C<:Color{T}} = reinterpretc(C, A)
__ccolorview(::Type{C}, A::AbstractArray{T}) where {T<:Number,C<:ColorAlpha} =
    _colorviewalpha(base_color_type(C), C, eltype(C), A)
__ccolorview(::Type{C}, A::AbstractArray{T}) where {T<:Number,C<:AlphaColor} =
    _colorview_reorder(C, A)
_colorviewalpha(::Type{C}, ::Type{CA}, ::Type{T}, A::AbstractArray{T}) where {C<:RGB,CA,T} =
    reinterpretc(CA, A)
_colorviewalpha(::Type{C}, ::Type{CA}, ::Type{T}, A::AbstractArray{T}) where {C<:AbstractRGB,CA,T} =
    _colorview_reorder(CA, A)
_colorviewalpha(::Type{C}, ::Type{CA}, ::Type{T}, A::AbstractArray{T}) where {C<:Color,CA,T} =
    reinterpretc(CA, A)

_colorview_reorder(::Type{C}, A) where C = reinterpretc(C, view(A, dimorder(C), Base.tail(colons(A))...))

colorview(::Type{ARGB32}, A::AbstractArray{BGRA{N0f8}}) = reinterpret(ARGB32, A)

colorview(::Type{C1}, A::AbstractArray{C2}) where {C1<:Colorant,C2<:Colorant} =
    colorview(C1, channelview(A))

colons(A::AbstractArray{T,N}) where {T,N} = ntuple(d->Colon(), Val(N))

"""
    colorview(C, gray1, gray2, ...) -> imgC

Combine numeric/grayscale images `gray1`, `gray2`, etc., into the
separate color channels of an array `imgC` with element type
`C<:Colorant`.

As a convenience, the constant `zeroarray` fills in an array of
matched size with all zeros.

# Example
```julia
imgC = colorview(RGB, r, zeroarray, b)
```

creates an image with `r` in the red chanel, `b` in the blue channel,
and nothing in the green channel.

See also: [`StackedView`](@ref).
"""
function colorview(::Type{C}, gray1, gray2, grays...) where C<:Colorant
    T = _colorview_type(eltype(C), promote_eleltype_all(gray1, gray2, grays...))
    CT = base_colorant_type(C){T}
    axs = firstinds(gray1, gray2, grays...)
    mappedarray(CT, extractchannels, take_zeros(eltype(CT), axs, gray1, gray2, grays...)...)
end

_colorview_type(::Type{Any}, ::Type{T}) where {T} = T
_colorview_type(::Type{T1}, ::Type{T2}) where {T1,T2} = T1

Base.@pure promote_eleltype_all(gray, grays...) = _promote_eleltype_all(beltype(eltype(gray)), grays...)
@inline function _promote_eleltype_all(::Type{T}, gray, grays...) where T
    _promote_eleltype_all(promote_type(T, beltype(eltype(gray))), grays...)
end
_promote_eleltype_all(::Type{T}) where {T} = T

beltype(::Type{T}) where {T} = eltype(T)
beltype(::Type{Union{}}) = Union{}

extractchannels(c::AbstractGray)    = (gray(c),)
extractchannels(c::TransparentGray) = (gray(c), alpha(c))
extractchannels(c::Color3)          = (comp1(c), comp2(c), comp3(c))
extractchannels(c::Transparent3)    = (comp1(c), comp2(c), comp3(c), alpha(c))

## Tuple & indexing utilities

_size(A::AbstractArray) = map(length, axes(A))

# color->number
@inline channelview_size(parent::AbstractArray{C}) where {C<:Colorant} = (length(C), _size(parent)...)
@inline channelview_axes(parent::AbstractArray{C}) where {C<:Colorant} =
    _cvi(Base.OneTo(length(C)), axes(parent))
_cvi(rc, ::Tuple{}) = (rc,)
_cvi(rc, inds::Tuple{R,Vararg{R}}) where {R<:AbstractUnitRange} = (convert(R, rc), inds...)
@inline channelview_size(parent::AbstractArray{C}) where {C<:Color1} = _size(parent)
@inline channelview_axes(parent::AbstractArray{C}) where {C<:Color1} = axes(parent)

function check_ncolorchan(::AbstractArray{C}, dims) where C<:Colorant
    dims[1] == length(C) || throw(DimensionMismatch("new array has $(dims[1]) color channels, must have $(length(C))"))
end
chanparentsize(::AbstractArray{C}, dims) where {C<:Colorant} = tail(dims)
@inline colparentsize(::Type{C}, dims) where {C<:Colorant} = (length(C), dims...)

channelview_dims_offset(parent::AbstractArray{C}) where {C<:Colorant} = 1

check_ncolorchan(::AbstractArray{C}, dims) where {C<:Color1} = nothing
chanparentsize(::AbstractArray{C}, dims) where {C<:Color1} = dims
colparentsize(::Type{C}, dims) where {C<:Color1} = dims
channelview_dims_offset(parent::AbstractArray{C}) where {C<:Color1} = 0

@inline indexsplit(A::AbstractArray{C}, I) where {C<:Colorant} = I[1], tail(I)
@inline indexsplit(A::AbstractArray{C}, I) where {C<:Color1} = 1, I

# number->color
@inline colorview_size(::Type{C}, parent::AbstractArray) where {C<:Colorant} = tail(_size(parent))
@inline colorview_axes(::Type{C}, parent::AbstractArray) where {C<:Colorant} = tail(axes(parent))
@inline colorview_size(::Type{C}, parent::AbstractArray) where {C<:Color1} = _size(parent)
@inline colorview_axes(::Type{C}, parent::AbstractArray) where {C<:Color1} = axes(parent)

function checkdim1(::Type{C}, inds) where C<:Colorant
    inds[1] == (1:length(C)) || throw(DimensionMismatch("dimension 1 must have indices 1:$(length(C)), got $(inds[1])"))
    nothing
end
checkdim1(::Type{C}, dims) where {C<:Color1} = nothing

parentaxes(::Type, inds) = tail(inds)
parentaxes(::Type{C}, inds) where {C<:Color1} = inds

celtype(::Type{Any}, ::Type{T}) where {T} = T
celtype(::Type{T1}, ::Type{T2}) where {T1,T2} = T1

## Low-level color utilities

tuplify(c::Color1) = (comp1(c),)
tuplify(c::Color3) = (comp1(c), comp2(c), comp3(c))
tuplify(c::Color2) = (comp1(c), alpha(c))
tuplify(c::Color4) = (comp1(c), comp2(c), comp3(c), alpha(c))

"""
    getchannels(P, C::Type, I)

Get a tuple of all channels needed to construct a Colorant of type `C`
from an `P::AbstractArray{<:Number}`.
"""
getchannels
@inline getchannels(P, ::Type{C}, I) where {C<:Color1} = (@inbounds ret = (P[I...],); ret)
@inline getchannels(P, ::Type{C}, I::Real) where {C<:Color1} = (@inbounds ret = (P[I],); ret)
@inline function getchannels(P, ::Type{C}, I) where C<:Color2
    @inbounds ret = (P[1,I...], P[2,I...])
    ret
end
@inline function getchannels(P, ::Type{C}, I) where C<:Color3
    @inbounds ret = (P[1,I...], P[2,I...],P[3,I...])
    ret
end
@inline function getchannels(P, ::Type{C}, I) where C<:Color4
    @inbounds ret = (P[1,I...], P[2,I...], P[3,I...], P[4,I...])
    ret
end

# setchannel (similar to setfield!)
# These don't check bounds since that's already done
"""
    setchannel(c, val, idx)

Equivalent to:

    cc = copy(c)
    cc[idx] = val
    cc

for immutable colors. `idx` is interpreted in the sense of constructor
arguments, so `setchannel(c, 0.5, 1)` would set red color channel for
any `c::AbstractRGB`, even if red isn't the first field in the type.
"""
setchannel(c::Colorant{T,1}, val, Ic::Int) where {T} = typeof(c)(val)

setchannel(c::TransparentColor{C,T,2}, val, Ic::Int) where {C,T} =
    typeof(c)(ifelse(Ic==1,val,comp1(c)),
              ifelse(Ic==2,val,alpha(c)))

setchannel(c::Colorant{T,3}, val, Ic::Int) where {T} = typeof(c)(ifelse(Ic==1,val,comp1(c)),
                                                          ifelse(Ic==2,val,comp2(c)),
                                                          ifelse(Ic==3,val,comp3(c)))
setchannel(c::TransparentColor{C,T,4}, val, Ic::Int) where {C,T} =
    typeof(c)(ifelse(Ic==1,val,comp1(c)),
              ifelse(Ic==2,val,comp2(c)),
              ifelse(Ic==3,val,comp3(c)),
              ifelse(Ic==4,val,alpha(c)))

"""
    setchannels!(P, val, I)

For a color `val`, distribute its channels along `P[:, I...]` for
`P::AbstractArray{<:Number}`.
"""
setchannels!
@inline setchannels!(P, val::Color1, I) = (@inbounds P[I...] = comp1(val); val)
@inline function setchannels!(P, val::Color2, I)
    @inbounds P[1,I...] = comp1(val)
    @inbounds P[2,I...] = alpha(val)
    val
end
@inline function setchannels!(P, val::Color3, I)
    @inbounds P[1,I...] = comp1(val)
    @inbounds P[2,I...] = comp2(val)
    @inbounds P[3,I...] = comp3(val)
    val
end
@inline function setchannels!(P, val::Color4, I)
    @inbounds P[1,I...] = comp1(val)
    @inbounds P[2,I...] = comp2(val)
    @inbounds P[3,I...] = comp3(val)
    @inbounds P[4,I...] = alpha(val)
    val
end


text:
WIDTH=600
HEIGHT=500
BACKGROUND=colorant"#c54245"
tree_green=colorant"#45c542"

tree_lines=[((340,300),(340,260)),
((340,260),(500,260)),
((5000,260),(340,150)),
((340,150),(450,150)),
((450,150),(340,075)),
((340,075),(400,075)),
((400,075),(300,000)),
((300,000),(200,075)),
((200,075),(260,075)),
((260,075),(150,150)),
((150,150),(260,150)),
((260,150),(100,260)),
((100,260),(260,260)),
((260,260),(260,300)),
((260,300),(340,300))
]
word_list=readlines(joinpath(@__DIR__,"words.txt"))
word_list=strip.(word_list)
word=rand(word_list)
word_array=string.(collect(word))
num_chars=length(word)
answers=fill(" ",num_chars)
letters=string.(collect('a':'z'))
score=15
state="playing"

function draw()
    for i in 1:(15-score)
        draw(Line(tree_lines[i][1]...,tree_lines[i][2]...),tree_green)
    end

    for i in 1:num_chars
        draw(Line(50i,450,50i+40,450),colorant"white")
        draw(Line(50i,451,50i+40,451),colorant"white")
        t=TextActor(answers[i],"moonhouse";pos=(50i+10,420))
        draw(t)
    end
    sc=TextActor("Tries: $score","moonhouse";pos=(450,20))
    draw(sc)

end

function on_key_down(g,k)
    global score
    key_pressed=lowercase(string(k))
    if key_pressed in letters && key_pressed in word_array
        @show "found $key_pressed in  word"
        answers[findall(x->x==key_pressed,word_array)] .= key_pressed
    else
        score =score-1
    end

    global state
    if score<1
        state="lost"
    end

    if !(" " in answers)
        state ="won"
    end

end

max_stars_repo_name: stepinski/julia-in

max_issues_repo_name: stepinski/julia-init

max_forks_repo_name: stepinski/julia-init


text:
@testset "Continuous space" begin

  # Basic model initialization
  space1 = ContinuousSpace(2; periodic = true, extend = (1, 1))
  space2 = ContinuousSpace(2; periodic = false, extend = (1, 1))
  space3 = ContinuousSpace(2; periodic = true, extend = (2, 1))
  space4 = ContinuousSpace(2; periodic = true, extend = (1, 2))
  space5 = ContinuousSpace(2; periodic = true, extend = (2, 2))
  space6 = ContinuousSpace(1; periodic = true, extend = (1,))
  space7 = ContinuousSpace(3; periodic = true, extend = (1,1,1))
  space8 = ContinuousSpace(3; periodic = false, extend = (1,1,1))

  @test space1.D == 2
  @test space2.D == 2
  @test space3.D == 2
  @test space4.D == 2
  @test space5.D == 2
  @test space6.D == 1
  @test space7.D == 3
  @test space8.D == 3

  model1 = ABM(Agent6, space1)
  model2 = ABM(Agent6, space2)
  model3 = ABM(Agent6, space3)
  model4 = ABM(Agent6, space4)
  model5 = ABM(Agent6, space5)
  model6 = ABM(Agent6, space6)

  @test nagents(model1) == 0
  @test Agents.collect_ids(DBInterface.execute(model1.space.db, "select id from tab")) == []

  # add_agent! with no existing agent (the agent is created)
  pos = (0.5, 0.5)
  vel = (0.2, 0.1)
  dia = 0.01
  add_agent!(pos, model1, vel, dia)
  @test Agents.collect_ids(DBInterface.execute(model1.space.db, "select id from tab")) == [1]
  dbrow = DBInterface.execute(model1.space.db, "select * from tab") |> DataFrame;
  @test dbrow[1, :a] == 0.5
  @test dbrow[1, :b] == 0.5

  # move_agent! without provided update_vel! function
  move_agent!(model1.agents[1], model1)
  dbrow = DBInterface.execute(model1.space.db, "select * from tab") |> DataFrame;
  @test dbrow[1, :a] == 0.7
  @test dbrow[1, :b] == 0.6
  @test dbrow[1, :a] == model1.agents[1].pos[1]
  @test dbrow[1, :b] == model1.agents[1].pos[2]

  kill_agent!(model1.agents[1], model1)
  dbrow = DBInterface.execute(model1.space.db, "select * from tab") |> DataFrame
  @test size(dbrow) == (0,0)

  # add_agent! with an existing agent
  agent = Agent6(2, pos, vel, dia)
  add_agent!(agent, model1)
  @test Agents.collect_ids(DBInterface.execute(model1.space.db, "select id from tab")) == [2]

  # agents within some range are found correctly (once this is implemented)
  agent2 = model1.agents[2]
  agent3 = Agent6(3, agent2.pos .+ 0.005, vel, dia)
  add_agent_pos!(agent3, model1)
  n_ids = space_neighbors(agent2, model1, agent2.weight)
  @test length(n_ids) == 1
  @test n_ids[1] == 3
  n_ids = space_neighbors(agent2, model1, agent2.weight/10)
  @test length(n_ids) == 0

  # test that it finds both
  n_ids = space_neighbors(agent2.pos, model1, agent2.weight)
  @test sort!(n_ids) == [2, 3]
end


text:
using OffsetArrays

setprecision(2000)
T = BigFloat

function rs_find_λ(T, n)
    # https://oeis.org/A087617
    λ = OffsetArray(zeros(T, n + 1), 0:n)
    λ[0] = 1
    eulers = rs_find_eulers(T, n)
    for i in 0:n-1
        c = 0
        for k in 0:i
            c += 2^(4k + 1) * eulers[k + 1] * λ[i - k]
        end
        λ[i + 1] = c / (i + 1)
    end

    λ
end

function rs_find_eulers(T, n)
    # https://oeis.org/A000364
    a = OffsetArray(zeros(T, n + 1), 0:n)
    a[0] = 1
    for i in 1:n
        for j in 0:i-1
            a[i] -= (-1)^(i + j) * a[j] * T(binomial(BigInt(2 * i), BigInt(2 * j)))
        end
    end
    a
end

function rs_find_d(T, m)
    d = OffsetArray(zeros(T, m + 1, m + 1), 0:m, 0:m)
    λ = rs_find_λ(T, m)
    for i in 0:m÷4
        d[i * 4, i * 3] = λ[i]
    end
    for n in 0:m-1
        for k in 0:(3n+2)÷4  # k < 3(n+1)/4  <=> k <= (3n+2)/4
            d[n + 1, k] = (3n + 1 - 4k) * (3n + 2 - 4k) * d[n, k] + (k >= 1 ? d[n, k - 1] : 0)
        end
    end

    d
end

function get_last_nonzero(ps)
    n = length(ps) - 1
    for i in n:-1:0
        if ps[i] != 0
            return i
        end
    end
    0
end

function ps_mul(a, b)
    n = length(a) - 1
    c = OffsetArray(zeros(typeof(a[0]), n+1), 0:n)
    max_i = get_last_nonzero(a)
    max_j = get_last_nonzero(b)
    for i in 0:min(n, max_i)
        for j in 0:min(n-i, max_j)
            c[i + j] += a[i] * b[j]
        end
    end
    c
end

function ps_div(a, b)
    n = length(a) - 1
    c = OffsetArray(zeros(typeof(a[0]), n+1), 0:n)
    for i in 0:n
        v = a[i]
        for j in 0:i-1
            v = v - b[i - j] * c[j]
        end
        c[i] = v / b[0]
    end
    c
end

function ps_compose(coeffs, ps)
    x = ps[0]
    n = length(ps) - 1
    a = copy(ps)
    a[0] = 0

    b = OffsetArray(zeros(typeof(a[0]), n+1), 0:n)
    for i in n:-1:0
        b = ps_mul(b, a)
        b[0] += coeffs[i]
    end
    b
end

function ps_apply(derivatives, ps)
    x = ps[0]
    n = length(ps) - 1
    a = copy(ps)
    a[0] = 0

    b = OffsetArray(zeros(typeof(a[0]), n+1), 0:n)
    factorial = typeof(a[0])(1.0)
    for i in 0:n
        derivatives[i] /= factorial
        factorial *= i + 1
    end
    for i in n:-1:0
        b = ps_mul(b, a)
        b[0] += derivatives[i]
    end
    b
end

function ps_cos(ps)
    n = length(ps) - 1
    x = ps[0]

    cos_x = cos(x)
    sin_x = sin(x)

    derivatives = [-sin_x, -cos_x, sin_x, cos_x]
    fn = OffsetArray(zeros(typeof(x), n+1), 0:n)
    for i in 0:n
        fn[i] = derivatives[mod1(i, 4)]
    end
    fn

    ps_apply(fn, ps)
end

function rs_find_ps(z, n)
    @assert n >= 2

    T = typeof(z)
    numer = OffsetArray(zeros(T, n + 1), 0:n)
    n_pi = T(π)
    numer[0] = z^2 * n_pi / 2 + 3 * n_pi / 8
    numer[1] = n_pi * z
    numer[2] = n_pi / 2
    numer = ps_cos(numer)

    denom = OffsetArray(zeros(T, n + 1), 0:n)
    denom[0] = n_pi * z
    denom[1] = n_pi
    denom = ps_cos(denom)

    ps = ps_div(numer, denom)
#     println(ps_mul(denom, ps) - numer)
    ps
end

function diff_ps(ps)
    n = length(ps) - 1
    ps = copy(ps)
    for i in 1:n
        ps[i - 1] = ps[i] * i
    end
    ps[n] = 0
    ps
end

n = 10
m = 1000

d = rs_find_d(T, n)
ps = rs_find_ps(T(0), m)
derivatives = OffsetArray(zeros(T, 3n+1, m + 1), 0:3n, 0:m)
derivatives[0, :] = ps
for i in 1:3n
    derivatives[i, :] = diff_ps(derivatives[i-1, :])
end

C = OffsetArray(zeros(T, n+1, m + 1), 0:n, 0:m)
for j in 0:n
    for k in 0:(3j ÷ 4)
        C[j, :] = C[j, :] + derivatives[3j - 4k, :] / factorial(T(3j - 4k)) * d[j, k] / T(π)^(2j - 2k)
    end
    C[j, :] /= T(2)^(2j)
end

fac = OffsetArray(zeros(T, 1001), 0:1000)
for i in 0:1000
    fac[i] = factorial(T(i))
end

function binom(n::Int, m::Int)
    fac[n] / fac[m] / fac[n - m]
end

function transform_to_chebyshev(ps) # transform a power series to Chebyshev's
    γ = OffsetArray(zeros(T, m+1), 0:m)

    # even case
    if ps[0] == 0
        ps = OffsetArray([ps[1:m]; 0], 0:m)
    end
    for k in 0:m÷2
        for l in k:m÷2
            γ[2k] += binom(2l, l-k) * ps[2l] / T(2)^(2l - 1)
        end
    end
    γ
end

γ = OffsetArray(zeros(T, n+1, m+1), 0:n, 0:m)
for j in 0:n
    γ[j, :] = transform_to_chebyshev(C[j, :])
end

function repr(x::BigFloat)
    hi = Float64(x)
    lo = Float64(x - hi)
    "f64x2::new($hi, $lo)"
    # "$hi"
end

output_m = 50
println("[")
for i in 0:n
    println("  [")
    for j in 0:2:output_m
        println("    $(repr(j == 0 ? γ[i, 0] / 2 : γ[i, j])),")
    end
    println("  ],")
end
println("]")

# println("[")
# for i in 0:n
#     println("  [")
#     for j in 0:2:output_m
#         println("    $(repr(j == 0 ? γ[i, 0] / 2 : γ[i, j])),")
#     end
#     println("  ], ")
# end
# println("]")

max_stars_repo_name: roosephu/apc

max_issues_repo_name: roosephu/apc

max_forks_repo_name: roosephu/apc


text:
module TunabilityOscillations

using Catalyst
using DifferentialEquations
using Statistics
using DSP
using Peaks
using LatinHypercubeSampling
using LSODA

export create_model, calc_main_freq, find_oscillations

include("create_model_equations.jl")
include("detect_oscillations.jl")

end


text:
@testset "6.6.7 (d hyper)^m (a+b (c csch)^n)^p" begin
    (a, b, c, d, x, ) = @variables a b c d x

    #= ::Package:: =#

    #= ::Section:: =#
    #=Integrands*of*the*form*Hyper(c+d*x)^m*(a+b*csch(c+d*x)^2)^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sinh(c+d*x)^m*(a+b*csch(c+d*x)^2)^n=#

    #= ::Subsubsection::Closed:: =#
    #=n=#

    @test_int [(a + b*csch(c + d*x)^2)^4, x, 4, a^4*x - ((2*a - b)*b*(2*a^2 - 2*a*b + b^2)*coth(c + d*x))/d - (b^2*(6*a^2 - 8*a*b + 3*b^2)*coth(c + d*x)^3)/(3*d) - ((4*a - 3*b)*b^3*coth(c + d*x)^5)/(5*d) - (b^4*coth(c + d*x)^7)/(7*d)]
    @test_int [(a + b*csch(c + d*x)^2)^3, x, 4, a^3*x - (b*(3*a^2 - 3*a*b + b^2)*coth(c + d*x))/d - ((3*a - 2*b)*b^2*coth(c + d*x)^3)/(3*d) - (b^3*coth(c + d*x)^5)/(5*d)]
    @test_int [(a + b*csch(c + d*x)^2)^2, x, 4, a^2*x - ((2*a - b)*b*coth(c + d*x))/d - (b^2*coth(c + d*x)^3)/(3*d)]
    @test_int [(a + b*csch(c + d*x)^2)^1, x, 3, a*x - (b*coth(c + d*x))/d]
    @test_int [1/(a + b*csch(c + d*x)^2)^1, x, 3, x/a - (sqrt(b)*atan((sqrt(a - b)*tanh(c + d*x))/sqrt(b)))/(a*sqrt(a - b)*d)]
    @test_int [1/(a + b*csch(c + d*x)^2)^2, x, 5, x/a^2 - ((3*a - 2*b)*sqrt(b)*atan((sqrt(a - b)*tanh(c + d*x))/sqrt(b)))/(2*a^2*(a - b)^(3/2)*d) + (b*coth(c + d*x))/(2*a*(a - b)*d*(a - b + b*coth(c + d*x)^2))]
    @test_int [1/(a + b*csch(c + d*x)^2)^3, x, 6, x/a^3 - (sqrt(b)*(15*a^2 - 20*a*b + 8*b^2)*atan((sqrt(a - b)*tanh(c + d*x))/sqrt(b)))/(8*a^3*(a - b)^(5/2)*d) + (b*coth(c + d*x))/(4*a*(a - b)*d*(a - b + b*coth(c + d*x)^2)^2) + ((7*a - 4*b)*b*coth(c + d*x))/(8*a^2*(a - b)^2*d*(a - b + b*coth(c + d*x)^2))]
    @test_int [1/(a + b*csch(c + d*x)^2)^4, x, 7, x/a^4 - (sqrt(b)*(35*a^3 - 70*a^2*b + 56*a*b^2 - 16*b^3)*atan((sqrt(a - b)*tanh(c + d*x))/sqrt(b)))/(16*a^4*(a - b)^(7/2)*d) + (b*coth(c + d*x))/(6*a*(a - b)*d*(a - b + b*coth(c + d*x)^2)^3) + ((11*a - 6*b)*b*coth(c + d*x))/(24*a^2*(a - b)^2*d*(a - b + b*coth(c + d*x)^2)^2) + (b*(19*a^2 - 22*a*b + 8*b^2)*coth(c + d*x))/(16*a^3*(a - b)^3*d*(a - b + b*coth(c + d*x)^2))]

    #= ::Subsubsection::Closed:: =#
    #=n/2=#

    @test_int [(a + b*csch(c + d*x)^2)^(5/2), x, 8, (a^(5/2)*atanh((sqrt(a)*coth(c + d*x))/sqrt(a - b + b*coth(c + d*x)^2)))/d - (sqrt(b)*(15*a^2 - 10*a*b + 3*b^2)*atanh((sqrt(b)*coth(c + d*x))/sqrt(a - b + b*coth(c + d*x)^2)))/(8*d) - ((7*a - 3*b)*b*coth(c + d*x)*sqrt(a - b + b*coth(c + d*x)^2))/(8*d) - (b*coth(c + d*x)*(a - b + b*coth(c + d*x)^2)^(3/2))/(4*d)]
    @test_int [(a + b*csch(c + d*x)^2)^(3/2), x, 7, (a^(3/2)*atanh((sqrt(a)*coth(c + d*x))/sqrt(a - b + b*coth(c + d*x)^2)))/d - ((3*a - b)*sqrt(b)*atanh((sqrt(b)*coth(c + d*x))/sqrt(a - b + b*coth(c + d*x)^2)))/(2*d) - (b*coth(c + d*x)*sqrt(a - b + b*coth(c + d*x)^2))/(2*d)]
    @test_int [(a + b*csch(c + d*x)^2)^(1/2), x, 6, (sqrt(a)*atanh((sqrt(a)*coth(c + d*x))/sqrt(a - b + b*coth(c + d*x)^2)))/d - (sqrt(b)*atanh((sqrt(b)*coth(c + d*x))/sqrt(a - b + b*coth(c + d*x)^2)))/d]
    @test_int [1/(a + b*csch(c + d*x)^2)^(1/2), x, 3, atanh((sqrt(a)*coth(c + d*x))/sqrt(a + b*csch(c + d*x)^2))/(sqrt(a)*d), atanh((sqrt(a)*coth(c + d*x))/sqrt(a - b + b*coth(c + d*x)^2))/(sqrt(a)*d)]
    @test_int [1/(a + b*csch(c + d*x)^2)^(3/2), x, 4, atanh((sqrt(a)*coth(c + d*x))/sqrt(a - b + b*coth(c + d*x)^2))/(a^(3/2)*d) + (b*coth(c + d*x))/(a*(a - b)*d*sqrt(a - b + b*coth(c + d*x)^2))]
    @test_int [1/(a + b*csch(c + d*x)^2)^(5/2), x, 6, atanh((sqrt(a)*coth(c + d*x))/sqrt(a - b + b*coth(c + d*x)^2))/(a^(5/2)*d) + (b*coth(c + d*x))/(3*a*(a - b)*d*(a - b + b*coth(c + d*x)^2)^(3/2)) + ((5*a - 3*b)*b*coth(c + d*x))/(3*a^2*(a - b)^2*d*sqrt(a - b + b*coth(c + d*x)^2))]
    @test_int [1/(a + b*csch(c + d*x)^2)^(7/2), x, 7, atanh((sqrt(a)*coth(c + d*x))/sqrt(a - b + b*coth(c + d*x)^2))/(a^(7/2)*d) + (b*coth(c + d*x))/(5*a*(a - b)*d*(a - b + b*coth(c + d*x)^2)^(5/2)) + ((9*a - 5*b)*b*coth(c + d*x))/(15*a^2*(a - b)^2*d*(a - b + b*coth(c + d*x)^2)^(3/2)) + (b*(33*a^2 - 40*a*b + 15*b^2)*coth(c + d*x))/(15*a^3*(a - b)^3*d*sqrt(a - b + b*coth(c + d*x)^2))]

    @test_int [(1 + csch(x)^2)^(3/2), x, 4, (-(1/2))*(coth(x)^2)^(3/2)*tanh(x) + sqrt(coth(x)^2)*log(sinh(x))*tanh(x)]
    @test_int [sqrt(1 + csch(x)^2), x, 3, sqrt(coth(x)^2)*log(sinh(x))*tanh(x)]
    @test_int [1/sqrt(1 + csch(x)^2), x, 3, (coth(x)*log(cosh(x)))/sqrt(coth(x)^2)]

    @test_int [(1 - csch(x)^2)^(3/2), x, 6, 2*asin(coth(x)/sqrt(2)) + atanh(coth(x)/sqrt(2 - coth(x)^2)) + (1/2)*coth(x)*sqrt(2 - coth(x)^2)]
    @test_int [sqrt(1 - csch(x)^2), x, 5, asin(coth(x)/sqrt(2)) + atanh(coth(x)/sqrt(2 - coth(x)^2))]
    @test_int [1/sqrt(1 - csch(x)^2), x, 3, atanh(coth(x)/sqrt(2 - coth(x)^2))]

    @test_int [(-1 + csch(x)^2)^(3/2), x, 7, atan(coth(x)/sqrt(-2 + coth(x)^2)) + 2*atanh(coth(x)/sqrt(-2 + coth(x)^2)) - (1/2)*coth(x)*sqrt(-2 + coth(x)^2)]
    @test_int [sqrt(-1 + csch(x)^2), x, 6, -atan(coth(x)/sqrt(-2 + coth(x)^2)) - atanh(coth(x)/sqrt(-2 + coth(x)^2))]
    @test_int [1/sqrt(-1 + csch(x)^2), x, 3, atan(coth(x)/sqrt(-2 + coth(x)^2))]

    @test_int [(-1 - csch(x)^2)^(3/2), x, 4, (1/2)*coth(x)*sqrt(-coth(x)^2) - sqrt(-coth(x)^2)*log(sinh(x))*tanh(x)]
    @test_int [sqrt(-1 - csch(x)^2), x, 3, sqrt(-coth(x)^2)*log(sinh(x))*tanh(x)]
    @test_int [1/sqrt(-1 - csch(x)^2), x, 3, (coth(x)*log(cosh(x)))/sqrt(-coth(x)^2)]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*cosh(c+d*x)^m*(a+b*csch(c+d*x)^2)^n=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*tanh(c+d*x)^m*(a+b*csch(c+d*x)^2)^n=#
    end

max_stars_repo_path: test/6 Hyperbolic functions/6.6 Hyperbolic cosecant/6.6.7 (d hyper)^m (a+b (c csch)^n)^p.jl


text:
export SG_Filter, maxDerivativeOrder, polynomialOrder, apply_SG_filter, apply_SG_filter2D

import Base: filter, length

function _Vandermonde(T::DataType=Float64;halfWidth::Int=5,degree::Int=2)::Array{T,2}

    @assert halfWidth>=0
    @assert degree>=0

    x=T[i for i in -halfWidth:halfWidth]

    n = degree+1
    m = length(x)

    V = Array{T}(undef,m, n)

    for i = 1:m
        V[i,1] = T(1)
    end
    for j = 2:n
        for i = 1:m
            V[i,j] = x[i] * V[i,j-1]
        end
    end

    return V
end

#+SG_Filters
#
# A structure to store Savitzky-Golay filters.
#
struct SG_Filter{T<:AbstractFloat,N}
    _filter_set::Array{LinearFilter_DefaultCentered{T,N},1}
end
#+SG_Filters
#
# Returns the filter to be used to compute the  smoothed derivatives of order *derivativeOrder*.
#
function filter(sg::SG_Filter{T,N};derivativeOrder::Int=0) where {T<:AbstractFloat,N}
    @assert 0<= derivativeOrder <= maxDerivativeOrder(sg)
    return sg._filter_set[derivativeOrder+1]
end
#+SG_Filters
#
# Returns filter length, this is an odd number, see [[SG_Filters_Constructor][]]
length(sg::SG_Filter{T,N}) where {T<:AbstractFloat,N} = length(filter(sg))
#+SG_Filters
#
# Maximum order of the smoothed derivatives we can compute with *sg*
#
maxDerivativeOrder(sg::SG_Filter{T,N}) where {T<:AbstractFloat,N} = size(sg._filter_set,1)-1
#+SG_Filters
# Returns the degree of the polynomial used to construct the Savitzky-Golay filters, see [[SG_Filters_Constructor][]].
polynomialOrder(sg::SG_Filter{T,N}) where {T<:AbstractFloat,N} = maxDerivativeOrder(sg)

#+SG_Filters L:SG_Filters_Constructor
#
# Creates a set of Savitzky-Golay filters
#
# - filter length is 2*halfWidth+1
# - polynomial degree is degree
#
function SG_Filter(T::DataType=Float64;halfWidth::Int=5,degree::Int=2)::SG_Filter
    @assert degree>=0
    @assert halfWidth>=1
    @assert 2*halfWidth>degree

    V=_Vandermonde(T,halfWidth=halfWidth,degree=degree)
    Q,R=qr(V)
    # breaking change in Julia V1.0,
    # see https://github.com/JuliaLang/julia/issues/27397
    #
    # before Q was a "plain" matrix, now stored in compact form
    #
    # SG=R\Q'
    #
    # must be replaced by
    #
    # Q=Q*Matrix(I, size(V))
    # SG=R\Q'
    #
    Q=Q*Matrix(I, size(V))
    SG=R\Q'

    n_filter,n_coef = size(SG)

    buffer=Array{LinearFilter_DefaultCentered{T,n_coef},1}()

    for i in 1:n_filter
        SG[i,:]*=factorial(i-1)
        push!(buffer,LinearFilter(SG[i,:]))
    end

# Returns filters set
    return SG_Filter{T,n_coef}(buffer)
end

# +SG_Filters
#
# Applies SG filter to 1D signal
#
# *Returns:*
# - the smoothed signal
#
function apply_SG_filter(signal::Array{T,1},
                         sg::SG_Filter{T};
                         derivativeOrder::Int=0,
                         left_BE::Type{<:BoundaryExtension}=ConstantBE,
                         right_BE::Type{<:BoundaryExtension}=ConstantBE) where {T<:AbstractFloat}

    return directCrossCorrelation(filter(sg,derivativeOrder=derivativeOrder),
                                  signal,
                                  left_BE,
                                  right_BE)
end

# +SG_Filters
#
# Applies SG filter to 2D signal
#
# *Returns:*
# - the smoothed signal
#
function apply_SG_filter2D(signal::Array{T,2},
                           sg_I::SG_Filter{T},
                           sg_J::SG_Filter{T};
                           derivativeOrder_I::Int=0,
                           derivativeOrder_J::Int=0,
                           min_I_BE::Type{<:BoundaryExtension}=ConstantBE,
                           max_I_BE::Type{<:BoundaryExtension}=ConstantBE,
                           min_J_BE::Type{<:BoundaryExtension}=ConstantBE,
                           max_J_BE::Type{<:BoundaryExtension}=ConstantBE) where {T<:AbstractFloat}

    return directCrossCorrelation2D(filter(sg_I,derivativeOrder=derivativeOrder_I),
                                    filter(sg_J,derivativeOrder=derivativeOrder_J),
                                    signal,
                                    min_I_BE,
                                    max_I_BE,
                                    min_J_BE,
                                    max_J_BE)
end


text:
using BasisChange
using Documenter

DocMeta.setdocmeta!(BasisChange, :DocTestSetup, :(using BasisChange); recursive=true)

makedocs(;
    modules=[BasisChange],
    authors="Sebastian Micluța-Câmpeanu <m.c.sebastian95@gmail.com>, Petru-Vlad Toma <tomapv@gmail.com>",
    repo="https://github.com/ctp-fpub/BasisChange.jl/blob/{commit}{path}#{line}",
    sitename="BasisChange.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ctp-fpub.github.io/BasisChange.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/ctp-fpub/BasisChange.jl",
)


text:
# user functions

"""
    identitypose6fg()

Initialize a and exiting or new factor graph with a first pose and prior as specified by default
keywords or user.
"""
function identitypose6fg!(;fg=nothing,
      N::Int=100,
      initCov::Array{Float64,2}=0.001*eye(6),
      initpose::SE3=SE3(0) )
  #
  fg = fg != nothing ? fg : initfg()
  N = 100

  println("Adding PriorPose3 and :x0 to graph...")
  # pts = 0.001*randn(6,N)
  initPosePrior = PriorPose3( MvNormal(veeEuler(initpose), initCov)  )
  pts = getSample(initPosePrior,N)[1]
  v0 = addVariable!(fg, :x0,  pts,  N=N, tags=[:POSE;])
  f1  = addFactor!(fg,[v0], initPosePrior)
  return fg
end

function identitypose6fg(;fg=nothing,
      N::Int=100,
      initCov::Array{Float64,2}=0.001*eye(6),
      initpose::SE3=SE3(0) )
   @warn "Caesar.identitypose6fg deprecated, use identitypose6fg! instead"
   identitypose6fg!(fg=fg, N=N, initCov=initCov, initpose=initpose)
end

function hasval(d::Dict, va)
  for (k,v) in d
    if v == va
      return true
    end
  end
  return false
end

function projectrbe(fgl::G, from::Symbol, to::Vector{Float64}) where G <: AbstractDFG
  x1 = getKDEMax(getBelief(fgl,from))
  x2 = to

  wTc = SE3(x1[1:3],Euler(x1[4:6]...))

  cX2 = (matrix(wTc) \ [x2[1:3];1.0])[1:3]

  # Dx = x2[1:3]-x1[1:3]
  Dx = cX2
  range = norm(Dx)
  bearing = atan2(Dx[2],Dx[1])
  elev = -atan2(Dx[3],norm(Dx[1:2]))
  return range, bearing, elev
end

function projectrbe(fgl::G, from::Symbol, to::Symbol) where G <: AbstractDFG
  # x1 = getKDEMax(getBelief(fg,from))
  x2 = getKDEMax(getBelief(fgl,to))
  return projectrbe(fgl, from, x2[1:3])

  # wTc = SE3(x1[1:3],Euler(x1[4:6]...))
  #
  # cX2 = (matrix(wTc) \ [x2;1.0])[1:3]
  #
  # # Dx = x2[1:3]-x1[1:3]
  # Dx = cX2
  # range = norm(Dx)
  # bearing = atan2(Dx[2],Dx[1])
  # elev = -atan2(Dx[3],norm(Dx[1:2]))
  # return range, bearing, elev
end

#