text:
# Julia wrapper for header: /usr/include/libavutil/opt.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_opt_show2,
    av_opt_set_defaults,
    av_opt_set_defaults2,
    av_set_options_string,
    av_opt_set_from_string,
    av_opt_free,
    av_opt_flag_is_set,
    av_opt_set_dict,
    av_opt_set_dict2,
    av_opt_get_key_value,
    av_opt_eval_flags,
    av_opt_eval_int,
    av_opt_eval_int64,
    av_opt_eval_float,
    av_opt_eval_double,
    av_opt_eval_q,
    av_opt_find,
    av_opt_find2,
    av_opt_next,
    av_opt_child_next,
    av_opt_child_class_next,
    av_opt_set,
    av_opt_set_int,
    av_opt_set_double,
    av_opt_set_q,
    av_opt_set_bin,
    av_opt_set_image_size,
    av_opt_set_pixel_fmt,
    av_opt_set_sample_fmt,
    av_opt_set_video_rate,
    av_opt_set_channel_layout,
    av_opt_set_dict_val,
    av_opt_get,
    av_opt_get_int,
    av_opt_get_double,
    av_opt_get_q,
    av_opt_get_image_size,
    av_opt_get_pixel_fmt,
    av_opt_get_sample_fmt,
    av_opt_get_video_rate,
    av_opt_get_channel_layout,
    av_opt_get_dict_val,
    av_opt_ptr,
    av_opt_freep_ranges,
    av_opt_query_ranges,
    av_opt_copy,
    av_opt_query_ranges_default,
    av_opt_is_set_to_default,
    av_opt_is_set_to_default_by_name,
    av_opt_serialize

function av_opt_show2(obj, av_log_obj, req_flags::Integer, rej_flags::Integer)
    ccall((:av_opt_show2, libavutil), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Cint), obj, av_log_obj, req_flags, rej_flags)
end

function av_opt_set_defaults(s)
    ccall((:av_opt_set_defaults, libavutil), Cvoid, (Ptr{Cvoid},), s)
end

function av_opt_set_defaults2(s, mask::Integer, flags::Integer)
    ccall((:av_opt_set_defaults2, libavutil), Cvoid, (Ptr{Cvoid}, Cint, Cint), s, mask, flags)
end

function av_set_options_string(ctx, opts, key_val_sep, pairs_sep)
    ccall((:av_set_options_string, libavutil), Cint, (Ptr{Cvoid}, Cstring, Cstring, Cstring), ctx, opts, key_val_sep, pairs_sep)
end

function av_opt_set_from_string(ctx, opts, shorthand, key_val_sep, pairs_sep)
    ccall((:av_opt_set_from_string, libavutil), Cint, (Ptr{Cvoid}, Cstring, Ptr{Cstring}, Cstring, Cstring), ctx, opts, shorthand, key_val_sep, pairs_sep)
end

function av_opt_free(obj)
    ccall((:av_opt_free, libavutil), Cvoid, (Ptr{Cvoid},), obj)
end

function av_opt_flag_is_set(obj, field_name, flag_name)
    ccall((:av_opt_flag_is_set, libavutil), Cint, (Ptr{Cvoid}, Cstring, Cstring), obj, field_name, flag_name)
end

function av_opt_set_dict(obj, options)
    ccall((:av_opt_set_dict, libavutil), Cint, (Ptr{Cvoid}, Ref{Ptr{AVDictionary}}), obj, options)
end

function av_opt_set_dict2(obj, options, search_flags::Integer)
    ccall((:av_opt_set_dict2, libavutil), Cint, (Ptr{Cvoid}, Ref{Ptr{AVDictionary}}, Cint), obj, options, search_flags)
end

function av_opt_get_key_value(ropts, key_val_sep, pairs_sep, flags::Integer, rkey, rval)
    ccall((:av_opt_get_key_value, libavutil), Cint, (Ptr{Cstring}, Cstring, Cstring, UInt32, Ptr{Cstring}, Ptr{Cstring}), ropts, key_val_sep, pairs_sep, flags, rkey, rval)
end

function av_opt_eval_flags(obj, o, val, flags_out)
    ccall((:av_opt_eval_flags, libavutil), Cint, (Ptr{Cvoid}, Ptr{AVOption}, Cstring, Ptr{Cint}), obj, o, val, flags_out)
end

function av_opt_eval_int(obj, o, val, int_out)
    ccall((:av_opt_eval_int, libavutil), Cint, (Ptr{Cvoid}, Ptr{AVOption}, Cstring, Ptr{Cint}), obj, o, val, int_out)
end

function av_opt_eval_int64(obj, o, val, int64_out)
    ccall((:av_opt_eval_int64, libavutil), Cint, (Ptr{Cvoid}, Ptr{AVOption}, Cstring, Ptr{Int64}), obj, o, val, int64_out)
end

function av_opt_eval_float(obj, o, val, float_out)
    ccall((:av_opt_eval_float, libavutil), Cint, (Ptr{Cvoid}, Ptr{AVOption}, Cstring, Ptr{Cfloat}), obj, o, val, float_out)
end

function av_opt_eval_double(obj, o, val, double_out)
    ccall((:av_opt_eval_double, libavutil), Cint, (Ptr{Cvoid}, Ptr{AVOption}, Cstring, Ptr{Cdouble}), obj, o, val, double_out)
end

function av_opt_eval_q(obj, o, val, q_out)
    ccall((:av_opt_eval_q, libavutil), Cint, (Ptr{Cvoid}, Ptr{AVOption}, Cstring, Ptr{AVRational}), obj, o, val, q_out)
end

function av_opt_find(obj, name, unit, opt_flags::Integer, search_flags::Integer)
    ccall((:av_opt_find, libavutil), Ptr{AVOption}, (Ptr{Cvoid}, Cstring, Cstring, Cint, Cint), obj, name, unit, opt_flags, search_flags)
end

function av_opt_find2(obj, name, unit, opt_flags::Integer, search_flags::Integer, target_obj)
    ccall((:av_opt_find2, libavutil), Ptr{AVOption}, (Ptr{Cvoid}, Cstring, Cstring, Cint, Cint, Ptr{Ptr{Cvoid}}), obj, name, unit, opt_flags, search_flags, target_obj)
end

function av_opt_next(obj, prev)
    ccall((:av_opt_next, libavutil), Ptr{AVOption}, (Ptr{Cvoid}, Ptr{AVOption}), obj, prev)
end

function av_opt_child_next(obj, prev)
    ccall((:av_opt_child_next, libavutil), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}), obj, prev)
end

function av_opt_child_class_next(parent, prev)
    ccall((:av_opt_child_class_next, libavutil), Ptr{AVClass}, (Ptr{AVClass}, Ptr{AVClass}), parent, prev)
end

function av_opt_set(obj, name, val, search_flags::Integer)
    ccall((:av_opt_set, libavutil), Cint, (Ptr{Cvoid}, Cstring, Cstring, Cint), obj, name, val, search_flags)
end

function av_opt_set_int(obj, name, val::Int64, search_flags::Integer)
    ccall((:av_opt_set_int, libavutil), Cint, (Ptr{Cvoid}, Cstring, Int64, Cint), obj, name, val, search_flags)
end

function av_opt_set_double(obj, name, val::Cdouble, search_flags::Integer)
    ccall((:av_opt_set_double, libavutil), Cint, (Ptr{Cvoid}, Cstring, Cdouble, Cint), obj, name, val, search_flags)
end

function av_opt_set_q(obj, name, val::AVRational, search_flags::Integer)
    ccall((:av_opt_set_q, libavutil), Cint, (Ptr{Cvoid}, Cstring, AVRational, Cint), obj, name, val, search_flags)
end

function av_opt_set_bin(obj, name, val, size::Integer, search_flags::Integer)
    ccall((:av_opt_set_bin, libavutil), Cint, (Ptr{Cvoid}, Cstring, Ptr{UInt8}, Cint, Cint), obj, name, val, size, search_flags)
end

function av_opt_set_image_size(obj, name, w::Integer, h::Integer, search_flags::Integer)
    ccall((:av_opt_set_image_size, libavutil), Cint, (Ptr{Cvoid}, Cstring, Cint, Cint, Cint), obj, name, w, h, search_flags)
end

function av_opt_set_pixel_fmt(obj, name, fmt::AVPixelFormat, search_flags::Integer)
    ccall((:av_opt_set_pixel_fmt, libavutil), Cint, (Ptr{Cvoid}, Cstring, AVPixelFormat, Cint), obj, name, fmt, search_flags)
end

function av_opt_set_sample_fmt(obj, name, fmt::AVSampleFormat, search_flags::Integer)
    ccall((:av_opt_set_sample_fmt, libavutil), Cint, (Ptr{Cvoid}, Cstring, AVSampleFormat, Cint), obj, name, fmt, search_flags)
end

function av_opt_set_video_rate(obj, name, val::AVRational, search_flags::Integer)
    ccall((:av_opt_set_video_rate, libavutil), Cint, (Ptr{Cvoid}, Cstring, AVRational, Cint), obj, name, val, search_flags)
end

function av_opt_set_channel_layout(obj, name, ch_layout::Int64, search_flags::Integer)
    ccall((:av_opt_set_channel_layout, libavutil), Cint, (Ptr{Cvoid}, Cstring, Int64, Cint), obj, name, ch_layout, search_flags)
end

function av_opt_set_dict_val(obj, name, val, search_flags::Integer)
    ccall((:av_opt_set_dict_val, libavutil), Cint, (Ptr{Cvoid}, Cstring, Ptr{AVDictionary}, Cint), obj, name, val, search_flags)
end

function av_opt_get(obj, name, search_flags::Integer, out_val)
    ccall((:av_opt_get, libavutil), Cint, (Ptr{Cvoid}, Cstring, Cint, Ptr{Ptr{UInt8}}), obj, name, search_flags, out_val)
end

function av_opt_get_int(obj, name, search_flags::Integer, out_val)
    ccall((:av_opt_get_int, libavutil), Cint, (Ptr{Cvoid}, Cstring, Cint, Ptr{Int64}), obj, name, search_flags, out_val)
end

function av_opt_get_double(obj, name, search_flags::Integer, out_val)
    ccall((:av_opt_get_double, libavutil), Cint, (Ptr{Cvoid}, Cstring, Cint, Ptr{Cdouble}), obj, name, search_flags, out_val)
end

function av_opt_get_q(obj, name, search_flags::Integer, out_val)
    ccall((:av_opt_get_q, libavutil), Cint, (Ptr{Cvoid}, Cstring, Cint, Ptr{AVRational}), obj, name, search_flags, out_val)
end

function av_opt_get_image_size(obj, name, search_flags::Integer, w_out, h_out)
    ccall((:av_opt_get_image_size, libavutil), Cint, (Ptr{Cvoid}, Cstring, Cint, Ptr{Cint}, Ptr{Cint}), obj, name, search_flags, w_out, h_out)
end

function av_opt_get_pixel_fmt(obj, name, search_flags::Integer, out_fmt)
    ccall((:av_opt_get_pixel_fmt, libavutil), Cint, (Ptr{Cvoid}, Cstring, Cint, Ptr{AVPixelFormat}), obj, name, search_flags, out_fmt)
end

function av_opt_get_sample_fmt(obj, name, search_flags::Integer, out_fmt)
    ccall((:av_opt_get_sample_fmt, libavutil), Cint, (Ptr{Cvoid}, Cstring, Cint, Ptr{AVSampleFormat}), obj, name, search_flags, out_fmt)
end

function av_opt_get_video_rate(obj, name, search_flags::Integer, out_val)
    ccall((:av_opt_get_video_rate, libavutil), Cint, (Ptr{Cvoid}, Cstring, Cint, Ptr{AVRational}), obj, name, search_flags, out_val)
end

function av_opt_get_channel_layout(obj, name, search_flags::Integer, ch_layout)
    ccall((:av_opt_get_channel_layout, libavutil), Cint, (Ptr{Cvoid}, Cstring, Cint, Ptr{Int64}), obj, name, search_flags, ch_layout)
end

function av_opt_get_dict_val(obj, name, search_flags::Integer, out_val)
    ccall((:av_opt_get_dict_val, libavutil), Cint, (Ptr{Cvoid}, Cstring, Cint, Ref{Ptr{AVDictionary}}), obj, name, search_flags, out_val)
end

function av_opt_ptr(avclass, obj, name)
    ccall((:av_opt_ptr, libavutil), Ptr{Cvoid}, (Ptr{AVClass}, Ptr{Cvoid}, Cstring), avclass, obj, name)
end

function av_opt_freep_ranges(ranges)
    ccall((:av_opt_freep_ranges, libavutil), Cvoid, (Ptr{Ptr{AVOptionRanges}},), ranges)
end

function av_opt_query_ranges(arg1, obj, key, flags::Integer)
    ccall((:av_opt_query_ranges, libavutil), Cint, (Ptr{Ptr{AVOptionRanges}}, Ptr{Cvoid}, Cstring, Cint), arg1, obj, key, flags)
end

function av_opt_copy(dest, src)
    ccall((:av_opt_copy, libavutil), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), dest, src)
end

function av_opt_query_ranges_default(arg1, obj, key, flags::Integer)
    ccall((:av_opt_query_ranges_default, libavutil), Cint, (Ptr{Ptr{AVOptionRanges}}, Ptr{Cvoid}, Cstring, Cint), arg1, obj, key, flags)
end

function av_opt_is_set_to_default(obj, o)
    ccall((:av_opt_is_set_to_default, libavutil), Cint, (Ptr{Cvoid}, Ptr{AVOption}), obj, o)
end

function av_opt_is_set_to_default_by_name(obj, name, search_flags::Integer)
    ccall((:av_opt_is_set_to_default_by_name, libavutil), Cint, (Ptr{Cvoid}, Cstring, Cint), obj, name, search_flags)
end

function av_opt_serialize(obj, opt_flags::Integer, flags::Integer, buffer, key_val_sep::UInt8, pairs_sep::UInt8)
    ccall((:av_opt_serialize, libavutil), Cint, (Ptr{Cvoid}, Cint, Cint, Ptr{Cstring}, UInt8, UInt8), obj, opt_flags, flags, buffer, key_val_sep, pairs_sep)
end


text:
@testset "Aut(Σ₃.₀)" begin
    genus = 3

    π₁Σ = Groups.SurfaceGroup(genus, 0)

    Groups.PermRightAut(p::Perm) = Groups.PermRightAut(p.d)
    # Groups.PermLeftAut(p::Perm) = Groups.PermLeftAut(p.d)
    autπ₁Σ = let autπ₁Σ = AutomorphismGroup(π₁Σ)
        pauts = let p = perm"(1,3,5)(2,4,6)"
            [Groups.PermRightAut(p^i) for i in 0:2]
        end
        T = eltype(KnuthBendix.letters(alphabet(autπ₁Σ)))
        S = eltype(pauts)

        A = Alphabet(Union{T,S}[KnuthBendix.letters(alphabet(autπ₁Σ)); pauts])

        autG = AutomorphismGroup(
            π₁Σ,
            autπ₁Σ.gens,
            A,
            ntuple(i->inv(gens(π₁Σ, i)), 2Groups.genus(π₁Σ))
        )

        autG
    end

    Al = alphabet(autπ₁Σ)
    S = [gens(autπ₁Σ); inv.(gens(autπ₁Σ))]

    sautFn = let ltrs = KnuthBendix.letters(Al)
        parent(first(ltrs).autFn_word)
    end

    τ = Groups.rotation_element(sautFn)

    @testset "Twists" begin
        A = KnuthBendix.alphabet(sautFn)
        λ = Groups.ΡΛ(:λ, A, 2genus)
        ϱ = Groups.ΡΛ(:ϱ, A, 2genus)
        @test sautFn(Groups.Te_diagonal(λ, ϱ, 1)) ==
            conj(sautFn(Groups.Te_diagonal(λ, ϱ, 2)), τ)

        @test sautFn(Groups.Te_diagonal(λ, ϱ, 3)) == sautFn(Groups.Te(λ, ϱ, 3, 1))
    end

    z = let d = Groups.domain(τ)
        Groups.evaluate(τ^genus)
    end

    @test π₁Σ.(word.(z)) == Groups.domain(first(S))
    d = Groups.domain(first(S))
    p = perm"(1,3,5)(2,4,6)"
    @test Groups.evaluate!(deepcopy(d), τ) == d^inv(p)
    @test Groups.evaluate!(deepcopy(d), τ^2) == d^p

    E, sizes = Groups.wlmetric_ball(S, radius=3)
    @test sizes == [49, 1813, 62971]
    B2 = @view E[1:sizes[2]]

    σ = autπ₁Σ(Word([Al[Groups.PermRightAut(p)]]))

    @test conj(S[7], σ) == S[10]
    @test conj(S[7], σ^2) == S[11]
    @test conj(S[9], σ) == S[12]
    @test conj(S[9], σ^2) == S[8]

    @test conj(S[1], σ) == S[4]
    @test conj(S[1], σ^2) == S[5]
    @test conj(S[3], σ) == S[6]
    @test conj(S[3], σ^2) == S[2]

    B2ᶜ = [conj(b, σ) for b in B2]
    @test B2ᶜ != B2
    @test Set(B2ᶜ) == Set(B2)
end

ext: jl


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using DensityInterface, ValueShapes
using ArraysOfArrays, Distributions, StatsBase, IntervalSets

@testset "truncated_density" begin
    prior_dist = unshaped(NamedTupleDist(
        a = truncated(Normal(), -2, 2),
        b = Exponential(),
        c = [1 2; 3 4],
        d = [-3..3, -4..4]
    ))
    prior = convert(AbstractDensity, prior_dist)

    likelihood = v -> (logval = 0,)

    posterior = PosteriorDensity(likelihood, prior)

    bounds = [-1..3, 1..4, -2..2, 0..4]

    @test (@inferred BAT.truncate_dist_hard(Exponential(), -1..5)).dist isa Truncated{<:Exponential}

    let
        orig_dist = Exponential()
        trunc_dist, logrenormf = @inferred BAT.truncate_dist_hard(orig_dist, -6..6)
        @test trunc_dist.lower == 0 && trunc_dist.upper == 6
        @test logpdf(trunc_dist, 3) + logrenormf ≈ logpdf(orig_dist, 3)
    end

    let
        orig_dist = truncated(Exponential(), 1, 5)
        trunc_dist, logrenormf = @inferred BAT.truncate_dist_hard(orig_dist, 0..4)
        @test trunc_dist.lower == 1 && trunc_dist.upper == 4
        @test logpdf(trunc_dist, 3) + logrenormf ≈ logpdf(orig_dist, 3)
    end

    let
        orig_dist = truncated(Exponential(), 1, 5)
        trunc_dist, logrenormf = @inferred BAT.truncate_dist_hard(orig_dist, 2..6)
        @test trunc_dist.lower == 2 && trunc_dist.upper == 5
        @test logpdf(trunc_dist, 3) + logrenormf ≈ logpdf(orig_dist, 3)
    end

    @test @inferred(BAT.truncate_dist_hard(prior_dist, bounds)).dist isa ValueShapes.UnshapedNTD
    let
        trunc_dist, logrenormf = @inferred BAT.truncate_dist_hard(prior_dist, bounds)
        @test logpdf(unshaped(trunc_dist), [1, 2, 0, 3]) + logrenormf ≈ logpdf(unshaped(prior_dist), [1, 2, 0, 3])
    end

    @test @inferred(truncate_density(prior, bounds)) isa BAT.RenormalizedDensity

    @test BAT.checked_logdensityof(unshaped(truncate_density(prior, bounds)), [1, 2, 0, 3]) ≈ BAT.checked_logdensityof(unshaped(prior), [1, 2, 0, 3])
    @test BAT.checked_logdensityof(truncate_density(prior, bounds), varshape(prior)([1, 2, 0, 3])) ≈ BAT.checked_logdensityof(prior, varshape(prior)([1, 2, 0, 3]))
    @test varshape(truncate_density(prior, bounds)) == varshape(prior)

    @test @inferred(truncate_density(posterior, bounds)) isa PosteriorDensity

    trunc_pstr = truncate_density(posterior, bounds)
    @test @inferred(BAT.checked_logdensityof(unshaped(trunc_pstr), [1, 2, 0, 3])) ≈ BAT.checked_logdensityof(unshaped(posterior), [1, 2, 0, 3])
    @test @inferred(BAT.checked_logdensityof(unshaped(trunc_pstr), [-1, -1, -1, -1])) ≈ -Inf
    @test varshape(trunc_pstr) == varshape(posterior)

    let
        trunc_prior_dist = parent(BAT.getprior(trunc_pstr)).dist
        s = bat_sample(trunc_pstr, MCMCSampling(mcalg = MetropolisHastings(), trafo = NoDensityTransform(), nsteps = 10^5)).result
        s_flat = flatview(unshaped.(s))
        @test all(minimum.(bounds) .< minimum(s_flat))
        @test all(maximum.(bounds) .> maximum(s_flat))
        cov_est = cov(unshaped.(s))
        @test isapprox(cov_est[1,1], var(trunc_prior_dist.shaped.a), rtol = 0.05)
        @test isapprox(cov_est[4,4], var(trunc_prior_dist.shaped.d.v[2]), rtol = 0.05)
    end
end


text:
type LSH_vars
    # LSHに用いる各種変数を格納しておくtype
    numhash::Int64
    numband::Int64
    bandwidth::Int64
    function LSH_vars(nb, bw)
        new(nb*bw, nb, bw)
    end
end

# 多重ディスパッチ
function minhash(words::Array{Int64, 1}, seed::UInt64)
    current_min = typemax(UInt64)
    minhash_word = nothing
    for word in words
        hash_ = hash(word, seed)
        if hash_ < current_min
            minhash_word = word
            current_min = hash_
        end
    end
    return(minhash_word)
end
function minhash(words::Set{Int64}, seed::UInt64)
    current_min = typemax(UInt64)
    minhash_word = nothing
    for word in words
        hash_ = hash(word, seed)
        if hash_ < current_min
            minhash_word = word
            current_min = hash_
        end
    end
    return(minhash_word)
end

function generate_signature_matrix(lsh_vars::LSH_vars, binarray_csc::SparseMatrixCSC{Bool,Int64}; window_num=100, slidewidth=50)
    numhash = lsh_vars.numhash
    extent = 1:slidewidth:(size(binarray_csc)[2] - 2window_num)
    signature_matrix = zeros((numhash, length(extent)))
    for row in 1:numhash
        for (idx, col) in enumerate(extent)
            # 発火が存在するIDのみからなるSetを作成する
            idsets = Set{Int64}(findn(binarray_csc[:, col:(col+window_num-1)] .== 1)[1])
            if length(idsets) > 0
                signature_matrix[row, idx] = minhash(idsets, UInt64(row))
            else
                # idsetsが空の場合 = 発火が存在しない場合も衝突を回避するためランダムな値を詰めておく
                signature_matrix[row, idx] = minhash(rand(1:typemax(Int64), 1), UInt64(row))
            end
         end
    end
    return(signature_matrix)
end

# b: numband, r:bandwidth(rows)
cal_prob(s, b, r) = 1 - ( 1- s^r )^b
function ret_prob_table(b, r, step=0.1)
    for s in 0:step:1
        println(s, ' ', cal_prob(s, b, r))
    end
end

function generate_buckets(lsh_vars::LSH_vars, signature_matrix::Array{Float64,2})
    # signature matrixの各列をバケツに放り込む関数
    # 似ているものは同じバケツに放り込まれやすい．具体的な確率はバンドの設計による
    numhash = lsh_vars.numhash
    numband = lsh_vars.numband
    bandwidth = lsh_vars.bandwidth
    buckets_list = []
    for band in 1:bandwidth:numhash
        buckets = Dict{UInt64, Set{Int64}}()
        for col in 1:size(signature_matrix)[2]
            hash_ = hash(signature_matrix[band:(band+bandwidth-1), col], UInt64(band))
            if !haskey(buckets, hash_)
                buckets[hash_] = Set{Int64}()
            end
            push!(buckets[hash_], col)
        end
        push!(buckets_list, buckets)
    end
    return(buckets_list)
end

function find_similar(lsh_vars::LSH_vars, signature_matrix::Array{Float64,2}, buckets_list::Array{Any,1}, col::Int64)
    # buckets_listの情報を元にcolに類似した（≃ Jaccard係数が高い）要素を抽出するための関数
    # Python版とは異なり，編集類似度を用いたスクリーニングはここでは行わない
    numhash = lsh_vars.numhash
    numband = lsh_vars.numband
    bandwidth = lsh_vars.bandwidth
    candidates = Set{Int64}()
    for (idx, band) in enumerate(1:numhash:bandwidth)
        hash_ = hash(signature_matrix[band:(band+bandwidth-1), col], UInt64(band))
        candidates = union(candidates, buckets_list[idx][hash_])
    end
    return(candidates)
end

max_stars_repo_name: KeitaW/MinLSH

max_issues_repo_name: KeitaW/MinLSH


text:
# Contiguous on-device arrays

export MtlDeviceArray, MtlDeviceVector, MtlDeviceMatrix

## construction

"""
    MtlDeviceArray(dims, ptr)
    MtlDeviceArray{T}(dims, ptr)
    MtlDeviceArray{T,A}(dims, ptr)
    MtlDeviceArray{T,A,N}(dims, ptr)

Construct an `N`-dimensional dense Metal device array with element type `T` wrapping a
pointer, where `N` is determined from the length of `dims` and `T` is determined from the
type of `ptr`.

`dims` may be a single scalar, or a tuple of integers corresponding to the
lengths in each dimension). If the rank `N` is supplied explicitly as in `Array{T,N}(dims)`,
then it must match the length of `dims`. The same applies to the element type `T`, which
should match the type of the pointer `ptr`.
"""
MtlDeviceArray

# NOTE: we can't support the typical `tuple or series of integer` style construction,
#       because we're Mtlrrently requiring a trailing pointer argument.

struct MtlDeviceArray{T,N,A} <: DenseArray{T,N}
    ptr::Core.LLVMPtr{T,A}
    shape::Dims{N}
    # inner constructors, fully parameterized, exact types (ie. Int not <:Integer)
    MtlDeviceArray{T,N,A}(shape::Dims{N}, ptr::Core.LLVMPtr{T,A}) where {T,A,N} = new(ptr,shape)
end

const MtlDeviceVector = MtlDeviceArray{T,1,A} where {T,A}
const MtlDeviceMatrix = MtlDeviceArray{T,2,A} where {T,A}

# outer constructors, non-parameterized
MtlDeviceArray(dims::NTuple{N,<:Integer}, p::Core.LLVMPtr{T,A})                where {T,A,N} = MtlDeviceArray{T,N,A}(dims, p)
MtlDeviceArray(len::Integer,              p::Core.LLVMPtr{T,A})                where {T,A}   = MtlDeviceVector{T,A}((len,), p)

# outer constructors, partially parameterized
MtlDeviceArray{T}(dims::NTuple{N,<:Integer},   p::Core.LLVMPtr{T,A}) where {T,A,N} = MtlDeviceArray{T,N,A}(dims, p)
MtlDeviceArray{T}(len::Integer,                p::Core.LLVMPtr{T,A}) where {T,A}   = MtlDeviceVector{T,A}((len,), p)
MtlDeviceArray{T,N}(dims::NTuple{N,<:Integer}, p::Core.LLVMPtr{T,A}) where {T,A,N} = MtlDeviceArray{T,N,A}(dims, p)
#MtlDeviceVector{T}(len::Integer,               p::Core.LLVMPtr{T,1}) where {T,A}   = MtlDeviceVector{T,A}((len,), p)

# outer constructors, fully parameterized
MtlDeviceArray{T,N,A}(dims::NTuple{N,<:Integer}, p::Core.LLVMPtr{T,A}) where {T,A,N} = MtlDeviceArray{T,N,A}(Int.(dims), p)
MtlDeviceVector{T,A}(len::Integer,               p::Core.LLVMPtr{T,A}) where {T,A}   = MtlDeviceVector{T,A}((Int(len),), p)
# MtlDeviceVector{T,A}(len::NTuple{N,<:Integer},               p::Core.LLVMPtr{T,1}) where {T,A}   = MtlDeviceVector{T,A}((Int(len),), p)

## getters

Base.pointer(a::MtlDeviceArray) = a.ptr
Base.pointer(a::MtlDeviceArray, i::Integer) = pointer(a) + (i - 1) * Base.elsize(a)

Base.elsize(::Type{<:MtlDeviceArray{T}}) where {T} = sizeof(T)
Base.size(g::MtlDeviceArray) = g.shape
# Testing to fix argument encoding with the trailing , for vectors
Base.size(g::MtlDeviceVector) = length(g)
Base.length(g::MtlDeviceArray) = prod(g.shape)

## conversions

Base.unsafe_convert(::Type{Core.LLVMPtr{T,A}}, a::MtlDeviceArray{T,N,A}) where {T,A,N} = pointer(a)

## indexing intrinsics

# NOTE: these intrinsics are now implemented using plain and simple pointer operations;
#       when adding support for isbits union arrays we will need to implement that here.

# FIXME: Bounscheck

@inline function arrayref(A::MtlDeviceArray{T}, index::Integer) where {T}
    #@boundscheck checkbounds(A, index)
    align = alignment(pointer(A))
    unsafe_load(pointer(A), index, Val(align))
end

@inline function arrayset(A::MtlDeviceArray{T}, x::T, index::Integer) where {T}
    #@boundscheck checkbounds(A, index)
    align = alignment(pointer(A))
    unsafe_store!(pointer(A), x, index, Val(align))
    return A
end

@inline function const_arrayref(A::MtlDeviceArray{T}, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    align = Base.datatype_alignment(T)
    unsafe_cached_load(pointer(A), index, Val(align))
end

## indexing

Base.@propagate_inbounds Base.getindex(A::MtlDeviceArray{T}, i1::Integer) where {T} =
    arrayref(A, i1)
Base.@propagate_inbounds Base.setindex!(A::MtlDeviceArray{T}, x, i1::Integer) where {T} =
    arrayset(A, convert(T,x)::T, i1)

# preserve the specific integer type when indexing device arrays,
# to avoid extending 32-bit hardware indices to 64-bit.
Base.to_index(::MtlDeviceArray, i::Integer) = i

# Base doesn't like Integer indices, so we need our own ND get and setindex! routines.
# See also: https://github.com/JuliaLang/julia/pull/42289
Base.@propagate_inbounds Base.getindex(A::MtlDeviceArray,
                                       I::Union{Integer, CartesianIndex}...) =
    A[Base._to_linear_index(A, to_indices(A, I)...)]
Base.@propagate_inbounds Base.setindex!(A::MtlDeviceArray, x,
                                        I::Union{Integer, CartesianIndex}...) =
    A[Base._to_linear_index(A, to_indices(A, I)...)] = x

# TODO: Put this in pointer_ptr.jl?
Base.@propagate_inbounds Base.getindex(A::Core.LLVMPtr{T}, i1::Integer) where {T} =
    arrayref(A, i1)
Base.setindex!(A::Core.LLVMPtr{T}, x, i1::Integer) where {T} =
    arrayset(A, convert(T,x)::T, i1)

Base.IndexStyle(::Type{<:Core.LLVMPtr}) = Base.IndexLinear()

@generated function alignment(::Core.LLVMPtr{T}) where {T}
    if Base.isbitsunion(T)
        _, sz, al = Base.uniontype_layout(T)
        al
    else
        Base.datatype_alignment(T)
    end
end

@inline function arrayref(A::Core.LLVMPtr{T,AS}, index::Integer) where {T,AS}
    #@boundscheck checkbounds(A, index)
    align = alignment(A)
    unsafe_load(A, index, Val(align))
end

@inline function arrayset(A::Core.LLVMPtr{T,AS}, x::T, index::Integer) where {T,AS}
    #@boundscheck checkbounds(A, index)
    align = alignment(A)
    unsafe_store!(A, x, index, Val(align))
    return A
end

## const indexing

"""
    Const(A::MtlDeviceArray)

Mark a MtlDeviceArray as constant/read-only and to use the constant address space.
!!! warning
    Experimental API. Subject to change without deprecation.
"""
struct Const{T,N,AS} <: DenseArray{T,N}
    a::MtlDeviceArray{T,N,AS}
end
Base.Experimental.Const(A::MtlDeviceArray) = Const(A)

Base.IndexStyle(::Type{<:Const}) = IndexLinear()
Base.size(C::Const) = size(C.a)
Base.axes(C::Const) = axes(C.a)
Base.@propagate_inbounds Base.getindex(A::Const, i1::Integer) = const_arrayref(A.a, i1)

## other

Base.show(io::IO, a::MtlDeviceVector) =
    print(io, "$(length(a))-element device array at $(pointer(a))")
Base.show(io::IO, a::MtlDeviceArray) =
    print(io, "$(join(a.shape, '×')) device array at $(pointer(a))")

Base.show(io::IO, mime::MIME"text/plain", a::MtlDeviceArray) = show(io, a)

@inline function Base.unsafe_view(A::MtlDeviceVector{T}, I::Vararg{Base.ViewIndex,1}) where {T}
    ptr = pointer(A) + (I[1].start-1)*sizeof(T)
    len = I[1].stop - I[1].start + 1
    return MtlDeviceArray(len, ptr)
end

@inline function Base.iterate(A::MtlDeviceArray, i=1)
    if (i % UInt) - 1 < length(A)
        (@inbounds A[i], i + 1)
    else
        nothing
    end
end


text:
using Revise
module QFactors
using Revise

using  CSV, DataFrames, Serialization, Dates, Statistics,
  Distributed, Finometrics, Measures,
  Formatting, CodecLz4, CodecZlib

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

############################CONSTANTS

const OUT_PATH = pwd() * "\\output"
const DATA_PATH = pwd() * "\\data"

const CRSP_NAME = "CRSP-M"
const CRSP_DATE_FORMAT = dateformat"yyyymmdd"

const COMP_NAME = "COMP"
const COMP_A_NAME = "COMP-A"
const COMP_Q_NAME = "COMP-Q"
const COMP_DATE_FORMAT = dateformat"yyyymmdd"

const CCM_NAME = "CCM"

###########################IO Compression/Decompression routines
IN_CSV_STREAM(p::String) = LZ4DecompressorStream(open(p))
IN_CSV_STREAM(F::Function, p::String) = open(F, LZ4DecompressorStream, p)

IN_JLS_STREAM(p::String) = LZ4DecompressorStream(open(p))
IN_JLS_STREAM(F::Function, p::String) = open(F, LZ4DecompressorStream, p)

OUT_JLS_STREAM(p::String) = LZ4CompressorStream(open(p, "w"))
OUT_JLS_STREAM(F::Function, p::String) = open(F, LZ4CompressorStream, p, "w")

##############################EXPORTS

############################FILES
include("QUtilities.jl")
include("PreprocessCRSP.jl")
include("PreprocessCOMP.jl")
include("MergeCOMPCRSP.jl")

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
### Implementation via recursion ###
# recursion approach
function decode(m::Tacotron₂, values::T, keys::T, stepsleft::Integer) where T <: DenseArray{<:Real,3}
   #initialize dimensions
   _, time_in, batchsize = size(values)
   querydim = size(m.attention.dense.W, 2)
   nmelfeatures = length(m.frameproj.b)
   # initialize parameters
   query   = gpu(zeros(Float32, querydim, batchsize))
   weights = gpu(zeros(Float32, time_in, 1, batchsize))
   frame   = gpu(zeros(Float32, nmelfeatures, batchsize))

   context, weights = m.attention(values, keys, query, weights, weights)
   prenetoutput = m.prenet(frame)
   prenetoutput_context = [prenetoutput; context]
   query = m.lstms(prenetoutput_context)
   query_context = [query; context]
   frame = m.frameproj(query_context)
   decoding₁ = (frame, query_context)
   decode(m, values, keys, stepsleft-1, query, weights, weights, frame, decoding₁)
end

function decode(m::Tacotron₂, values::T, keys::T, stepsleft::Integer, query::M, weights::T, Σweights::T, frame::M, decodings) where {M <: DenseMatrix, T <: DenseArray{<:Real,3}}
   (stepsleft == 0) && (return decodings)
   context, weights = m.attention(values, keys, query, weights, Σweights)
   Σweights += weights
   prenetoutput = m.prenet(frame)
   prenetoutput_context = [prenetoutput; context]
   query = m.lstms(prenetoutput_context)
   query_context = [query; context]
   frame = m.frameproj(query_context)
   decoding = (frame, query_context)
   decodings = [decodings; decoding]
   decode(m, values, keys, stepsleft-1, query, weights, Σweights, frame, decodings)
end

###
function (m::Tacotron₂)(textindices::DenseMatrix{<:Integer}, time_out::Integer)
   # dimensions
   batchsize = size(textindices, 2)
   nmelfeatures = length(m.frameproj.b)
   # encoding stage
   chex = m.che(textindices)
   convblock₃x = m.convblock₃(chex)
   values = m.blstm(convblock₃x)
   # @ein keys[a,t,b] := m.attention.V[a,d] * values[d,t,b] # dispatches to batched_contract (vetted)
   keys = einsum(EinCode{((1,2), (2,3,4)), (1,3,4)}(), (m.attention.V, values))
   # #=check=# Vh ≈ reduce(hcat, [reshape(m.V * values[:,t,:], size(m.V,1), 1, :) for t ∈ axes(values,2)])
   decodings = decode(m, values, keys, time_out)
   frames = first.(decodings)
   query_contexts = last.(decodings)
   σ⁻¹stoprobs = reshape(m.stopproj(reduce(hcat, query_contexts)), (batchsize, time_out))
   # #=check=# σ⁻¹stoprobs == reduce(hcat, reshape.(m.stopproj.(query_contexts), Val(1)))
   prediction = reshape(reduce(hcat, frames), (nmelfeatures, batchsize, time_out))
   # #=check=# prediction == cat(frames...; dims=3)
   melprediction = permutedims(prediction, (3,1,2))
   melprediction⁺residual = melprediction + m.postnet(melprediction)
   return melprediction, melprediction⁺residual, σ⁻¹stoprobs
end

### Implementation via iterator ###
# Warning: adjoint of the `collect(x::Decodings)` defined below is likely wrong
struct Decodings{T, M}
   m :: M
   values   :: T
   keys     :: T
   time_out :: Int
end

@adjoint function Decodings(m::Tacotron2, values::T, keys::T, time_out::Int) where T <: DenseArray{<:Real,3}
   return Decodings(m, values, keys, time_out), function (d̄)
      return d̄.m, d̄.values, d̄.keys, nothing
   end
end

# this adjoint is seems wrong
@adjoint function collect(x::Decodings)
   collect(x), function(ȳ)
      return (Decodings(first(ȳ)[3:end]..., length(ȳ)),)
   end
end

Flux.trainable(d::Decodings) = (d.m,)
@functor Decodings

Base.eltype(::Type{<:Decodings{T,M}}) where {T <: DenseArray{<:Real,3}, M <: Tacotron2} = Tuple{matrixof(T), vectorof(T), M, T, T}
Base.length(itr::Decodings) = itr.time_out

function Base.iterate(itr::Decodings)
   (itr.time_out <= 0) && (return nothing)
   #initialize dimensions
   _, time_in, batchsize = size(itr.values)
   querydim = size(m.attention.dense.W, 2)
   nmelfeatures = length(m.frameproj.b)
   # initialize parameters
   query   = gpu(zeros(Float32, querydim, batchsize))
   weights = gpu(zeros(Float32, time_in, 1, batchsize))
   frame   = gpu(zeros(Float32, nmelfeatures, batchsize))

   context, weights = itr.m.attention(itr.values, itr.keys, query, weights, weights)
   prenetoutput = itr.m.prenet(frame)
   prenetoutput_context = [prenetoutput; context]
   query = itr.m.lstms(prenetoutput_context)
   query_context = [query; context]
   frame = itr.m.frameproj(query_context)
   σ⁻¹pstopᵀ = itr.m.stopproj(query_context)
   σ⁻¹pstop = reshape(σ⁻¹pstopᵀ, Val(1))
   i = (frame, σ⁻¹pstop, m, values, keys)
   state = (1, query, weights, weights, frame)
   return i, state
end

function Base.iterate(itr::Decodings, state::Tuple{Int,M,T,T,M}) where {M <: DenseMatrix, T <: DenseArray{<:Real,3}}
   (t, query, weights, Σweights, frame) = state
   (t == itr.time_out) && (return nothing)
   t += 1
   context, weights = itr.m.attention(itr.values, itr.keys, query, weights, Σweights)
   Σweights += weights
   prenetoutput = itr.m.prenet(frame)
   prenetoutput_context = [prenetoutput; context]
   query = itr.m.lstms(prenetoutput_context)
   query_context = [query; context]
   frame = itr.m.frameproj(query_context)
   σ⁻¹pstopᵀ = itr.m.stopproj(query_context)
   σ⁻¹pstop = reshape(σ⁻¹pstopᵀ, Val(1))
   i = (frame, σ⁻¹pstop, m, values, keys)
   state = (t, query, weights, Σweights, frame)
   return i, state
end

function (m::Tacotron2)(textindices::DenseMatrix{<:Integer}, time_out::Integer)
   # dimensions
   batchsize = size(textindices, 2)
   nmelfeatures = length(m.frameproj.b)
   # encoding stage
   chex = m.che(textindices)
   convblock₃x = m.convblock₃(chex)
   values = m.blstm(convblock₃x)
   # @ein keys[a,t,b] := m.attention.V[a,d] * values[d,t,b] # dispatches to batched_contract (vetted)
   keys = einsum(EinCode{((1,2), (2,3,4)), (1,3,4)}(), (m.attention.V, values))
   # #=check=# Vh ≈ reduce(hcat, [reshape(m.V * values[:,t,:], size(m.V,1), 1, :) for t ∈ axes(values,2)])
   # decodings′ = Zygote.bufferfrom(Vector{Tuple{Matrix{Float32},Vector{Float32}}}(undef, time_out))
   # decodings′ .= Decodings(m, values, keys, time_out)
   # decodings = copy(decodings′)
   # decodings = [decoding for decoding ∈ decodings′]
   decodings = collect(Decodings(m, values, keys, time_out))
   frames = first.(decodings)
   σ⁻¹pstops = getindex.(decodings, 2)
   σ⁻¹stoprobs = reduce(hcat, σ⁻¹pstops)
   prediction = reshape(reduce(hcat, frames), nmelfeatures, batchsize, time_out)
   # #=check=# prediction == cat(frames...; dims=3)
   melprediction = permutedims(prediction, (3,1,2))
   melprediction⁺residual = melprediction + m.postnet(melprediction)
   return melprediction, melprediction⁺residual, σ⁻¹stoprobs
end


text:
export MomentConstraint

"""
Description:

    Constraints the marginal of the connected variable to an
    expectation ∫q(x)g(x)dx = G. The parameter η in the node function
    is actively adapted s.t. the marginal respects the above constraint.
    Implementation according to (van de Laar et al. "Chance-Constrained
    Active Inference", MIT Neural Computation, 2021).

    f(out) = exp(η g(out))

Interfaces:

    1. out

Construction:

    MomentConstraint(out; g=g, G=G, id=:my_node)
"""
mutable struct MomentConstraint <: SoftFactor
    id::Symbol
    interfaces::Array{Interface,1}
    i::Dict{Symbol, Interface}

    g::Function # Constraint function
    G::Float64 # Constraint value
    eta_init::Float64 # Initial value for eta optimization

    function MomentConstraint(out; g::Function, G::Float64, eta_init::Float64, id=ForneyLab.generateId(MomentConstraint))
        @ensureVariables(out)
        self = new(id, Vector{Interface}(undef, 1), Dict{Symbol,Interface}(), g, G, eta_init)
        ForneyLab.addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)

        return self
    end
end

slug(::Type{MomentConstraint}) = "E"

# A breaker message is required if interface is partnered with a moment constraint
requiresBreaker(interface::Interface, partner_interface::Interface, partner_node::MomentConstraint) = true

breakerParameters(interface::Interface, partner_interface::Interface, partner_node::MomentConstraint) = (Message{GaussianMeanVariance, Univariate}, ()) # Univariate only

# Constraints do not contribute to average energy
averageEnergy(::Type{MomentConstraint}, marg_out::ProbabilityDistribution) = 0.0


text:
importall Bukdu
importall Bukdu.Octo
importall Bukdu.Tag

type User
    name::String
    age::Int
    description::String
end

using Base.Test

form = change(User)
@test """<label for="user_name" />""" == label(form, :name)

f = change(User, name="foo bar")
@test """<label for="user_name" />""" == label(f, :name)
@test """<label for="user_name">Name</label>""" == label(f, :name, "Name")

@test """<input id="user_name" name="user[name]" type="text" value="foo bar" />""" == text_input(f, :name)

@test """
<select id="user_age" name="user[age]">
    <option value="18">18</option>
    <option value="19">19</option>
    <option value="20">20</option>
</select>""" == select(f, :age, 18:20)

user = User("tom", 20, "")
f = change(user)
@test """
<select id="user_age" name="user[age]">
    <option value="18">18</option>
    <option value="19">19</option>
    <option value="20" selected>20</option>
</select>""" == select(f, :age, 18:20)

@test """<textarea id="user_description" name="user[description]">
</textarea>""" == textarea(f, :description)


text:
using Test
using MPI
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.MPIStateArrays
using KernelAbstractions

device(A) = typeof(A) <: Array ? CPU() : CUDA()

function main()
    FT = Float64
    Nhorz = 3
    Nstack = 5
    N = 4
    DA = Array

    MPI.Initialized() || MPI.Init()

    comm = MPI.COMM_WORLD
    crank = MPI.Comm_rank(comm)
    csize = MPI.Comm_size(comm)

    Rrange = FT.(accumulate(+, 1:(Nstack + 1)))
    topology = StackedCubedSphereTopology(
        MPI.COMM_SELF,
        Nhorz,
        Rrange;
        boundary = (1, 2),
    )
    grid = DiscontinuousSpectralElementGrid(
        topology;
        FloatType = FT,
        DeviceArray = DA,
        polynomialorder = N,
        meshwarp = Topologies.cubedshellwarp,
    )

    let
        Np = (N + 1)^3
        activedofs = Array(grid.activedofs)
        vmaprecv = Array(grid.vmaprecv)
        active = union(1:(length(topology.realelems) * Np), vmaprecv)
        inactive = setdiff(1:(length(topology.elems) * Np), active)

        @test all(activedofs[active] .== true)
        @test all(activedofs[inactive] .== false)
    end

    #=
    @show elems       = topology.elems
    @show realelems   = topology.realelems
    @show ghostelems  = topology.ghostelems
    @show sendelems   = topology.sendelems
    @show elemtocoord = topology.elemtocoord
    @show elemtoelem  = topology.elemtoelem
    @show elemtoface  = topology.elemtoface
    @show elemtoordr  = topology.elemtoordr
    @show elemtobndy  = topology.elemtobndy
    @show nabrtorank  = topology.nabrtorank
    @show nabrtorecv  = topology.nabrtorecv
    @show nabrtosend  = topology.nabrtosend
    =#

    # Check x1x2x3 matches before comm
    x1 = @view grid.vgeo[:, Grids._x1, :]
    x2 = @view grid.vgeo[:, Grids._x2, :]
    x3 = @view grid.vgeo[:, Grids._x3, :]

    @test x1[grid.vmap⁻] ≈ x1[grid.vmap⁺]
    @test x2[grid.vmap⁻] ≈ x2[grid.vmap⁺]
    @test x3[grid.vmap⁻] ≈ x3[grid.vmap⁺]

    Np = (N + 1)^3
    x1x2x3 = MPIStateArray{FT}(
        topology.mpicomm,
        DA,
        Np,
        3,
        length(topology.elems),
        realelems = topology.realelems,
        ghostelems = topology.ghostelems,
        vmaprecv = grid.vmaprecv,
        vmapsend = grid.vmapsend,
        nabrtorank = topology.nabrtorank,
        nabrtovmaprecv = grid.nabrtovmaprecv,
        nabrtovmapsend = grid.nabrtovmapsend,
    )
    x1x2x3.data[:, :, topology.realelems] .= @view grid.vgeo[
        :,
        [Grids._x1, Grids._x2, Grids._x3],
        topology.realelems,
    ]

    event = Event(device(x1x2x3.data))
    event = MPIStateArrays.begin_ghost_exchange!(x1x2x3, dependencies = event)
    event = MPIStateArrays.end_ghost_exchange!(x1x2x3, dependencies = event)
    wait(device(x1x2x3.data), event)

    # Check x1x2x3 matches after
    x1 = @view x1x2x3.data[:, 1, :]
    x2 = @view x1x2x3.data[:, 2, :]
    x3 = @view x1x2x3.data[:, 3, :]

    @test x1[grid.vmap⁻] ≈ x1[grid.vmap⁺]
    @test x2[grid.vmap⁻] ≈ x2[grid.vmap⁺]
    @test x3[grid.vmap⁻] ≈ x3[grid.vmap⁺]

    nothing
end
isinteractive() || main()


text:
export BlackHole,
    compute_eddington_luminosity,
    compute_bolometric_luminosity,
    compute_mass_accretion_rate

struct BlackHole{T<:AbstractFloat}
    M::T
    mdot::T
    spin::T
    isco::T
    efficiency::T
    Rg::T
    function BlackHole(M, mdot, spin)
        Rg = G * M / C^2
        isco = compute_isco(spin)
        eff = compute_efficiency(spin)
        new{typeof(M)}(M, mdot, spin, isco, eff, Rg)
    end
end

function BlackHole(parameters::Parameters)
    M = parameters.M * M_SUN
    return BlackHole(M, parameters.mdot, parameters.spin)
end

"""
Computes the Eddington luminosity for the given black hole mass.
"""
function compute_eddington_luminosity(bh_mass; mu_electron = 1.17)
    return 4π * G * bh_mass * C * mu_electron / SIGMA_E
end
compute_eddington_luminosity(bh::BlackHole; mu_electron = 1.17) =
    compute_eddington_luminosity(bh.M, mu_electron = mu_electron)

"""
Computes the bolometric luminosity of the AGN given the black
hole mass and the normalized eddington rate.
"""
compute_bolometric_luminosity(bh_mass, eddington_rate; mu_electron = 1.17) =
    eddington_rate * compute_eddington_luminosity(bh_mass; mu_electron = mu_electron)
compute_bolometric_luminosity(bh::BlackHole; mu_electron = 1.17) =
    compute_bolometric_luminosity(bh.M, bh.mdot, mu_electron = 1.17)

"""
Computes the Innermost Stable Circular Orbit (ISCO) in units of Rg.

# Parameters
- spin : normalized spin parameter ∈ (-1, 1)
"""
function compute_isco(spin)
    z1 = 1 + (1 - spin^2)^(1 / 3) * ((1 + abs(spin))^(1 / 3) + (1 - abs(spin))^(1 / 3))
    z2 = sqrt(3 * spin^2 + z1^2)
    rms = 3 + z2 - sign(spin) * sqrt((3 - z1) * (3 + z1 + 2 * z2))
    return rms
end

"""
Computes the black hole accretion efficiency parameter from the isco.

# Parameters
- spin : black hole spin parameter ∈ (-1, 1)
"""
function compute_efficiency(spin)
    isco = compute_isco(spin)
    return 1 - sqrt(1 - 2 / (3 * isco))
end

"""
Computes the mass accretion rate to the black hole,
assuming no wind losses.

# Parameters
- M : black hole mass [in grams]
- mdot : normalized eddington rate ∈ [0,1]
- spin : black hole spin parameter ∈ (-1, 1)
"""
function compute_mass_accretion_rate(M, mdot, spin; mu_electron = 1.17)
    efficiency = compute_efficiency(spin)
    return compute_bolometric_luminosity(M, mdot, mu_electron = mu_electron) /
           (efficiency * C^2)
end

function compute_mass_accretion_rate(bh; mu_electron = 1.17)
    return compute_mass_accretion_rate(bh.M, bh.mdot, bh.spin, mu_electron = mu_electron)
end


text:
###################################################
#   Test inicialización de pesos
###################################################
using Test
using Random
Random.seed!(2);

include("single-input-single-output.test.jl")
include("multiple-input-single-output.test.jl")
include("multiple-input-multiple-output.test.jl")


text:
"""
smoother of solution

ref:
- Dechter, R., Meiri, I., and Pearl, J. (1991).
  Temporal constraint networks.
  Artificial Intelligence (AIJ).
- Honing, W. et al. (2016).
  Multi-agent path finding with kinematic constraints.
  In ICAPS.
"""

mutable struct Action{State<:AbsState}
    id::String
    from::Node{State}
    to::Node{State}
    t::Int64
    agent::Int64
    predecessors::Vector{Tuple{Int64,String}}
    successors::Vector{Tuple{Int64,String}}
end

function get_action_id(
    from::Node{State},
    to::Node{State},
    t::Int64,
)::String where {State<:AbsState}
    return string(from.id) * "_" * string(to.id) * "_" * string(t)
end

"""create temporal plan graph temporal plan graph"""
function get_temporal_plan_graph(
    solution::Vector{Vector{Node{State}}},
    collide::Function,
    connect::Function;
    skip_connection::Bool = true,
)::Vector{Vector{Action{State}}} where {State<:AbsState}

    N = length(solution[1])
    T = length(solution)

    # temporal plan graph
    TPG = [Vector{Action{State}}() for i = 1:N]

    # type-1 dependency
    for i = 1:N
        v_last = solution[1][i]
        for (t, Q) in enumerate(solution[1:end])
            v_current = Q[i]

            if v_current.id != v_last.id
                action = Action(
                    get_action_id(v_last, v_current, t),
                    v_last,
                    v_current,
                    t,
                    i,
                    Vector{Tuple{Int64,String}}(),
                    Vector{Tuple{Int64,String}}(),
                )

                # add type-1 dependency
                if length(TPG[i]) > 0
                    action_pre = TPG[i][end]
                    push!(action.predecessors, (i, action_pre.id))
                end

                # insert
                push!(TPG[i], action)

                # update
                v_last = v_current
            end
        end
    end

    # type-2 dependency
    for i = 1:N
        for action_self in TPG[i]
            for j = 1:N
                if j == i
                    continue
                end

                # exclude ealier actions
                for action_other in filter(a -> a.t > action_self.t, TPG[j])

                    # create type-2 dependency
                    if collide(
                        action_self.from.q,
                        action_self.to.q,
                        action_other.from.q,
                        action_other.to.q,
                        i,
                        j,
                    )
                        push!(
                            action_other.predecessors,
                            (action_self.agent, action_self.id),
                        )
                        break
                    end
                end
            end
        end
    end

    # remove redundant dependencies
    for i = 1:N
        for action in TPG[i]
            latest_actions = Dict{Int64,Action}()
            for (j, id) in action.predecessors
                action_other = TPG[j][findfirst(a -> a.id == id, TPG[j])]
                if !haskey(latest_actions, j) || latest_actions[j].t < action_other.t
                    latest_actions[j] = action_other
                end
            end
            action.predecessors = map(a -> (a.agent, a.id), values(latest_actions))

            # update successors
            for (j, id) in action.predecessors
                action_other = TPG[j][findfirst(a -> a.id == id, TPG[j])]
                push!(action_other.successors, (i, action.id))
            end
        end
    end

    skip_connection && try_skip_connection!(TPG, collide, connect)
    return TPG
end

"""remove redundant actions in temporal plan graph"""
function try_skip_connection!(
    TPG::Vector{Vector{Action{State}}},
    collide::Function,
    connect::Function,
) where {State<:AbsState}

    N = length(TPG)
    get_action = (j, id) -> TPG[j][findfirst(a -> a.id == id, TPG[j])]

    for i = 1:N
        k = 1
        while k < length(TPG[i])
            k += 1
            a1 = TPG[i][k-1]
            a2 = TPG[i][k]

            # check connection
            !connect(a1.from.q, a2.to.q, i) && continue

            # check dependencies
            any([j != i for (j, _) in a2.predecessors]) && continue
            any([j != i for (j, _) in a1.successors]) && continue

            # new action candidate
            a3 = Action(
                get_action_id(a1.from, a2.to, a1.t),
                a1.from,
                a2.to,
                a1.t,
                a1.agent,
                a1.predecessors,
                a2.successors,
            )

            # get causalities
            causal_actions = get_ancestors(a1, TPG)
            union!(causal_actions, get_descendants(a2, TPG))
            causal_actions = filter(val -> val[1] != i, causal_actions)

            # check collisions
            conflicted = false
            for j = 1:N
                if j == i
                    continue
                end
                for a4 in filter(a -> !((j, a.id) in causal_actions), TPG[j])
                    if collide(a3.from.q, a3.to.q, a4.from.q, a4.to.q, i, j)
                        conflicted = true
                        break
                    end
                end
                if conflicted
                    continue
                end
                # check last location
                if collide(a3.from.q, a3.to.q, TPG[j][end].to.q, TPG[j][end].to.q, i, j)
                    conflicted = true
                    break
                end
            end
            if conflicted
                continue
            end

            # update action orders
            deleteat!(TPG[i], k-1:k)
            insert!(TPG[i], k - 1, a3)

            # update dependencies
            for (j, id) in a3.predecessors
                a4 = get_action(j, id)
                m = findfirst(tpl -> tpl[1] == i && tpl[2] == a1.id, a4.successors)
                deleteat!(a4.successors, m)
                push!(a4.successors, (i, a3.id))
            end
            for (j, id) in a3.successors
                a4 = get_action(j, id)
                m = findfirst(tpl -> tpl[1] == i && tpl[2] == a2.id, a4.predecessors)
                deleteat!(a4.predecessors, m)
                push!(a4.predecessors, (i, a3.id))
            end
            k -= 1
        end
    end
end

"""obtain actions with causal dependencies"""
function get_causal_actions(
    action::Action{State},
    TPG::Vector{Vector{Action{State}}},
    for_ancestors::Bool,
) where {State<:AbsState}

    N = length(TPG)
    tables = [Dict{String,Set{Tuple{Int64,String}}}() for i = 1:N]
    function f(i, id)
        action = TPG[i][findfirst(a -> a.id == id, TPG[i])]
        if haskey(tables[i], id)
            return tables[i][id]
        end
        causal_actions = Set{Tuple{Int64,String}}()
        for (j, id_j) in ((for_ancestors) ? action.predecessors : action.successors)
            push!(causal_actions, (j, id_j))
            foreach(val -> push!(causal_actions, val), f(j, id_j))
        end
        tables[i][id] = causal_actions
        return causal_actions
    end
    return f(action.agent, action.id)
end

function get_ancestors(
    action::Action{State},
    TPG::Vector{Vector{Action{State}}},
) where {State<:AbsState}
    return get_causal_actions(action, TPG, true)
end

function get_descendants(
    action::Action{State},
    TPG::Vector{Vector{Action{State}}},
) where {State<:AbsState}
    return get_causal_actions(action, TPG, false)
end

"""sampling one solution from temporal plan graph"""
function get_greedy_solution(
    TPG::Vector{Vector{Action{State}}},
    config_goal::Vector{State},
)::Vector{Vector{Node{State}}} where {State<:AbsState}

    N = length(TPG)
    indexes = [1 for i = 1:N]  # internal clock of agents
    solution = [[
        (length(arr) > 0 ? arr[1].from : Node(config_goal[i], 0, Vector{Int64}())) for
        (i, arr) in enumerate(TPG)
    ]]

    while !all([k > length(TPG[i]) for (i, k) in enumerate(indexes)])
        config = []
        indexes_last = copy(indexes)
        for i = 1:N
            if length(TPG[i]) > 0
                action = TPG[i][min(indexes[i], length(TPG[i]))]
                if all([
                    findfirst(a -> a.id == id, TPG[j]) < indexes_last[j] for
                    (j, id) in action.predecessors
                ])
                    # increment
                    indexes[i] = min(indexes[i] + 1, length(TPG[i]) + 1)
                    push!(config, action.to)
                else
                    # stay
                    push!(config, action.from)
                end
            else
                push!(config, Node(config_goal[i], 0, Vector{Int64}()))
            end
        end
        push!(solution, config)
    end
    return solution
end

"""
    get_solution_cost(
        solution::Union{Nothing,Vector{Vector{Node{State}}}},
    )::Union{Nothing,Dict{Symbol,Float64}} where {State<:AbsState}

compute costs (sum_of_costs & makespan) of solution
"""
function get_solution_cost(
    solution::Union{Nothing,Vector{Vector{Node{State}}}},
)::Union{Nothing,Dict{Symbol,Float64}} where {State<:AbsState}

    isnothing(solution) && return nothing
    N = length(solution[1])
    T = length(solution)

    # compute last timesteps
    last_timesteps = fill(T, N)
    for i = 1:N
        for t in reverse(collect(1:T-1))
            solution[t][i] != solution[t+1][i] && break
            last_timesteps[i] = t
        end
    end

    makespan = 0
    sum_of_cost = 0
    for t = 2:T
        c = maximum(i -> dist(solution[t-1][i], solution[t][i]), 1:N)
        sum_of_cost += c * length(filter(e -> t <= e, last_timesteps))
        makespan += c
    end

    return Dict(:sum_of_cost => sum_of_cost, :makespan => makespan)
end

"""
    get_tpg_cost(
        TPG::Union{Nothing,Vector{Vector{Action{State}}}},
    )::Union{Nothing,Dict{Symbol,Float64}} where {State<:AbsState}

compute costs (sum_of_costs & makespan) of temporal plan graph (by dynamic programming)
"""
function get_tpg_cost(
    TPG::Union{Nothing,Vector{Vector{Action{State}}}},
)::Union{Nothing,Dict{Symbol,Float64}} where {State<:AbsState}

    isnothing(TPG) && return nothing

    N = length(TPG)
    cost_tables = [Dict{String,Float64}() for i = 1:N]

    function f(i, id)
        action = TPG[i][findfirst(a -> a.id == id, TPG[i])]
        c = dist(action.from.q, action.to.q)
        if !isempty(action.predecessors)
            c += maximum([
                (haskey(cost_tables[j], id_j) ? cost_tables[j][id_j] : f(j, id_j)) for
                (j, id_j) in action.predecessors
            ])
        end
        cost_tables[i][id] = c
        return c
    end

    arr = [(length(TPG[i]) > 0 ? f(i, TPG[i][end].id) : 0) for i = 1:N]
    return Dict(:sum_of_cost => sum(arr), :makespan => maximum(arr))
end

"""
    smoothing(
        solution::Vector{Vector{Node{State}}},
        connect::Function,
        collide::Function;
        VERBOSE::Int64 = 0,
    )::Tuple{
        Vector{Vector{Action{State}}},  # temporal plan graph
        Vector{Vector{Node{State}}},    # solution
        Dict{Symbol,Float64},           # final cost
    } where {State<:AbsState}

smoothing by temporal plan graph
"""
function smoothing(
    solution::Vector{Vector{Node{State}}},
    connect::Function,
    collide::Function;
    VERBOSE::Int64 = 0,
)::Tuple{
    Vector{Vector{Action{State}}},  # temporal plan graph
    Vector{Vector{Node{State}}},  # solution
    Dict{Symbol,Float64},  # final cost
} where {State<:AbsState}

    isnothing(solution) && return nothing

    solution_tmp = solution
    config_goal = map(v -> v.q, solution[end])
    cost_last = nothing
    sum_of_cost_last = Inf
    while true
        # 1. create temporal plan graph
        TPG = get_temporal_plan_graph(solution_tmp, collide, connect)
        # 2. sampling from temporal plan graph
        solution_tmp = get_greedy_solution(TPG, config_goal)
        cost = get_tpg_cost(TPG)
        if sum_of_cost_last >= cost[:sum_of_cost]
            return (TPG, solution_tmp, cost)
        else
            # 3. update solution
            VERBOSE > 0 && @info @sprintf("cost is updated: %f -> %f", cost_last, cost)
            cost_last = cost
            sum_of_cost_last = get(cost, :sum_of_cost)
        end
    end
end

function smoothing(solution::Nothing, args...; kwargs...)::Nothing
    nothing
end


text:
using Documenter
using LearnConvection

makedocs(
    sitename = "LearnConvection",
    format = Documenter.HTML(),
    modules = [LearnConvection],
    pages = [
        "Home" => "index.md",
        "Library" => "main_function_index.md"
        # "Library" => Any[
        #     "Main" => "main_function_index.md",
        #     "GaussianProcess module" => "gaussian_process_function_index.md",
        #     "Data module" => "data_function_index.md"
        # ]
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/adelinehillier/LearnConvection.git"
)

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
using Distributions
using SampleDistributions
using BenchmarkTools

d = SampleDistribution(['a', 'a', 'b', 'c'])
d2 = SampleDistribution([1, 1, 2, 3])

v = rand(Char['a', 'b', 'c', 'd'], 1000)
D = SampleDistribution(v)


text:
#!/usr/bin/env julia

using ZMQ

println("Current ZMQ version is $(ZMQ.version)")

max_stars_repo_name: JailbreakFox/LightWeightRepository

max_issues_repo_name: JailbreakFox/LightWeightRepository

max_forks_repo_name: JailbreakFox/LightWeightRepository


text:
##### specific distributions #####

const discrete_distributions = [
    # "bernoulli",
    # "betabinomial",
    # "binomial",
    # "dirac",
    # "discreteuniform",
    # "discretenonparametric",
    # "categorical",
    # "geometric",
    # "hypergeometric",
    # "negativebinomial",
    # "noncentralhypergeometric",
    # "poisson",
    # "skellam",
    # "soliton",
    # "poissonbinomial"
]

const continuous_distributions = [
    # "arcsine",
    # "beta",
    # "betaprime",
    # "biweight",
    # "cauchy",
    # "chernoff",
    # "chisq",    # Chi depends on Chisq
    # "chi",
    # "cosine",
    # "epanechnikov",
     "exponential",
    # "fdist",
    # "frechet",
    "gamma",
    # "erlang",
    # "pgeneralizedgaussian", # GeneralizedGaussian depends on Gamma
    # "generalizedpareto",
    # "generalizedextremevalue",
    # "gumbel",
    # "inversegamma",
    # "inversegaussian",
    # "kolmogorov",
    # "ksdist",
    # "ksonesided",
    # "laplace",
    # "levy",
    # "locationscale",
    # "logistic",
    # "noncentralbeta",
    # "noncentralchisq",
    # "noncentralf",
    # "noncentralt",
     "normal",
    # "normalcanon",
    # "normalinversegaussian",
     "lognormal",    # LogNormal depends on Normal
     "logitnormal",    # LogitNormal depends on Normal
    # "pareto",
    # "rayleigh",
    # "semicircle",
    # "skewnormal",
    # "studentizedrange",
    # "symtriangular",
    # "tdist",
    # "triangular",
    # "triweight",
    # "uniform",
    # "vonmises",
    "weibull"
]

for dname in discrete_distributions
    include(joinpath("univariate", "discrete", "$(dname).jl"))
end

include("univariate/continuous/estimateMoments.jl")
for dname in continuous_distributions
    include(joinpath("univariate", "continuous", "$(dname).jl"))
end


text:
using YAML
using Streamfall

network = YAML.load_file("../test/data/campaspe/campaspe_network.yml")
sn = create_network("Example Network", network)

inlets, outlets = find_inlets_and_outlets(sn)

@info "Network has the following inlets and outlets:" inlets outlets


text:
# tests for the various forms of stParallel for the ScatteringTransform
# exit()
# using Revise
using Distributed
addprocs(min((Sys.CPU_THREADS) - 2 - nprocs(), 2))
@everywhere using Interpolations, ContinuousWavelets, Wavelets
@everywhere using FFTW
@everywhere using ParallelScattering
@everywhere using SharedArrays
using Test
using AbstractFFTs, LinearAlgebra, Statistics, Random

include("planTests.jl") # all from Parallel
include("cwtTests.jl") # also from Parallel


text:
export Surgery

module Surgery

using CSV
using ..Empirikos: BinomialSample

const DATA = joinpath(@__DIR__, "surgery.csv")

function load_table()
    CSV.File(DATA)
end

function ebayes_samples()
    tbl = load_table()
    Zs = BinomialSample.(tbl.s, tbl.n)
end

end


text:
@testset "ladderize" begin
    start_tree = ParseNewick("((A,(B,(C,(D,E)))),(F,(G,H)));")
    descending_tree = ParseNewick("(((((D,E),C),B),A),((G,H),F));")
    ascending_tree = ParseNewick("((F,(G,H)),(A,(B,(C,(D,E)))));")

    start_newick = newick(start_tree)
    desc_newick = newick(descending_tree)
    asc_newick = newick(ascending_tree)

    ladderized_asc_tree = ladderize_tree(start_tree)
    @test newick(ladderized_asc_tree) == asc_newick

    ladderized_desc_tree = ladderize_tree(start_tree,false)
    @test newick(ladderized_desc_tree) == desc_newick

    ladderize_tree!(start_tree)
    @test newick(start_tree) == asc_newick

    start_tree = ParseNewick("((A,(B,(C,(D,E)))),(F,(G,H)));")
    ladderize_tree!(start_tree, false)
    @test newick(start_tree) == desc_newick
end


text:
using Cxx
using CxxStd: WrappedCppPrimArray

cxx"""
#include <tf/transform_datatypes.h>
"""

const tfVector3    = Union{cxxt"tf::Vector3",    cxxt"tf::Vector3&"}
const tfQuaternion = Union{cxxt"tf::Quaternion", cxxt"tf::Quaternion&"}
const tfMatrix3x3  = Union{cxxt"tf::Matrix3x3",  cxxt"tf::Matrix3x3&"}

function Base.unsafe_wrap(::Type{DenseArray}, v::T) where {T<:tfVector3}
    WrappedCppPrimArray(icxx"&$v[0];",3)
end

function Base.unsafe_wrap(::Type{DenseArray}, q::T) where {T<:tfQuaternion}
    WrappedCppPrimArray(icxx"&$q[0];",4)
end

function Base.convert(::Type{CT}, v::T) where {CT<:AbstractArray,T<:Union{tfVector3, tfQuaternion}}
    convert(CT, unsafe_wrap(DenseArray, v))
end

function Base.convert(::Type{CT}, m::T) where {CT<:AbstractArray, T<:tfMatrix3x3}
    [DenseArray(icxx"$m[0];"), DenseArray(icxx"$m[1];"), DenseArray(icxx"$m[2];")]
end

max_stars_repo_name: moyhig/CxxROS

max_issues_repo_name: moyhig/CxxROS

max_forks_repo_name: moyhig/CxxROS


text:
#=
DataIterator
=#

import Random.shuffle!

struct Dataset
    len::Int
    data
end

Base.length(dataset::Dataset) = dataset.len

struct SerialIterator
    len::Int
    get_example
    batchsize::Int
    indices::Vector
    function SerialIterator(dataset::Dataset, batchsize; shuffle=true)
        indices=ifelse(shuffle, Vector(1:dataset.len),shuffle!(Vector(1:dataset.len)))
        new(length(dataset), dataset.get_example, batchsize, indices)
    end
end

function Base.iterate(diter::SerialIterator, state=(1, 0))
    idx_start, count = state
    if idx_start + diter.batchsize > diter.len
        return nothing
    else
        indices = diter.indices[idx_start:idx_start+diter.batchsize]
        element = diter.get_example.(indices)
        return (element, (idx_start + diter.batchsize, count + 1))
    end
end

function main()
    len=10000000
    dataset = Dataset(len, i->rand(30,30))
    dataset_iter = SerialIterator(dataset, 32, shuffle=true)
    for i in dataset_iter
        println(i)
    end
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
function parse_multipart_image(::Val{:jpeg}, body::AbstractVector{UInt8})
    body_length = length(body)
    body_length_minus_one = body_length - 1
    image_starts = Vector{Int}(undef, 0)
    image_ends = Vector{Int}(undef, 0)
    for i = 1:body_length_minus_one
        body_i = body[i]
        body_i_plus_one = body[i+1]
        if body_i == 0xff && body_i_plus_one == 0xd8
            push!(image_starts, i)
        end
        if body_i == 0xff && body_i_plus_one == 0xd9
            push!(image_ends, i+1)
        end
    end
    num_image_starts = length(image_starts)
    num_image_ends = length(image_ends)
    if num_image_starts !== num_image_ends
        throw(ArgumentError("There are $(num_image_starts) starts but $(num_image_ends) ends"))
    end
    image_bytes = Vector{Vector{UInt8}}(undef, num_image_starts)
    for i = 1:num_image_starts
        image_start = image_starts[i]
        image_end = image_ends[i]
        bytes = body[image_start:image_end]
        image_bytes[i] = bytes
    end
    images = [_bytes_to_image("jpeg", bytes) for bytes in image_bytes]
    return images
end


text:
push!(LOAD_PATH, pwd()*"/../../src/")
using Distributions
using PyPlot
using BVB

x = linspace(-0.5,1.5,500)
max_iter = 200
init_approx_d = Normal(-1.0, 2)

# example 1
@show true_d = Exponential(1.0)
ll(x) = logpdf(true_d, x)

include("run_laplacian_gradient_boosting.jl")
show()

println("Done.")
readline(STDIN)

max_stars_repo_name: richardkwo/boosting-VB

max_issues_repo_name: richardkwo/boosting-VB

max_forks_repo_name: richardkwo/boosting-VB


text:
export basin_fractions, tipping_probabilities

"""
    basin_fractions(basins::Array) → fs::Dict
Calculate the fraction of the basins of attraction encoded in `basins`.
The elements of `basins` are integers, enumerating the attractor that the entry of `basins`
converges to. Return a dictionary that maps attractor IDs to their relative fractions.

In[^Menck2013] the authors use these fractions to quantify the stability of a basin of
attraction, and specifically how it changes when a parameter is changed.

[^Menck2013]: Menck, Heitzig, Marwan & Kurths. How basin stability complements the linear stability paradigm. [Nature Physics, 9(2), 89–92](https://doi.org/10.1038/nphys2516)
"""
function basin_fractions(basins::AbstractArray)
    fs = Dict{eltype(basins), Float64}()
    ids = unique(basins)
    N = length(basins)
    for ξ in ids
        B = count(isequal(ξ), basins)
        fs[ξ] = B/N
    end
    return fs
end

"""
    tipping_probabilities(basins_before, basins_after) → P
Return the tipping probabilities of the computed basins before and after a change
in the system parameters (or time forcing), according to the definition of[^Kaszás2019].

The input `basins` are integer-valued arrays, where the integers enumerate the attractor,
e.g. the output of [`basins_of_attraction`](@ref).

## Description
Let ``\\mathcal{B}_i(p)`` denote the basin of attraction of attractor ``A_i`` at
parameter(s) ``p``. Kaszás et al[^Kaszás2019] define the tipping probability
from ``A_i`` to ``A_j``, given a parameter change in the system of ``p_- \\to p_+``, as

```math
P(A_i \\to A_j | p_- \\to p_+) =
\\frac{|\\mathcal{B}_j(p_+) \\cap \\mathcal{B}_i(p_-)|}{|\\mathcal{B}_i(p_-)|}
```
where ``|\\cdot|`` is simply the volume of the enclosed set.
The value of `` P(A_i \\to A_j | p_- \\to p_+)`` is `P[i, j]`.
The equation describes something quite simple:
what is the overlap of the basin of attraction of ``A_i`` at ``p_-`` with that of the
attractor ``A_j`` at ``p_+``.
If `basins_before, basins_after` contain values of `-1`, corresponding to trajectories
that diverge, this is considered as the last attractor of the system in `P`.

[^Kaszás2019]: Kaszás, Feudel & Tél. Tipping phenomena in typical dynamical systems subjected to parameter drift. [Scientific Reports, 9(1)](https://doi.org/10.1038/s41598-019-44863-3)
"""
function tipping_probabilities(basins_before::AbstractArray, basins_after::AbstractArray)
    @assert size(basins_before) == size(basins_after)

    bid, aid = unique.((basins_before, basins_after))
    P = zeros(length(bid), length(aid))
    N = length(basins_before)
    # Make -1 last entry in bid, aid, if it exists
    put_minus_1_at_end!(bid); put_minus_1_at_end!(aid)

    # Notice: the following loops could be optimized with smarter boolean operations,
    # however they are so fast that everything should be done within milliseconds even
    # on a potato
    for (i, ι) in enumerate(bid)
        B_i = findall(isequal(ι), basins_before)
        μ_B_i = length(B_i) # μ = measure
        for (j, ξ) in enumerate(aid)
            B_j = findall(isequal(ξ), basins_after)
            μ_overlap = length(B_i ∩ B_j)
            P[i, j] = μ_overlap/μ_B_i
        end
    end
    return P
end

function put_minus_1_at_end!(bid)
    if -1 ∈ bid
        sort!(bid)
        popfirst!(bid)
        push!(bid, -1)
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 3b86934d-dbef-4ff5-b555-47d6a9456efa
begin
	using Futbolista: DataLoaders
	const sb = DataLoaders.Statsbomb
end

# ╔═╡ cd78a60c-56db-4f35-a82c-26e9b5335852
begin
	rldr = sb.SBLoader()
	lldr = sb.SBLoader("local", raw"C:\repository\open-data\data")
end

# ╔═╡ 3d945a95-e6ff-44a7-8f42-46064e2d76d1
begin
	comp_df = sb.load_competitions(lldr)
	matches_df = sb.load_matches(lldr, 72, 30)
	events_df = sb.load_events(lldr, 7298);
end

# ╔═╡ 6aab4a57-ba7c-4c9b-837f-e51164b991ee
matches_df

# ╔═╡ 6daad45f-e62c-4963-ab09-4fa307bd914d
comp_df

# ╔═╡ Cell order:
# ╠═3b86934d-dbef-4ff5-b555-47d6a9456efa
# ╠═cd78a60c-56db-4f35-a82c-26e9b5335852
# ╠═3d945a95-e6ff-44a7-8f42-46064e2d76d1
# ╠═6aab4a57-ba7c-4c9b-837f-e51164b991ee
# ╠═6daad45f-e62c-4963-ab09-4fa307bd914d


text:
export
ruleSPMultiplicationOutNPG,
ruleSPMultiplicationOutNGP,
ruleSPMultiplicationOutNPP,
ruleSPMultiplicationIn1GNP,
ruleSPMultiplicationIn1PNP,
ruleSPMultiplicationAGPN,
ruleSPMultiplicationAPPN,
ruleSPMultiplicationOutNPΓ,
ruleSPMultiplicationOutNΓP,
ruleSPMultiplicationIn1ΓNP,
ruleSPMultiplicationAΓPN

#-------------------------------
# Univariate (in1 and a commute)
#-------------------------------

function ruleSPMultiplicationOutNPG(msg_out::Nothing,
                                    msg_in1::Message{PointMass, Univariate},
                                    msg_a::Message{F, Univariate}) where F<:Gaussian

    d_a = convert(ProbabilityDistribution{Univariate, GaussianMeanVariance}, msg_a.dist)

    Message(Univariate, GaussianMeanVariance, m=msg_in1.dist.params[:m]*d_a.params[:m], v=d_a.params[:v]*msg_in1.dist.params[:m]^2)
end

ruleSPMultiplicationOutNGP(msg_out::Nothing, msg_in1::Message{F, Univariate}, msg_a::Message{PointMass, Univariate}) where F<:Gaussian =
    ruleSPMultiplicationOutNPG(nothing, msg_a, msg_in1)

function ruleSPMultiplicationIn1GNP(msg_out::Message{F, Univariate},
                                    msg_in1::Nothing,
                                    msg_a::Message{PointMass, Univariate}) where F<:Gaussian

    d_out = convert(ProbabilityDistribution{Univariate, GaussianWeightedMeanPrecision}, msg_out.dist)

    a = msg_a.dist.params[:m]

    Message(Univariate, GaussianWeightedMeanPrecision, xi=a*d_out.params[:xi], w=d_out.params[:w]*a^2)
end

ruleSPMultiplicationAGPN(msg_out::Message{F, Univariate}, msg_in1::Message{PointMass, Univariate}, msg_a::Nothing) where F<:Gaussian =
    ruleSPMultiplicationIn1GNP(msg_out, nothing, msg_in1)

ruleSPMultiplicationOutNPP(msg_out::Nothing, msg_in1::Message{PointMass, Univariate}, msg_a::Message{PointMass, Univariate}) =
    Message(Univariate, PointMass, m=msg_in1.dist.params[:m]*msg_a.dist.params[:m])

ruleSPMultiplicationIn1PNP(msg_out::Message{PointMass, Univariate}, msg_in1::Nothing, msg_a::Message{PointMass, Univariate}) =
    Message(Univariate, PointMass, m=msg_out.dist.params[:m]/msg_a.dist.params[:m])

ruleSPMultiplicationAPPN(msg_out::Message{PointMass, Univariate}, msg_in1::Message{PointMass, Univariate}, msg_a::Nothing) =
    Message(Univariate, PointMass, m=msg_out.dist.params[:m]/msg_in1.dist.params[:m])

function ruleSPMultiplicationOutNPΓ(msg_out::Nothing,
                                    msg_in1::Message{PointMass, Univariate},
                                    msg_a::Message{Gamma})

    Message(Univariate, Gamma, a=msg_a.dist.params[:a], b=msg_a.dist.params[:b]/msg_in1.dist.params[:m])
end

ruleSPMultiplicationOutNΓP(msg_out::Nothing, msg_in1::Message{Gamma}, msg_a::Message{PointMass, Univariate}) =
    ruleSPMultiplicationOutNPΓ(nothing, msg_a, msg_in1)

function ruleSPMultiplicationIn1ΓNP(msg_out::Message{Gamma},
                                    msg_in1::Nothing,
                                    msg_a::Message{PointMass, Univariate})

    Message(Univariate, Gamma, a=msg_out.dist.params[:a], b=msg_out.dist.params[:b]*msg_a.dist.params[:m])
end

ruleSPMultiplicationAΓPN(msg_out::Message{Gamma}, msg_in1::Message{PointMass, Univariate}, msg_a::Nothing) =
    ruleSPMultiplicationIn1ΓNP(msg_out, nothing, msg_in1)

#------------------------------------------------------
# MatrixVariate*multivariate (in1 and a do NOT commute)
#------------------------------------------------------

function ruleSPMultiplicationOutNGP(msg_out::Nothing,
                                    msg_in1::Message{F, Multivariate},
                                    msg_a::Message{PointMass, MatrixVariate}) where F<:Gaussian

    d_in1 = convert(ProbabilityDistribution{Multivariate, GaussianMeanVariance}, msg_in1.dist)

    A = msg_a.dist.params[:m]

    Message(Multivariate, GaussianMeanVariance, m=A*d_in1.params[:m], v=A*d_in1.params[:v]*A')
end

function ruleSPMultiplicationIn1GNP(msg_out::Message{F, Multivariate},
                                    msg_in1::Nothing,
                                    msg_a::Message{PointMass, MatrixVariate}) where F<:Gaussian

    d_out = convert(ProbabilityDistribution{Multivariate, GaussianWeightedMeanPrecision}, msg_out.dist)

    A = msg_a.dist.params[:m]
    W = A'*d_out.params[:w]*A
    W = W + tiny*diageye(size(W)[1]) # Ensure precision is always invertible

    Message(Multivariate, GaussianWeightedMeanPrecision, xi=A'*d_out.params[:xi], w=W)
end

ruleSPMultiplicationOutNPP(msg_out::Nothing, msg_in1::Message{PointMass, Multivariate}, msg_a::Message{PointMass, MatrixVariate}) =
    Message(Multivariate, PointMass, m=msg_a.dist.params[:m]*msg_in1.dist.params[:m])

ruleSPMultiplicationIn1PNP(msg_out::Message{PointMass, Multivariate}, msg_in1::Nothing, msg_a::Message{PointMass, MatrixVariate}) =
    Message(Multivariate, PointMass, m=pinv(msg_a.dist.params[:m])*msg_out.dist.params[:m])

#------------------------
# Univariate*multivariate
#------------------------

# We consider the following updates as a special case of the MatrixVariate*Multivariate updates.
# Namely, Ax = y, where A ∈ R^{nx1}, x ∈ R^1, and y ∈ R^n. In this case, the matrix A
# can be represented by a n-dimensional vector, and x by a scalar. Before computation,
# quantities are converted to their proper dimensions (see situational sketch below).
#
#     | a ~ Multivariate -> R^{nx1}
#     v  out ~ Multivariate -> R^n
# -->[x]-->
# in1 ~ Univariate -> R^1

function ruleSPMultiplicationOutNGP(msg_out::Nothing,
                                    msg_in1::Message{F, Univariate},
                                    msg_a::Message{PointMass, Multivariate}) where F<:Gaussian

    dist_in1_mult = convert(ProbabilityDistribution{Multivariate, F}, msg_in1.dist)
    dist_a_matr = convert(ProbabilityDistribution{MatrixVariate, PointMass}, msg_a.dist)

    return ruleSPMultiplicationOutNGP(nothing, Message(dist_in1_mult), Message(dist_a_matr))
end

function ruleSPMultiplicationIn1GNP(msg_out::Message{F, Multivariate},
                                    msg_in1::Nothing,
                                    msg_a::Message{PointMass, Multivariate}) where F<:Gaussian

    dist_a_matr = convert(ProbabilityDistribution{MatrixVariate, PointMass}, msg_a.dist)
    msg_in1_mult = ruleSPMultiplicationIn1GNP(msg_out, nothing, Message(dist_a_matr))
    (dims(msg_in1_mult.dist) == 1) || error("Implicit conversion to Univariate failed for $(msg_in1_mult.dist)")

    return Message(Univariate, GaussianWeightedMeanPrecision, xi=msg_in1_mult.dist.params[:xi][1], w=msg_in1_mult.dist.params[:w][1,1])
end


text:
function get_nbins(alg, data)

    n = length(data)

    if alg == :sqrt
        nbins = ceil(Int, sqrt(n))
    elseif alg == :sturges
        nbins = ceil(Int, log(2,n)) + 1
    elseif alg == :rice
        nbins = ceil(Int, 2*cbrt(n))
    elseif alg == :doane
        g₁ = moment(data, 3)
        σ = sqrt((6*(n-2))/((n+1)*(n+3)))
        nbins = ceil(Int, 1 + log(2,n) + log(2, 1+abs(g₁)/σ))
    elseif alg == :scott
        σ = std(data)
        binwidth = 3.5*σ/cbrt(n)
        lo, hi = extrema(data)
        nbins = ceil(Int, (hi - lo)/binwidth)
    elseif alg == :fd
        binwidth = 2*iqr(data)/cbrt(n)
        lo, hi = extrema(data)
        nbins = ceil((hi - lo)/binwidth)
    else
        binwidth = 2*iqr(data)/cbrt(n)
        lo, hi = extrema(data)
        nbins_fd = ceil(Int, (hi - lo)/binwidth)
        nbins_sturges = ceil(Int, log(2,n)) + 1
        nbins = max(nbins_fd, nbins_sturges)
    end

    return(round(nbins))
end

# M: data
# ndim: number of dimensions
# bintype: type of binning, select :manual for manual setting of number of bins
# bins: number of bins, only used if bintype=:manual
function DoaneRule(M,ndim)

    results = zeros(size(M))

    for i in 1:ndim

        if length(unique(M[i,:]))<3
            res_temp = zeros(length(M[i,:]))

            for u in 1:length(unique(M[i,:]))
                res_temp[M[i,:] .== unique(M[i,:])[u]] .= u
            end

            results[i,:] = res_temp
        else

            nbins = get_nbins(:doane, M[i,:])
            disc = LinearDiscretizer(binedges(DiscretizeUniformWidth(nbins), M[i,:]))
            results[i,:] = [encode(disc,one) for one in M[i,:]]

        end

    end

    return(DataFrame(results))

    return(results)

end

# function DoaneRule_cutpoints(M,ndim)
#     results = DataFrame()
#
#     for i in 1:ndim
#
#         if length(unique(M[i,:]))>5
#
#             nbins = get_nbins(:doane, M[i,:])
#             cuts = binedges(DiscretizeUniformWidth(nbins), M[i,:])
#             disc = LinearDiscretizer(cuts)
#             results[!,i] = [encode(disc,one) for one in M[i,:]]
#             results[!,i] = float(results[!,i])
#
#             for j in 1:nbins
#                 results[!,i][results[!,i].==j] .= cuts[j]
#             end
#
#
#         end
#
#     end
#
#     return(results)
#
# end


text:
# ---
# title: 1433. Check If a String Can Break Another String
# id: problem1433
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Greedy
# link: <https://leetcode.com/problems/check-if-a-string-can-break-another-string/description/>
# hidden: true
# ---
#
# Given two strings: `s1` and `s2` with the same size, check if some permutation
# of string `s1` can break some permutation of string `s2` or vice-versa. In
# other words `s2` can break `s1` or vice-versa.
#
# A string `x` can break string `y` (both of size `n`) if `x[i] >= y[i]` (in
# alphabetical order) for all `i` between `0` and `n-1`.
#
#
#
# **Example 1:**
#
#
#
#     Input: s1 = "abc", s2 = "xya"
#     Output: true
#     Explanation: "ayx" is a permutation of s2="xya" which can break to string "abc" which is a permutation of s1="abc".
#
#
# **Example 2:**
#
#
#
#     Input: s1 = "abe", s2 = "acd"
#     Output: false
#     Explanation: All permutations for s1="abe" are: "abe", "aeb", "bae", "bea", "eab" and "eba" and all permutation for s2="acd" are: "acd", "adc", "cad", "cda", "dac" and "dca". However, there is not any permutation from s1 which can break some permutation from s2 and vice-versa.
#
#
# **Example 3:**
#
#
#
#     Input: s1 = "leetcodee", s2 = "interview"
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `s1.length == n`
#   * `s2.length == n`
#   * `1 <= n <= 10^5`
#   * All strings consist of lowercase English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module Common
using Compat: String
using PiecewiseIncreasingRanges

export File, Channel, ContinuousChannel, SpikeChannel, EventChannel,
       AbstractChannels, Channels, GenericSpikeChannel, validchannels, data, datavecs,
       times, hasdata, samplerate, bitspersample, voltagemultiplier, unitnumbers, label

abstract File
abstract NeuroChannel
abstract ContinuousChannel <: NeuroChannel
abstract SpikeChannel <: NeuroChannel
abstract EventChannel <: NeuroChannel
abstract AbstractChannels{C<:NeuroChannel}

immutable Unit{W,T}
    number::Int
    times::Vector{Float64}
    data::W
    multiplier::T
end

"""
   Channels{C<:NeuroChannel}()

A container for channels of type C that is iterable but indexed by channel number
"""
immutable Channels{C<:NeuroChannel} <: AbstractChannels{C}
    channels::Vector{C}           # all channels that exist
    number2idx::Vector{Int}       # mapping between channel number and index
    idx2number::Vector{Int}       # mapping between index and channel number
    Channels() = new(C[], Int[], Int[])
end

# Iteration returns only channels that exist
Base.start(c::AbstractChannels) = 1
Base.next(c::AbstractChannels, i) = (c.channels[i], i+1)
Base.done(c::AbstractChannels, i) = i == length(c)+1
Base.eltype{C}(c::AbstractChannels{C}) = C
Base.length(c::AbstractChannels) = length(c.channels)

immutable ChannelEnumerator{C<:AbstractChannels}
    x::C
end
Base.start(c::ChannelEnumerator) = 1
Base.next(c::ChannelEnumerator, i) = ((c.x.idx2number[i], c.x.channels[i]), i+1)
Base.done(c::ChannelEnumerator, i) = i == length(c.x)+1
Base.eltype(c::ChannelEnumerator) = Tuple{Int, eltype(c.x)}
Base.length(c::ChannelEnumerator) = length(c.x)
Base.enumerate(c::AbstractChannels) = ChannelEnumerator(c)

function Base.getindex{C}(c::AbstractChannels{C}, x::AbstractVector)
    newch = Channels{C}()
    for i in x
        newch[i] = c[i]
    end
    newch
end

Base.getindex{C}(c::AbstractChannels{C}, i) = c.channels[c.number2idx[i]]
function Base.setindex!(c::AbstractChannels, v, i)
    n = length(c.number2idx)

    if i > n
        resize!(c.number2idx, i)
        c.number2idx[n+1:i] = 0
    end

    if c.number2idx[i] == 0
        push!(c.channels, v)
        push!(c.idx2number, i)
        c.number2idx[i] = length(c.channels)
    else
        c.channels[c.number2idx[i]] = v
    end
    c
end
Base.endof(c::AbstractChannels) = isempty(c.channels) ? 0 : endof(c.number2idx)
Base.haskey(c::AbstractChannels, i::Integer) =
    0 < i <= length(c.number2idx) && c.number2idx[i] != 0

Base.isassigned(c::AbstractChannels, i) =
    i <= length(c.number2idx) && c.number2idx[i] != 0

function Base.show(io::IO, c::AbstractChannels)
    write(io, "$(length(c))-channel $(typeof(c))")
    for (i, ch) in enumerate(c)
        write(io, "\n  $i => $ch")
    end
end

"""
    validchannels(c::AbstractChannels)

Get the channel numbers for all defined channels
"""
validchannels(c::AbstractChannels) = c.idx2number

"""
    NeuroIO.times(c::Union{Unit,NeuroChannel})

Get the times at which spikes, events, or recorded samples were acquired.
"""
times(c::Union{Unit,NeuroChannel}) = c.times

"""
    NeuroIO.times(c::AbstractChannels{C<:ContinuousChannel})

Get the times corresponding to recorded samples on all continuous channels, if
they had a common time basis.
"""
function times{C<:ContinuousChannel}(c::AbstractChannels{C})
    length(c) == 0 && throw(ArgumentError("no channels defined"))
    times = c.channels[1].times
    for i = 2:length(c.channels)
        if c.channels[i].times != times
            throw(ArgumentError("channels do not have a common time basis"))
        end
    end
    return times
end

"""
    NeuroIO.data(c::Union{Unit,NeuroChannel})

Get the spike waveforms (for `SpikeChannel` or `Unit`), event codes (for
`EventChannel`), or continuous data (for `ContinuousChannel`).
"""
data(c::NeuroChannel) = c.data

"""
    NeuroIO.data(c::ContinuousChannel, samples=1:length(NeuroIO.times(c)))

Get the recorded samples for the given channel. If scaled is false, the raw data
is returned. If scaled is true, samples are returned in microvolts.
"""
data(c::NeuroChannel, samples::AbstractVector{Int}) = sub(data(c), samples)

"""
    NeuroIO.data(c::AbstractChannels{C<:ContinuousChannel},
                 channels=validchannels(c),
                 samples=1:length(chtimes(c)))

Get data from continuous channels recorded with a common time basis as an
nchannels × nsamples matrix. `samples` specifies the sample indexes. `channels`
specifies the specific channels to be returned. If set to `nothing`, all
channels are returned.
"""
function data{C<:ContinuousChannel}(c::AbstractChannels{C},
                                    channels::Union{Int,AbstractVector{Int}}=validchannels(c),
                                    samples::Range{Int}=1:length(chtimes(c)))
    length(channels) == 0 && throw(ArgumentError("no channels specified"))

    # Ensure common time basis
    times = c[channels[1]].times
    for i = 2:length(channels)
        if c[channels[i]].times != times
            throw(ArgumentError("channels do not have a common time basis"))
        end
    end

    data = Array(eltype(c[1]), length(samples), length(channels))
    for i = 1:length(channels)
        copy!(data, size(data, 1)*(i-1)+1,
              sub(chdata(c[channels[i]]), samples))
    end
    data.'
end

"""
    NeuroIO.datavecs(c::AbstractChannels, channels=nothing)

Get data from multiple channels simultaneously and return a vector of vectors.
This can be more efficient than calling NeuroIO.data on each channel,
depending on the on-disk layout.
"""
datavecs(c::AbstractChannels, channels::AbstractVector{Int}=validchannels(c)) =
    [data(c[x]) for x in channels]

"""
    NeuroIO.hasdata(c::NeuroChannel)

Whether data (waveforms, event codes, or samples) is available for this channel.
"""
hasdata(c::NeuroChannel) = isdefined(c, :data)

"""
    samplerate(c::SpikeChannel)

Get the sample rate for the given spike channel.
"""
samplerate(c::SpikeChannel) = c.sample_rate

"""
    samplerate(c::ContinuousChannel)

Get the sample rate for the given continuous channel.
"""
samplerate(c::ContinuousChannel) = _samplerate(times(c))
_samplerate{T,R<:UnitRange}(r::PiecewiseIncreasingRange{T,R}) = r.divisor
_samplerate(r::PiecewiseIncreasingRange) = r.divisor//step(r.ranges[1])

"""
    bitspersample(c::Union{SpikeChannel,ContinuousChannel})

Get the number of bits encoded for a waveform or continuous sample.
"""
bitspersample(c::Union{SpikeChannel,ContinuousChannel}) = c.bits_per_sample

"""
    voltagemultiplier(c::Union{SpikeChannel,ContinuousChannel})

Get multiplier to convert raw AD values to volts, so that
`NeuroIO.data(c)*voltagemultiplier(c)` gives the recored data voltages.
"""
voltagemultiplier(c::Union{SpikeChannel,ContinuousChannel}) = c.voltage_multiplier

"""
    NeuroIO.label(c::NeuroChannel)

Get the channel label, if one exists.
"""
label(c::NeuroChannel) = c.label

"""
    unitnumbers(c::SpikeChannel)

Unit numbers corresponding to spikes on the given channel. 0 if unsorted.
"""
unitnumbers(c::SpikeChannel) = c.unit_numbers

immutable GenericSpikeChannel{T,S} <: SpikeChannel
    label::String
    sample_rate::S
    voltage_multiplier::Float64
    times::Vector{Float64}
    unit_numbers::Vector{Int16}
    data::Matrix{T}
end
bitspersample{T}(::GenericSpikeChannel{T}) = sizeof(T)*8

Base.show(io::IO, c::SpikeChannel) =
    @printf io "%s %s %d spikes" typeof(c) label(c) length(times(c))
Base.show(io::IO, c::EventChannel) =
    @printf io "%s %s %d spikes" typeof(c) label(c) length(times(c))
Base.show(io::IO, c::ContinuousChannel) =
    @printf io "%s %s %d samples @ %.1f Hz" typeof(c) label(c) length(times(c)) samplerate(c)
end


text:
function _parseString(str)
	if length(str) > 1 && str[1] == '"' && str[end] == '"'
		try
			temp::String = Meta.parse(replace(str, '$' => "\$"))
			return temp
		catch
		end
	end
	return nothing
end

function _parseChar(str)
	if length(str) > 1 && str[1] == '\'' && str[end] == '\''
		try
			temp::Char = Meta.parse(str)
			return temp
		catch
		end
	end
	return nothing
end

function _parseInteger(str)
	m = match(r"^([-])?(0|0x|0b)?([0-9A-Fa-f]+)([ulUL]*)$", lowercase(str))
	if m !== nothing
		us = count(isequal('u'), m.captures[4])
		ls = count(isequal('l'), m.captures[4])
		if us in (0, 1) && ls in (0, 1 ,2)
			typ = ls == 0 ? Cint : (ls == 1 ? Clong : Clonglong)
			typ = us == 1 || m.captures[2] !== nothing ? unsigned(typ) : typ

			prefix = m.captures[2] === nothing ? "" : (m.captures[2] == "0" ? "0o" : m.captures[2])
			try
				literal = parse(typ, prefix*m.captures[3])
				if m.captures[1] == "-"
					literal = -literal
				end
				return literal
			catch
			end
		end
	end
	return nothing
end

function _parseFloat(str)
	if !isempty(str)
		typ = str[end] in ('l', 'L') ? nothing #= TODO: use 128-bit float when available =# : str[end] in ('f', 'F') ? Cfloat : Cdouble
		value = str[end] in ('l', 'L', 'f', 'F') ? str[1:end-1] : str
		value = value[1] == '.' ? "0"*value : value
		try
			return parse(typ, value)
		catch
		end
	end
	return nothing
end

function _parseSymbol(str)
	m = match(r"^[A-Za-z0-9_]+$", str)
	if m !== nothing
		return str
	end
	return nothing
end

function _parseExpr(str)
	try
		temp = Meta.parse(str)
		Base.is_expr(temp, :incomplete) || return temp
	catch
	end
	return nothing
end

function _convertLiteral(ctx, decl::Clang.CLIntegerLiteral)
	tokens = Clang.tokenize(decl)
	length(tokens) == 1 || error("Unable to handle more than 1 token for an integer literal")
	result = _parseInteger(tokens[1])
	isnothing(result) && error("Failed to parse integer literal")
	return result
end

#=
	function _prejuliaize(ctx, decl)
		if haskey(ctx.macroUses, _extent(decl))
			name = ctx.macroUses[_extent(decl)]
			tokens = map(t -> t.text, tokenize(decl))
			args = []
			for i in 3:length(tokens)-1
				push!(args, tokens[i] == "," ? ", " : tokens[i])
			end
			return Expr(:macrocall, Symbol('@', name), LineNumberNode(@__LINE__, Symbol(@__FILE__)), (isempty(args) ? () : (:(error($(join(args)))),))...)
		else
			return _juliaize(ctx, decl)
		end
	end

	function _juliaize(ctx, decl)
		error(_where(decl)*"   "*string(string(decl) => map(string, Clang.children(decl))))
	end

	function _juliaize(ctx, decl::Clang.CLUnaryExpr)
		c = Clang.children(decl)
		length(c) == 1 || error("Expected exactly 1 children of a unary expression.")
		tokens = Clang.tokenize(decl)
		length(tokens) > 1 || error("Not enough tokens in unary expression.")

		return :($(Symbol(first(tokens).text))($(_prejuliaize(ctx, first(c)))))
	end

	function _juliaize(ctx, decl::Clang.CLUnaryOperator)
		c = Clang.children(decl)
		length(c) == 1 || error("Expected exactly 1 children of a unary operator.")
		tokens = Clang.tokenize(decl)
		length(tokens) == (length(Clang.tokenize(first(c))) + 1) || error("Tokens don't map correctly in unary operator.")

		op = Symbol(strip(first(tokens).text))
		if op === :!
			op = :cnot
		elseif op === :*
			op = :unsafe_load
		elseif op === :&
			op = :addr_of
		elseif !(op in (:-, :~))
			error("Unhandled unary operator '$(op)'.")
		end

		return :($(op)($(_prejuliaize(ctx, first(c)))))
	end

	function _juliaize(ctx, decl::Clang.CLBinaryOperator)
		c = Clang.children(decl)
		length(c) == 2 || error("Expected exactly 2 children of a binary operator.")

		tokens = Clang.tokenize(decl)
		length(tokens) == (length(Clang.tokenize(first(c))) + length(Clang.tokenize(last(c))) + 1) || error("Tokens don't map correctly in binary operator.  $(map(t -> strip(t.text), Clang.tokenize(decl)))  $(map(string, Clang.children(decl)))")

		op = Symbol(strip(tokens[length(Clang.tokenize(first(c)))+1].text))
		return Expr(:call, op, _prejuliaize(ctx, first(c)), _prejuliaize(ctx, last(c)))
	end

	function _juliaize(ctx, decl::Clang.CLCStyleCastExpr)
		c = Clang.children(decl)
		length(c) == 2 || error("Expected exactly 2 children of a cast expression.")
		return :(reinterpret($(_prejuliaize(ctx, first(c))), $(_prejuliaize(ctx, last(c)))))
	end

	function _juliaize(ctx, decl::Clang.CLStmtExpr)
		c = Clang.children(decl)
		length(c) == 1 || error("Expected exactly 1 child of a statement expression.")
		return _prejuliaize(ctx, first(c))
	end

	function _juliaize(ctx, decl::Clang.CLCallExpr)
		c = Clang.children(decl)
		length(c) >= 1 || error("Expected a function name to be used in a call expression.")
		return :($(_prejuliaize(ctx, first(c)))($(map(x -> _prejuliaize(ctx, x), c[2:end])...)))
	end

	_juliaize(ctx, decl::Clang.CLVarDecl) = nothing #:($(_c2julia(ctx, decl)) = $(_c2julia(ctx, Clang.type(decl))))  # TODO: need better way of constructing the variable
	_juliaize(ctx, decl::Clang.CLTypeRef) = _c2julia(ctx, Clang.type(decl))

	function _juliaize(ctx, decl::Clang.CLDeclStmt)
		c = Clang.children(decl)
		length(c) == 1 || error("Unable to handle declaration statements with more than 1 child.")
		return _prejuliaize(ctx, first(c))
	end

	function _juliaize(ctx, decl::Clang.CLGotoStmt)
		c = Clang.children(decl)
		length(c) == 1 || error("Unable to handle goto statements without 1 child.")
		return :(@goto $(_prejuliaize(ctx, first(c))))
	end

	function _juliaize(ctx, decl::Clang.CLLabelStmt)
		c = Clang.children(decl)
		length(c) == 1 || error("Unable to handle label statements without 1 child.")
		return :(@label $(_c2julia(ctx, decl)) ; $(_prejuliaize(ctx, first(c))))
	end

	_juliaize(ctx, decl::Clang.CLLabelRef) = _c2julia(ctx, decl)
	_juliaize(ctx, decl::Clang.CLDeclRefExpr) = _c2julia(ctx, decl)

	function _juliaize(ctx, decl::Clang.CLReturnStmt)
		length(Clang.children(decl)) == 1 || error("Wrong number of children in return statement.")
		return :(return $(_prejuliaize(ctx, first(Clang.children(decl)))))
	end

	function _juliaize(ctx, decl::Clang.CLMemberRefExpr)
		c = Clang.children(decl)
		length(c) == 1 || error("Unable to handle member reference expressions with anything other than 1 child.")
		return :($(_prejuliaize(ctx, first(c))).$(_c2julia(ctx, decl)))
	end

	function _juliaize(ctx, decl::Clang.CLUnexposedExpr)
		c = Clang.children(decl)
		length(c) == 1 || error("Unable to handle unexposed expressions with anything other than 1 child.")
		return _prejuliaize(ctx, first(c))
	end

	function _juliaize(ctx, decl::Clang.CLCompoundStmt)
		return :(let ; $(filter(!isnothing, map(x -> _prejuliaize(ctx, x), Clang.children(decl)))...) end)
	end

	function _juliaize(ctx, decl::Clang.CLIfStmt)
		c = Clang.children(decl)
		length(c) == 2 || error("Expected a function name to be used in an if statement.")
		return :(if ($(_prejuliaize(ctx, first(c)))) ; $(_prejuliaize(ctx, last(c))) end)
	end

	function _juliaize(ctx, decl::Clang.CLIntegerLiteral)
		tokens = Clang.tokenize(decl)
		length(tokens) == 1 || error("Unable to handle more than 1 token for an integer literal.  $(map(t -> strip(t.text), tokens))")
		result = _parseInteger(tokens[1])
		return result === nothing ? Symbol(tokens[1]) : result
	end

	function _juliaize(ctx, decl::Clang.CLConditionalOperator)
		c = Clang.children(decl)
		length(c) == 3 || error("Expected exactly 3 children of a ternary operator.")
		return :(($(_prejuliaize(ctx, first(c)))) ? ($(_prejuliaize(ctx, c[2]))) : ($(_prejuliaize(ctx, last(c)))))
	end

	function _juliaize(ctx, decl::Clang.CLParenExpr)
		c = Clang.children(decl)
		length(c) == 1 || error("Expected 1 child in a paren expression.")
		return :(($(_prejuliaize(ctx, first(c)))))
	end

	function _juliaize(ctx, decl::Clang.CLFunctionDecl)
		c = Clang.children(decl)
		length(c) > 1 || error("Unable to handle inline function declarations without more than 1 child.")  # NOTE:  the last child is the contents of the function
		last(c) isa Clang.CLCompoundStmt || error("Expected inline function declaration to have a compount statement for its body.")
		return filter(!isnothing, map(x -> _prejuliaize(ctx, x), Clang.children(last(c))))
	end
=#


text:
using Test, LinearAlgebra, HierarchicalMatrices

# testing BBox
include("bbox_tests.jl")

# testing BinTree
include("bintree_tests.jl")

# A-Bt decomposition
include("abtdecomposition_tests.jl")

# Hierarchical Matrices
include("hmatrix_tests.jl")

max_stars_repo_name: marcbakry/HierarchicalMatrices

max_issues_repo_name: marcbakry/HierarchicalMatrices

max_forks_repo_name: marcbakry/HierarchicalMatrices