text:
# TODO make consistent with Graphulo's degree count, which only adds 1 (like logical)
function adjbfs(A::Assoc, v0::AbstractString, numsteps::Number, minDegree = 1::Number, maxDegree = 2^31 - 1::Number; takeunion = false, nodiag = false)
    # compute the outdegrees.
    # If nodiag is true, then the diagonal of the adjacency matrix will be ignored. (Handy when the diagonal represents edge counts, rather than self-loops)
    if nodiag
        Anodiag = A-diag(A)
        Anodiag = Assoc(Anodiag.row,Anodiag.col,Anodiag.val,dropzeros!(Anodiag.A))
        Adeg = sum(Anodiag, 1)
    else
        Anodiag = A
        Adeg = sum(A, 1)
    end

    return adjbfs(Anodiag, Adeg, v0, numsteps, minDegree, maxDegree; takeunion = takeunion)

end

function adjbfs(A::Assoc, Adeg::Assoc, v0::AbstractString, numsteps::Number, minDegree = 1::Number, maxDegree = 2^31 - 1::Number; takeunion = false)

    vk = v0
    Ak = Assoc("", "", "")

    if vk == ""
        println("No start nodes.")
    end

    if takeunion
        Akall = Assoc("", "", "")
    end

    for i = 1:numsteps
        # println("step " * string((i)))
        # println("vk " * string((vk)))

        if isa(vk, Array) && ndims(vk) == 0 # meaning no nodes reachable from previous step's uk (or v0 is '').
            Ak = Assoc("", "", "")
            println("Stopped BFS after step " * string(i) *  "; no nodes in search set.")
            break
        end

        ukAssoc = bounded(Adeg[:,vk], minDegree, maxDegree)
        uk = ukAssoc.col  # filtered for degree

        # println("uk " * string(uk))
        if ndims(uk) == 0
            Ak = Assoc("", "", "")
        else
            Ak = A[uk,:]   # submatrix within reach of filtered neighbors
            if takeunion
                Akall = Ak + Akall
            end
        end

        vk = Ak.col   # nodes exactly k hops away
    end
    if takeunion
        Ak = Akall
    end

    return Ak

end


text:
include("./params.jl")

function pprocess_inhomopoisson(lambda)
    # Example usage:
    # L = 100000;
    # ts = linspace(0, 10, L);
    # f = 1.0;
    # lambda = 20*(1+sin(2.0 * pi * f * ts));
    # x, rs = pprocess_inhomopoisson(lambda)
    #
    # Input argument
    # lambda: A vector that specifies a time-dependent spike-rate.
    #
    # Output argument
    # x:    Spike time.
    # rs:   Intensity function.

    #Initialize parameters

    x = [0.0]; # push!(x, num)
    rs = zeros(L);

    R = 0;

    rs[1] = lambda[1];

    eta = - log(rand());
    for i = 2:L
        tau = i*dt - x[end]

        r = lambda[i]; # Intensity
        R = R + r*dt;

        if eta <= R
            push!(x,i*dt) # Spike time

            R = 0.0;
            eta = - log(rand());
        end
        rs[i] = r;
    end

    return x, rs
end

ext: jl

max_stars_repo_name: ganow/julia-point-process

max_issues_repo_name: ganow/julia-point-process

max_forks_repo_name: ganow/julia-point-process


text:
using FluxProjects
using Test

@testset "FluxProjects.jl" begin
    # Write your tests here.
end


text:
"""
    gpinv(p, k, sigma)

Inverse Generalised Pareto distribution function.
"""
function gpinv(p, k, sigma)
    x = similar(p)
    if sigma <= 0
        return x
    end
    ok = (p .> 0) .& (p .< 1)
    pv = view(p, ok)
    xv = view(x,ok)
    if abs(k) < eps(Float64)
        xv .= -log1p.(-pv)
    else
        xv .= expm1.(log1p.(-pv) .* (-k)) ./ k
    end
    x .*= sigma
    if ~all(ok)
        x[p .== 0] = 0
        if k >= 0
            x[p .== 1] = Inf
        else
            x[p .== 1] = -sigma / k
        end
    end
    x
end


text:
struct LazyTree
    head::String
    icon::String
    isempty::Bool
    children
end

LazyTree(head, children) = LazyTree(head, "", false, children)
LazyTree(head, icon::String, children) = LazyTree(head, icon, false, children)

struct SubTree
    head::String
    icon::String
    child
end

SubTree(head, child) = LazyTree(head, "", child)

struct Leaf
    val
    icon::String
end

Leaf(val) = Leaf(val, "")

const TREES = Dict{Int,LazyTree}()
const ID = Ref(0)

const MAX_PARTITION_LENGTH = 20

treeid() = (ID[] += 1)

pluralize(n::Int, one, more = one) = string(n, " ", n == 1 ? one : more)
pluralize(::Tuple{}, one, more = one) = string(0, " ", more)
pluralize(n, one, more = one) = string(length(n) > 1 ? join(n, '×') : first(n), " ", prod(n) == 1 ? one : more)

function treerender(x::LazyTree)
    id = treeid()
    TREES[id] = x

    return ReplWorkspaceItem(
        x.head,
        id,
        !(x.isempty),
        true,
        x.icon,
        "",
        false,
        ""
    )
end

function treerender(x::SubTree)
    child = treerender(x.child)

    return ReplWorkspaceItem(
        x.head,
        child.id,
        child.haschildren,
        child.lazy,
        child.icon,
        child.head,
        false,
        ""
    )
end

function treerender(x::Leaf)
    return ReplWorkspaceItem(
        x.val,
        -1,
        false,
        false,
        x.icon,
        "",
        false,
        ""
    )
end

getfield_safe(x, f, default = UNDEF) = isdefined(x, f) ? getfield(x, f) : default

function treerender(x)
    fields = fieldnames(typeof(x))

    if isempty(fields)
        treerender(Text(string(typeof(x), "()")))
    else
        treerender(LazyTree(string(typeof(x)), wsicon(x), function ()
            collect([SubTree(string(f), wsicon(getfield_safe(x, f)), getfield_safe(x, f)) for f in fields])
        end))
    end
end

function treerender(x::AbstractDict{K,V}) where {K,V}
    treerender(LazyTree(string(nameof(typeof(x)), "{$(K), $(V)} with $(pluralize(length(keys(x)), "element", "elements"))"), wsicon(x), length(keys(x)) == 0, function ()
        if length(keys(x)) > MAX_PARTITION_LENGTH
            partition_by_keys(x, sz = MAX_PARTITION_LENGTH)
        else
            collect([SubTree(repr(k), wsicon(v), v) for (k, v) in x])
        end
    end))
end

function treerender(x::Module)
    treerender(LazyTree(string(x), wsicon(x), function ()
        ns = names(x, all = true)
        out = []
        for n in ns
            isdefined(x, n) || continue
            Base.isdeprecated(x, n) && continue
            startswith(string(n), '#') && continue

            v = getfield(x, n)
            v === x && continue

            push!(out, SubTree(string(n), wsicon(v), v))
        end

        out
    end))
end

struct Undef end

const UNDEF = Undef()

function assign_undefs(xs)
    xs′ = similar(xs, Any)
    for i in eachindex(xs)
        xs′[i] = isassigned(xs, i) ? xs[i] : UNDEF
    end
    return xs′
end

function treerender(x::AbstractArray{T,N}) where {T,N}
    treerender(LazyTree(string(typeof(x), " with $(pluralize(size(x), "element", "elements"))"), wsicon(x), length(x) == 0, function ()
        if length(x) > MAX_PARTITION_LENGTH
            partition_by_keys(x, sz = MAX_PARTITION_LENGTH)
        else
            collect([SubTree(repr(k), wsicon(v), v) for (k, v) in zip(keys(x), vec(assign_undefs(x)))])
        end
    end))
end

treerender(x::Number) = treerender(Leaf(strlimit(repr(x), limit = 100), wsicon(x)))
treerender(x::AbstractString) = treerender(Leaf(strlimit(repr(x), limit = 100), wsicon(x)))
treerender(x::AbstractChar) = treerender(Leaf(strlimit(repr(x), limit = 100), wsicon(x)))
treerender(x::Symbol) = treerender(Leaf(strlimit(repr(x), limit = 100), wsicon(x)))
treerender(x::Nothing) = treerender(Leaf(strlimit(repr(x), limit = 100), wsicon(x)))
treerender(x::Missing) = treerender(Leaf(strlimit(repr(x), limit = 100), wsicon(x)))
treerender(x::Ptr) = treerender(Leaf(string(typeof(x), ": 0x", string(UInt(x), base = 16, pad = Sys.WORD_SIZE >> 2)), wsicon(x)))
treerender(x::Text) = treerender(Leaf(x.content, wsicon(x)))
treerender(x::Function) = treerender(Leaf(strlimit(string(x), limit = 100), wsicon(x)))
treerender(x::Type) = treerender(Leaf(strlimit(string(x), limit = 100), wsicon(x)))
treerender(x::Undef) = treerender(Leaf("#undef", wsicon(x)))

function partition_by_keys(x, _keys = keys(x); sz = 20, maxparts = 100)
    partitions = Iterators.partition(_keys, max(sz, length(_keys) ÷ maxparts))
    out = []
    for part in partitions
        head = string(repr(first(part)), " ... ", repr(last(part)))
        if length(part) > sz
            push!(out, LazyTree(head, function ()
                partition_by_keys(x, collect(part), sz = sz, maxparts = maxparts)
            end))
        else
            push!(out, LazyTree(head, function ()
                collect([SubTree(repr(k), wsicon(v), v) for (k, v) in zip(part, getindex.(Ref(x), assign_undefs(part)))])
            end))
        end
    end
    return out
end

# workspace

repl_getvariables_request(conn, params::Nothing) = Base.invokelatest(getvariables)

function getvariables()
    M = Main
    variables = []
    clear_lazy()

    for n in names(M, all = true, imported = true)
        !isdefined(M, n) && continue
        Base.isdeprecated(M, n) && continue

        x = getfield(M, n)
        x === vscodedisplay && continue
        x === VSCodeServer && continue
        x === Main && continue

        s = string(n)
        startswith(s, "#") && continue
        try
            tree = treerender(SubTree(s, wsicon(x), x))
            tree.canshow = can_display(x)
            push!(variables, tree)
        catch err
            # FIXME: This should end up in the tree view as an "error child".
            # Ref: https://github.com/julia-vscode/julia-vscode/issues/1491
            #
            # printstyled("Internal Error: ", bold = true, color = Base.error_color())
            # Base.display_error(err, catch_backtrace())
        end
    end

    return variables
end

function clear_lazy(ids = [])
    if isempty(ids)
        empty!(TREES)
    else
        for id in ids
            delete!(TREES, id)
        end
    end
end

wsicon(::Any) = "symbol-variable"
wsicon(::Module) = "symbol-namespace"
wsicon(::Function) = "symbol-method"
wsicon(::Number) = "symbol-numeric"
wsicon(::Bool) = "symbol-boolean"
wsicon(::AbstractString) = "symbol-string"
wsicon(::AbstractArray) = "symbol-array"
wsicon(::Type) = "symbol-structure"
wsicon(::AbstractDict) = "symbol-enum"
wsicon(::Exception) = "warning"
wsicon(::Undef) = "question"

# handle lazy clicks

repl_getlazy_request(conn, id::Int) = Base.invokelatest(get_lazy, id)

function get_lazy(id::Int)
    try
        if haskey(TREES, id)
            x = [treerender(x) for x in pop!(TREES, id).children()]
            return x
        else
            return [treerender(Text("[out of date result]"))]
        end
    catch err
        # FIXME: This should end up in the tree view as an "error child".
        # Ref: https://github.com/julia-vscode/julia-vscode/issues/1491
        #
        # printstyled("Internal Error: ", bold = true, color = Base.error_color())
        # Base.display_error(err, catch_backtrace())
        return []
    end
end

max_stars_repo_name: fredrikekre/julia-vscode

max_issues_repo_name: fredrikekre/julia-vscode

max_forks_repo_name: fredrikekre/julia-vscode


text:
using PkgServer, Pkg, TOML, HTTP, JSON3, Tar, Test

# You can either perform the following setup:
#  - Already-running PkgServer, located at $JULIA_PKG_SERVER
#  - PkgServer cache directory located at $PKG_SERVER_CACHE_DIR
#
# Or you can leave those blank, and we'll start up a PkgServer for you,
# running on a background process.

# If these are not set, we will attempt to auto-initiate them.
server_process = nothing
code_dir = dirname(@__DIR__)
if isempty(get(ENV, "JULIA_PKG_SERVER", "")) || isempty(get(ENV, "JULIA_PKG_SERVER_STORAGE_ROOT", ""))
    # Start a background PkgServer as a separate process
    temp_dir = mktempdir()
    ENV["JULIA_PKG_SERVER"] = "http://127.0.0.1:8000"
    ENV["JULIA_PKG_SERVER_STORAGE_ROOT"] = temp_dir
    ENV["JULIA_PKG_SERVER_FQDN"] = "starfleet-central.pkg.julialang.org"

    @info("Automatically starting local PkgServer for testing at $(ENV["JULIA_PKG_SERVER"])")
    global server_process = run(`$(Base.julia_cmd()) --project=$(code_dir) $(code_dir)/bin/run_server.jl`; wait=false)
end

# Ensure that we don't end up spending time precompiling when we install stuff; it's useless
ENV["JULIA_PKG_PRECOMPILE_AUTO"]="0"

server_url = ENV["JULIA_PKG_SERVER"]
cache_dir = joinpath(ENV["JULIA_PKG_SERVER_STORAGE_ROOT"], "cache")

# Wrap all tests in a `try` block so that we can be certain we kill the
# server, even if we fail some tests.
try
    @info("Running test suite with configuration", server_url, cache_dir)
    include("tests.jl")
finally
    if server_process !== nothing
        @info("Reaping automatically-started local PkgServer...")
        kill(server_process)
        wait(server_process)
        @info("Outputting testing PkgServer logs:")
        for f in filter(f -> endswith(f, "-pkgserver.log"), readdir("$(temp_dir)/logs"; join=true))
            println("$(f):")
            write(stdout, read(f))
        end
    end
end

max_stars_repo_name: StefanKarpinski/PkgServer

max_issues_repo_name: StefanKarpinski/PkgServer

max_forks_repo_name: StefanKarpinski/PkgServer


text:
module Device
include("device/semaphore.jl")
include("device/mp.jl")
end


text:
function fib(n::Int=0)
fn1 = 1
fn2 = 0
if (n==0)||(n==1)
    fn = n
else
    for count in 2:n
        fn = fn2 + fn1
        fn2 = fn1
        fn1 = fn
    end
end
return fn
end


text:
module GradientTests

using Base.Test, Interpolations, DualNumbers

nx = 10
f1(x) = sin((x-3)*2pi/(nx-1) - 1)
g1(x) = 2pi/(nx-1) * cos((x-3)*2pi/(nx-1) - 1)

# Gradient of Constant should always be 0
itp1 = interpolate(Float64[f1(x) for x in 1:nx-1],
            BSpline(Constant()), OnGrid())

g = Array(Float64, 1)

for x in 1:nx
    @test gradient(itp1, x)[1] == 0
    @test gradient!(g, itp1, x)[1] == 0
    @test g[1] == 0
end

# Since Linear is OnGrid in the domain, check the gradients between grid points
itp1 = interpolate(Float64[f1(x) for x in 1:nx-1],
            BSpline(Linear()), OnGrid())
itp2 = interpolate((1:nx-1,), Float64[f1(x) for x in 1:nx-1],
            Gridded(Linear()))
for itp in (itp1, itp2)
    for x in 2.5:nx-1.5
        @test_approx_eq_eps g1(x) gradient(itp, x)[1] abs(.1*g1(x))
        @test_approx_eq_eps g1(x) gradient!(g, itp, x)[1] abs(.1*g1(x))
        @test_approx_eq_eps g1(x) g[1] abs(.1*g1(x))
    end

    for i = 1:10
        x = rand()*(nx-2)+1.5
        gtmp = gradient(itp, x)[1]
        xd = dual(x, 1)
        @test_approx_eq epsilon(itp[xd]) gtmp
    end
end

# test gridded on a non-uniform grid
knots = (1.0:0.3:nx-1,)
itp_grid = interpolate(knots, Float64[f1(x) for x in knots[1]],
                       Gridded(Linear()))

for x in 1.5:0.5:nx-1.5
    @test_approx_eq_eps g1(x) gradient(itp_grid, x)[1] abs(.5*g1(x))
    @test_approx_eq_eps g1(x) gradient!(g, itp_grid, x)[1] abs(.5*g1(x))
    @test_approx_eq_eps g1(x) g[1] abs(.5*g1(x))
end

# Since Quadratic is OnCell in the domain, check gradients at grid points
itp1 = interpolate(Float64[f1(x) for x in 1:nx-1],
            BSpline(Quadratic(Periodic())), OnCell())
for x in 2:nx-1
    @test_approx_eq_eps g1(x) gradient(itp1, x)[1] abs(.05*g1(x))
    @test_approx_eq_eps g1(x) gradient!(g, itp1, x)[1] abs(.05*g1(x))
    @test_approx_eq_eps g1(x) g[1] abs(.1*g1(x))
end

for i = 1:10
    x = rand()*(nx-2)+1.5
    gtmp = gradient(itp1, x)[1]
    xd = dual(x, 1)
    @test_approx_eq epsilon(itp1[xd]) gtmp
end

# For a quadratic function and quadratic interpolation, we expect an
# "exact" answer
# 1d
c = 2.3
a = 8.1
o = 1.6
qfunc = x -> a*(x-c).^2 + o
dqfunc = x -> 2*a*(x-c)
xg = Float64[1:5;]
y = qfunc(xg)

iq = interpolate(y, BSpline(Quadratic(Free())), OnCell())
x = 1.8
@test_approx_eq iq[x] qfunc(x)
@test_approx_eq gradient(iq, x)[1] dqfunc(x)

# 2d (biquadratic)
p = [(x-1.75)^2 for x = 1:7]
A = p*p'
iq = interpolate(A, BSpline(Quadratic(Free())), OnCell())
@test_approx_eq iq[4,4] (4-1.75)^4
@test_approx_eq iq[4,3] (4-1.75)^2*(3-1.75)^2
g = gradient(iq, 4, 3)
@test_approx_eq g[1] 2*(4-1.75)*(3-1.75)^2
@test_approx_eq g[2] 2*(4-1.75)^2*(3-1.75)

iq = interpolate!(copy(A), BSpline(Quadratic(InPlace())), OnCell())
@test_approx_eq iq[4,4] (4-1.75)^4
@test_approx_eq iq[4,3] (4-1.75)^2*(3-1.75)^2
g = gradient(iq, 4, 3)
@test_approx_eq_eps g[1] 2*(4-1.75)*(3-1.75)^2 0.03
@test_approx_eq_eps g[2] 2*(4-1.75)^2*(3-1.75) 0.2

# InPlaceQ is exact for an underlying quadratic
iq = interpolate!(copy(A), BSpline(Quadratic(InPlaceQ())), OnCell())
@test_approx_eq iq[4,4] (4-1.75)^4
@test_approx_eq iq[4,3] (4-1.75)^2*(3-1.75)^2
g = gradient(iq, 4, 3)
@test_approx_eq g[1] 2*(4-1.75)*(3-1.75)^2
@test_approx_eq g[2] 2*(4-1.75)^2*(3-1.75)

A2 = rand(Float64, nx, nx) * 100
for BC in (Flat,Line,Free,Periodic,Reflect,Natural), GT in (OnGrid, OnCell)
    itp_a = interpolate(A2, (BSpline(Linear()), BSpline(Quadratic(BC()))), GT())
    itp_b = interpolate(A2, (BSpline(Quadratic(BC())), BSpline(Linear())), GT())
    itp_c = interpolate(A2, (NoInterp(), BSpline(Quadratic(BC()))), GT())
    itp_d = interpolate(A2, (BSpline(Quadratic(BC())), NoInterp()), GT())

    for i = 1:10
        x = rand()*(nx-2)+1.5
        y = rand()*(nx-2)+1.5
        xd = dual(x, 1)
        yd = dual(y, 1)
        gtmp = gradient(itp_a, x, y)
        @test length(gtmp) == 2
        @test_approx_eq epsilon(itp_a[xd, y]) gtmp[1]
        @test_approx_eq epsilon(itp_a[x, yd]) gtmp[2]
        gtmp = gradient(itp_b, x, y)
        @test length(gtmp) == 2
        @test_approx_eq epsilon(itp_b[xd, y]) gtmp[1]
        @test_approx_eq epsilon(itp_b[x, yd]) gtmp[2]
        ix, iy = round(Int, x), round(Int, y)
        gtmp = gradient(itp_c, ix, y)
        @test length(gtmp) == 1
        @test_approx_eq epsilon(itp_c[ix, yd]) gtmp[1]
        gtmp = gradient(itp_d, x, iy)
        @test length(gtmp) == 1
        @test_approx_eq epsilon(itp_d[xd, iy]) gtmp[1]
    end
end

end


text:
export loc2cs,cs2loc

function loc2cs(loc::Array{Int64,1},n::Array{Int64,1})
@inbounds cs = loc[1] + (loc[2]-1)*n[1];
return cs;
end

function loc2cs(loc1::Int64,loc2::Int64,n::Array{Int64,1})
@inbounds cs = loc1 + (loc2-1)*n[1];
return cs;
end

function cs2loc(loc::Array{Int64,1},cs_loc::Int64,n::Array{Int64,1})
@inbounds loc[1] = mod(cs_loc-1,n[1]) + 1;
@inbounds loc[2] = div(cs_loc-1,n[1])+1;
return;
end

function cs2loc(cs_loc::Int64,n::Array{Int64,1})
@inbounds loc1 = mod(cs_loc-1,n[1]) + 1;
@inbounds loc2 = div(cs_loc-1,n[1])+1;
return loc1,loc2;
end

function get2DNeighbors(loc::Array{Int64,1},n::Array{Int64,1},ans::Array{Int64,2},relevant::Array{Bool,1})
# This function fills the neighbours array ans with neighbours the "C" style, so it doesn't allocate.
@inbounds relevant[1]= loc[1]>1;
@inbounds relevant[2]= loc[1]<n[1];
@inbounds relevant[3]= loc[2]>1;
@inbounds relevant[4]= loc[2]<n[2];

@inbounds ans[1,1] = loc[1]-1;
@inbounds ans[1,2] = loc[2];
@inbounds ans[2,1] = loc[1]+1;
@inbounds ans[2,2] = loc[2];
@inbounds ans[3,1] = loc[1];
@inbounds ans[3,2] = loc[2]-1;
@inbounds ans[4,1] = loc[1];
@inbounds ans[4,2] = loc[2]+1;
return;
end

function setPaddedArr2D(arr::Array{Bool,2},loc1::Int64,loc2::Int64,val::Bool)
# this function assumes that the boundaries are OK and the value is placed in a legal location.
@inbounds arr[loc1+2,loc2+2] = val
end

function getPaddedArr2D(arr::Array{Bool,2},loc1::Int64,loc2::Int64)
# this function assumes that the boundaries are OK and the value is placed in a legal location.
@inbounds return arr[loc1+2,loc2+2]
end

function ddxLongDiff(h::Float64,n::Int64)
	dm1 = (-1/(2*h)*ones(n-1));
	dm1[n-1] = -1/h;
	d   = zeros(n);
	d[1] = -1/h;
	d[n] = 1/h;
	dp1 = (1/(2*h))*ones(n-1);
	dp1[1] = 1/h;
	D = spdiagm((dm1,d,dp1),(-1,0,1),n,n);
	return D;
end


text:
@testset "converters" begin
    i = Bobby.cvt_to_int(
        "0000000000000000000000000000000000000000000000000000000000000001")
    @test i == 0x0000000000000001

    s = Bobby.cvt_to_binary_string(0x0000000000000001)
    @test (s ==
      "0000000000000000000000000000000000000000000000000000000000000001")

    s = Bobby.cvt_to_binary_string(1)
    @test (s ==
      "0000000000000000000000000000000000000000000000000000000000000001")

    b = Bobby.set_board()
    u = Bobby.cvt_to_uint(b.p)
    @test u == 0x00ff000000000000

    white_pawns = [0,0,0,0,0,0,0,0,
                   0,0,0,0,0,0,0,0,
                   0,0,0,0,0,0,0,0,
                   0,0,0,0,0,0,0,0,
                   0,0,0,0,0,0,0,0,
                   0,0,0,0,0,0,0,0,
                   1,1,1,1,1,1,1,1,
                   0,0,0,0,0,0,0,0]
    ba = Bobby.cvt_to_bitarray(b.P)
    @test Int.(ba) == white_pawns

    squares_uint = Bobby.gen_pgn_square_to_uint_dict()
    @test haskey(squares_uint, "e2")
    @test squares_uint["e2"] == 0x0000000000000800

    squares_int, squares_pgn = Bobby.gen_pgn_square_to_int_dict()
    @test haskey(squares_int, "e2")
    @test haskey(squares_pgn, 1)
    @test squares_int["a8"] == 1
    @test squares_pgn[1] == "a8"

    squares_int_uint = Bobby.gen_int_to_uint_dict()
    @test haskey(squares_int_uint, 1)
    @test squares_int_uint[2] == 0x4000000000000000

    color = "white"
    c = Bobby.change_color(color)
    @test c == "black"
end


text:
include("algorithm.jl")

graphNodes = 3
adjMatrix = [-1 50 100; 50 -1 60; 100 60 -1]
foodNode = 3
maxIteration = 200
antsCounts = 20
pheromoneExponent = 1
desirabilityExponent = 0.1
pheromoneDepositFactor = 1
evaporationRate = 0.1

graph = generateGraph(graphNodes,adjMatrix,foodNode)

bestPath, bestFitness, graph = algorithmAntColonyOptimization(graph, maxIteration, antsCounts, pheromoneExponent,
                                            desirabilityExponent, pheromoneDepositFactor, evaporationRate)

print("Best solution $(bestPath) with length $(bestFitness)")

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
#
#
# From base, but copied here to make sure we don't fail bacause base changed
function my_gcd(a::T, b::T) where T<:Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128}
    a == 0 && return abs(b)
    b == 0 && return abs(a)
    za = trailing_zeros(a)
    zb = trailing_zeros(b)
    k = min(za, zb)
    u = unsigned(abs(a >> za))
    v = unsigned(abs(b >> zb))
    Debugger.@bp
    while u != v
        if u > v
            u, v = v, u
        end
        v -= u
        v >>= trailing_zeros(v)
    end
    r = u << k
    # T(r) would throw InexactError; we want OverflowError instead
    r > typemax(T) && throw(OverflowError("gcd($a, $b) overflows"))
    r % T
end

# make sure source code is unavailable for this method:
include_string(@__MODULE__, """
function my_gcd_noinfo(a::T, b::T) where T<:Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128}
    a == 0 && return abs(b)
    Debugger.@bp
    b == 0 && return abs(a)
    za = trailing_zeros(a)
    zb = trailing_zeros(b)
    k = min(za, zb)
    u = unsigned(abs(a >> za))
    v = unsigned(abs(b >> zb))
    while u != v
        if u > v
            u, v = v, u
        end
        v -= u
        v >>= trailing_zeros(v)
    end
    r = u << k
    # T(r) would throw InexactError; we want OverflowError instead
    r > typemax(T) && throw(OverflowError("gcd(\$a, \$b) overflows"))
    r % T
end
""", "nope.jl")

function outer(a, b, c, d)
    inner_kw(a, b; c = c)
end

function inner_kw(a, b; c = 3, d = 10)
    return a + b + c + d
end

function f_end(x)
    x = 3 + 3
    return sin(x)
end
"we don't want to see this in the source code printing"

g_invoke(a, b, c, d, e) = a + b + c + d + e
function f_invoke(x, y, z)
    return g_invoke(x..., y, z...)
end

function f_until(x)
    s = 1 + 1
    s = 2 + 2
    s = 3 + 3
    for v in x
        s += v
    end
    return s
end

function f_up_down_1(x, y)
    f_up_down_2(x, y)
end
function f_up_down_2(x, y)
    f_up_down_3(x, y)
end
function f_up_down_3(x, y)
    return x + y
end

function mysum(x)
    s = 0
    for v in x
        s += v
    end
    return s
end

@testset "UI" begin
    if Sys.isunix() && (v"1.6.0-A" <= VERSION < v"1.7.0-A")
        Debugger._print_full_path[] = false
        using TerminalRegressionTests

        function run_terminal_test(frame, commands, validation)
            TerminalRegressionTests.automated_test(joinpath(@__DIR__, validation), commands) do emuterm
            #TerminalRegressionTests.create_automated_test(joinpath(@__DIR__, validation), commands) do emuterm
                repl = REPL.LineEditREPL(emuterm, true)
                repl.interface = REPL.setup_interface(repl)
                repl.specialdisplay = REPL.REPLDisplay(repl)
                RunDebugger(frame, repl, emuterm)
            end
        end

        CTRL_C = "\x3"
        EOT = "\x4"
        UP_ARROW = "\e[A"

        run_terminal_test(@make_frame(my_gcd(10, 20)),
                          ["n\n","`", "my_gc\t\n", "a\n", UP_ARROW, UP_ARROW, UP_ARROW, CTRL_C,
                           "w add a\n", "w add sin(a)\n", "w add b\n", "w\n", "w rm 1\n", "w\n",
                           "s\n", "fr 1\n", "fr 2\n", "f 2\n", "f 1\n",
                           "bt\n", "st\n", "C", "c\n", "C", "c\n"],
                          "ui/history_gcd.multiout")

        run_terminal_test(@make_frame(outer(1, 2, 5, 20)),
                          ["s\n", "c\n"],
                          "ui/history_kw.multiout")

        run_terminal_test(@make_frame(f_end(2)),
                          ["n\n", "n\n", "n\n"],
                          "ui/history_floor.multiout")

        run_terminal_test(@make_frame(f_invoke((1,2), 3, [4,5])),
                          ["nc\n", "s\n", "c\n"],
                          "ui/history_apply.multiout")

        run_terminal_test(@make_frame(f_until([1,2,3,4,5,6,7,8])),
                          ["+", "-", "u 76\n", "u\n", "u\n", "u\n", "c\n"],
                          "ui/history_until.multiout")

        run_terminal_test(@make_frame(f_up_down_1(5, 3)),
                          ["s\n", "s\n", "up\n", "down\n", "up 2\n", "down 2\n", "c\n"],
                          "ui/history_updown.multiout")

        run_terminal_test(@make_frame(mysum([1,2,3,4,5])),
                          ["bp add mysum:97 v > 3\n", "c\n", "bp rm\n", "c\n"],
                          "ui/bp_ui.multiout")

        run_terminal_test(@make_frame(mysum(collect(1:10000))),
                          ["bt\n", "c\n"],
                          "ui/big_repr_ui.multiout")

        Debugger.break_on(:error)
        run_terminal_test(@make_frame(error("foo")),
                          ["c\n", "bt\n", "q\n"],
                          "ui/history_break_error.multiout")
        Debugger.break_off(:error)

        run_terminal_test(@make_frame(my_gcd_noinfo(10, 20)),
                        ["n\n","`", "a\n", UP_ARROW, UP_ARROW, CTRL_C, EOT],
                         "ui/history_noinfo.multiout")
    else
        @warn "Skipping UI tests"
    end
end

# Completions
test_complete(c, s) = Debugger.completions(c, s, s)

module F
    local_var = 1
    f(f_args) = f_args
end

@testset "REPL completions" begin
    frame = JuliaInterpreter.enter_call_expr(:($(F.f)(1)))
    state = dummy_state(frame)
    prov = Debugger.DebugCompletionProvider(state)

    c, _ = test_complete(prov, "local")
    @test "local_var" in c
    c, _ = test_complete(prov, "f_")
    @test "f_args" in c
end


text:
# NOTE: there's a `kmeans!` function that updates centers, maybe a candidate
# for the `update` machinery. Same for `kmedoids!`
# NOTE: if the prediction is done on the original array, just the assignment
# should be returned, unclear what's the best way of doing this.

module MLJClusteringInterface

# ===================================================================
# IMPORTS
import Clustering
import MLJModelInterface
import MLJModelInterface: Continuous, Count, Finite, Multiclass, Table, OrderedFactor,
    @mlj_model, metadata_model, metadata_pkg

using Distances

# ===================================================================
## EXPORTS
export KMeans, KMedoids, HierarchicalClustering

# ===================================================================
## CONSTANTS
# Define constants for easy referencing of packages
const MMI = MLJModelInterface
const Cl = Clustering

# Definitions of model descriptions for use in model doc-strings.
const KMeansDescription ="""
K-Means algorithm: find K centroids corresponding to K clusters in the data.
"""

const KMedoidsDescription ="""
K-Medoids algorithm: find K centroids corresponding to K clusters in the data.
Unlike K-Means, the centroids are found among data points themselves.
"""

const HierarchicalClusteringDescription ="""
Hierarchical clustering algorithms: build a dendrogram of nested clusters
by repeatedly merging or splitting clusters.
"""

const KMFields ="""
## Keywords

* `k=3`     : number of centroids
* `metric`  : distance metric to use
"""

const PKG = "MLJClusteringInterface"

####
#### KMeans
####
"""
    KMeans(; kwargs...)

$KMeansDescription

$KMFields

See also the
[package documentation](http://juliastats.github.io/Clustering.jl/latest/kmeans.html).
"""
@mlj_model mutable struct KMeans <: MMI.Unsupervised
    k::Int = 3::(_ ≥ 2)
    metric::SemiMetric = SqEuclidean()
    init = :kmpp
end

####
#### KMeans
####

function MMI.fit(model::KMeans, verbosity::Int, X)
    # NOTE: using transpose here to get a LinearAlgebra.Transpose object
    # which Kmeans can handle.
    Xarray = transpose(MMI.matrix(X))
    result = Cl.kmeans(Xarray, model.k; distance=model.metric, init = model.init)
    cluster_labels = MMI.categorical(1:model.k)
    fitresult = (result.centers, cluster_labels) # centers (p x k)
    cache = nothing
    report = (
        assignments=result.assignments, # size n
        cluster_labels=cluster_labels
    )
    return fitresult, cache, report
end

MMI.fitted_params(::KMeans, fitresult) = (centers=fitresult[1],)

function MMI.transform(model::KMeans, fitresult, X)
    # pairwise distance from samples to centers
    X̃ = pairwise(
        model.metric,
        transpose(MMI.matrix(X)),
        fitresult[1],
        dims=2
    )
    return MMI.table(X̃, prototype=X)
end

"""
    KMedoids(; kwargs...)

$KMedoidsDescription

$KMFields

See also the
[package documentation](http://juliastats.github.io/Clustering.jl/latest/kmedoids.html).
"""
@mlj_model mutable struct KMedoids <: MMI.Unsupervised
    k::Int = 3::(_ ≥ 2)
    metric::SemiMetric = SqEuclidean()
end

function MMI.fit(model::KMedoids, verbosity::Int, X)
    # NOTE: using transpose=true will materialize the transpose (~ permutedims), KMedoids
    # does not yet accept LinearAlgebra.Transpose
    Xarray = MMI.matrix(X, transpose=true)
    # cost matrix: all the pairwise distances
    cost_array = pairwise(model.metric, Xarray, dims=2) # n x n
    result = Cl.kmedoids(cost_array, model.k)
    cluster_labels = MMI.categorical(1:model.k)
    fitresult = (view(Xarray, :, result.medoids), cluster_labels) # medoids
    cache = nothing
    report = (
        assignments=result.assignments, # size n
        cluster_labels=cluster_labels
    )
    return fitresult, cache, report
end

MMI.fitted_params(::KMedoids, fitresult) = (medoids=fitresult[1],)

function MMI.transform(model::KMedoids, fitresult, X)
    # pairwise distance from samples to medoids
    X̃ = pairwise(
        model.metric,
        MMI.matrix(X, transpose=true),
        fitresult[1], dims=2
    )
    return MMI.table(X̃, prototype=X)
end

####
#### Predict methods
####

function MMI.predict(model::Union{KMeans,KMedoids}, fitresult, Xnew)
    locations, cluster_labels = fitresult
    Xarray = MMI.matrix(Xnew)
    (n, p), k = size(Xarray), model.k
    pred = zeros(Int, n)

    @inbounds for i in 1:n
        minv = Inf
        @inbounds @simd for j in 1:k
            curv = evaluate(
                model.metric, view(Xarray, i, :), view(locations, :, j)
            )
            P = curv < minv
            pred[i] = j * P + pred[i] * !P # if P is true --> j
            minv = curv * P + minv * !P # if P is true --> curvalue
        end
    end
    return cluster_labels[pred]
end

####
#### HierarchicalClustering
####

"""
    HierarchicalClustering(; kwargs...)

$HierarchicalClusteringDescription

## Keywords

* `linkage` : `:single` (default), `:average`, `:complete`, `:ward`, `:ward_presquared`
* `metric` : distance metric to use (`SqEuclidean()` by default).
* `branchorder` : algorithm to order leaves and brancges (`:r` by default, `:barjoseph` or `:optimal`)

See also
[package documentation](https://juliastats.org/Clustering.jl/stable/hclust.html).

## Example
```
using MLJ, MLJClusteringInterface, DataFrames, StatsPlots
data = DataFrame([randn(50, 2) .+ [2 2]; randn(50, 2) .+ [-3 0]], :auto)
scatter(data.x1, data.x2) # plot data
hc = MLJ.transform(HierarchicalClustering(linkage = :complete), data)
hc(h = 25) # show cluster assignments for cutting height h = 25
hc(k = 2) # show cluster assignments for number of clusters k = 2
plot(hc.dendrogram) # plot dendrogram
scatter(data.x1, data.x2, color = hc(k = 2)) # plot data with class labels
```
"""
@mlj_model mutable struct HierarchicalClustering <: MMI.Static
    linkage::Symbol = :single :: (_ ∈ (:single, :average, :complete, :ward, :ward_presquared))
    metric::SemiMetric = SqEuclidean()
    branchorder::Symbol = :r :: (_ ∈ (:r, :barjoseph, :optimal))
end
mutable struct HierarchicalClusteringResult{T}
    dendrogram::T
end
function (hc::HierarchicalClusteringResult)(; h = nothing, k = 3)
    Cl.cutree(hc.dendrogram, k = k, h = h)
end
function Base.show(io::IO, hc::HierarchicalClusteringResult)
    println(io, "Hierarchical clustering object")
    println(io, "Call it with `h = height` or `k = number of clusters` to get the cluster assignments")
    println(io, "(see also docstring of HierarchicalClustering).")
end

function MMI.transform(model::HierarchicalClustering, X)
    Xarray = MMI.matrix(X)
    # all the pairwise distances
    d = pairwise(model.metric, Xarray, dims = 1) # n x n
    dendrogram = Cl.hclust(d, linkage = model.linkage, branchorder = model.branchorder)
    HierarchicalClusteringResult(dendrogram)
end

####
#### METADATA
####

metadata_pkg.(
    (KMeans, KMedoids, HierarchicalClustering),
    name="Clustering",
    uuid="aaaa29a8-35af-508c-8bc3-b662a17a0fe5",
    url="https://github.com/JuliaStats/Clustering.jl",
    julia=true,
    license="MIT",
    is_wrapper=false
)

metadata_model(
    KMeans,
    input = MMI.Table(Continuous),
    output = MMI.Table(Continuous),
    weights = false,
    descr = KMeansDescription,
    path = "$(PKG).KMeans"
)

metadata_model(
    KMedoids,
    input = MMI.Table(Continuous),
    output = MMI.Table(Continuous),
    weights = false,
    descr = KMedoidsDescription,
    path = "$(PKG).KMedoids"
)

metadata_model(
    HierarchicalClustering,
    input = MMI.Table(Continuous),
    weights = false,
    descr = HierarchicalClusteringDescription,
    path = "$(PKG).HierarchicalClustering"
)

end # module


text:
using Gaugefields
using MPI
using LinearAlgebra
using Random

const pes = Tuple(parse.(Int64,ARGS[1:4]))
const mpi = parse(Bool,ARGS[5])

function MDtest!(snet,U,Dim,mpi=false)
    p = initialize_TA_Gaugefields(U)
    Uold = similar(U)
    substitute_U!(Uold,U)
    MDsteps = 200
    temp1 = similar(U[1])
    temp2 = similar(U[1])
    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)
    numaccepted = 0

    plaq_t = calculate_Plaquette(U,temp1,temp2)*factor

    #poly = calculate_Polyakov_loop(U,temp1,temp2)
    if get_myrank(U) == 0
        println("0 plaq_t = $plaq_t")
        #println("polyakov loop = $(real(poly)) $(imag(poly))")
    end

    numtrj = 10
    for itrj = 1:numtrj
        @time accepted = MDstep!(snet,U,p,MDsteps,Dim,Uold)
        numaccepted += ifelse(accepted,1,0)

        plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        #poly = calculate_Polyakov_loop(U,temp1,temp2)

        if get_myrank(U) == 0
            println("$itrj plaq_t = $plaq_t")
            println("acceptance ratio ",numaccepted/itrj)
            #println("polyakov loop = $(real(poly)) $(imag(poly))")
        end
    end
end

function calc_action(snet,U,p)
    NC = U[1].NC
    Sg = -evaluate_GaugeAction(snet,U)/NC
    Sp = p*p/2
    S = Sp + Sg
    return real(S)
end

function MDstep!(snet,U,p,MDsteps,Dim,Uold)
    Δτ = 1/MDsteps
    gauss_distribution!(p)
    Sold = calc_action(snet,U,p)
    substitute_U!(Uold,U)

    for itrj=1:MDsteps
        U_update!(U,p,0.5,Δτ,Dim,snet)
        #println(getvalue(U[1],1,1,1,1,1,1))

        P_update!(U,p,1.0,Δτ,Dim,snet)
        #if get_myrank(U) == 0
        #    println(getvalue(U[1],1,1,1,1,1,1))
        #    println("p = ",p[1][1,1,1,1,1])
        #    if isnan(p[1][1,1,1,1,1])
        #        error("p")
        #    end
        #end

        U_update!(U,p,0.5,Δτ,Dim,snet)
        #error("dd")
    end
    #error("end")

    Snew = calc_action(snet,U,p)
    if get_myrank(U) == 0
        println("Sold = $Sold, Snew = $Snew")
        println("Snew - Sold = $(Snew-Sold)")
    end
    ratio = min(1,exp(Snew-Sold))
    r = rand()
    if mpi
        r = MPI.bcast(r, 0, MPI.COMM_WORLD)
    end
    #println(r,"\t",ratio)

    if r > ratio
        substitute_U!(U,Uold)
        return false
    else
        return true
    end
end

function U_update!(U,p,ϵ,Δτ,Dim,snet)
    temps = get_temporary_gaugefields(snet)
    temp1 = temps[1]
    temp2 = temps[2]
    expU = temps[3]
    W = temps[4]

    for μ=1:Dim
        exptU!(expU,ϵ*Δτ,p[μ],[temp1,temp2])
        mul!(W,expU,U[μ])
        substitute_U!(U[μ],W)

    end
end

function P_update!(U,p,ϵ,Δτ,Dim,snet) # p -> p +factor*U*dSdUμ
    NC = U[1].NC
    temps = get_temporary_gaugefields(snet)
    dSdUμ = temps[end]
    factor =  -ϵ*Δτ/(NC)

    for μ=1:Dim
        calc_dSdUμ!(dSdUμ,snet,μ,U)
        #println("dSdU = ",getvalue(dSdUμ,1,1,1,1,1,1))
        mul!(temps[1],U[μ],dSdUμ) # U*dSdUμ
        Traceless_antihermitian_add!(p[μ],factor,temps[1])
    end
end

function test1()
    NX = 8
    NY = 8
    NZ = 8
    NT = 8
    Nwing = 0
    Dim = 4
    NC = 3

    if mpi
        PEs = pes#(1,1,1,2)
        U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold",mpi=true,PEs = PEs,mpiinit = false)

    else
        U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")

    end

    Random.seed!(123+get_myrank(U[1]))

    snet = GaugeAction(U)
    plaqloop = make_loops_fromname("plaquette")
    append!(plaqloop,plaqloop')
    β = 5.7/2
    push!(snet,β,plaqloop)

    #show(snet)

    MDtest!(snet,U,Dim,mpi)

    temp1 = similar(U[1])
    temp2 = similar(U[1])

    if Dim == 4
        comb = 6 #4*3/2
    elseif Dim == 3
        comb = 3
    elseif Dim == 2
        comb = 1
    else
        error("dimension $Dim is not supported")
    end

    factor = 1/(comb*U[1].NV*U[1].NC)

    listnames = [plaqloop]
    listvalues = [1]
    g = Gradientflow_general(U,listnames,listvalues,eps = 0.01)

    for itrj=1:100
        flow!(U,g)
        if itrj % 10 == 0
            @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
            println_verbose_level1(U[1],"$itrj plaq_t = $plaq_t")
            #poly = calculate_Polyakov_loop(U,temp1,temp2)
            #println_verbose_level1(U[1],"$itrj polyakov loop = $(real(poly)) $(imag(poly))")
        end
    end

    #show(snet)

    #MDtest!(snet,U,Dim,mpi)

end

test1()


text:
using Cascadia: parsehtml, Selector, nodeText
using DataFrames
using LibPQ: Connection, load!, execute

const sel_tr = Selector("body > table > tbody > tr")

"""
    cran_packages()::Vector{String}

Downloads the list of all available packages in CRAN.
"""
function cran_packages()
    download("https://cran.r-project.org/web/packages/available_packages_by_name.html",
             joinpath("data", "oss", "original", "cran", "all_pkgs.html"))
    html = parsehtml(String(read(joinpath("data", "oss", "original", "cran", "all_pkgs.html"))))
    s = Selector("body > table > tbody > tr > td > a")
    nodeText.(eachmatch(s, html.root))
end

# Get all the package names
cran_pkgs = cran_packages()

# Download the HTML files with the metadata
for pkg in cran_pkgs
    path = joinpath("data", "oss", "original", "cran", "index", "$pkg.html")
    isfile(path) || download("https://cran.r-project.org/web/packages/$pkg/index.html", path)
end

output = DataFrame([String, Union{Missing,String}, Union{Missing,String}], [:pkg, :maintainer, :email], 0)

for pkg in cran_pkgs
    html = parsehtml(String(read(joinpath("data", "oss", "original", "cran", "index", "$pkg.html"))))
    tbl_rows = eachmatch(sel_tr, html.root)
    pos = findfirst(elem -> nodeText(elem[1]) == "Maintainer:", tbl_rows)
    name_email = nodeText(tbl_rows[pos][2])
    if name_email ≠ "ORPHANED"
        name, email = split(name_email, '<')
        name = strip(name)
        email = replace(email, " at " => '@')
        email = replace(email, r"( |>)" => "")
    else
        name = missing
        email = missing
    end
    push!(output, (pkg = pkg, maintainer = name, email = email))
end
unique!(output)

conn = Connection("dbname = sdad")
execute(conn, "BEGIN;")
load!(output, conn, "INSERT INTO cran.cran_maintainers VALUES(\$1, \$2, \$3);")
execute(conn, "COMMIT;")
close(conn)

max_stars_repo_name: uva-bi-sdad/CRAN

max_issues_repo_name: uva-bi-sdad/CRAN

max_forks_repo_name: uva-bi-sdad/CRAN


text:
using PotentialLearning
using Test

include("GaN-SNAP-LAMMPS-Test.jl")


text:
module ClassifierCombiner

	using StatsBase: countmap
	using Distances
	using j4pr: countapp, countappw

	export AbstractCombiner, NoCombiner, LabelCombiner, ContinuousCombiner, VoteCombiner, WeightedVoteCombiner,
	       NaiveBayesCombiner, GeneralizedMeanCombiner, WeightedMeanCombiner, MedianCombiner, ProductCombiner,
	       combiner_train, combiner_exec

	abstract type AbstractCombiner end
	abstract type LabelCombiner <: AbstractCombiner end
	abstract type ContinuousCombiner <: AbstractCombiner end

	"""
	Combiner that does not combine, fake.
	"""
	struct NoCombiner <: AbstractCombiner end

	"""
	Majority voting label combiner.

	# Fields
	 * `L::Int` is the expected ensemble size on whose output it operates
	"""
	struct VoteCombiner <: LabelCombiner
		L::Int 						# Ensemble size
		C::Int						# Number of classes
	end

	"""
	Weighted majority label combiner.

	# Fields
	 * `L::Int` is the expected ensemble size on whose output it operates
	 * `weights::Vector{Float64}` weight vector; its length should be equal to the field `L`
	"""
	struct WeightedVoteCombiner <: LabelCombiner
		L::Int 						# Ensemble size
		C::Int						# Number of classes
		weights::Vector{Float64}			# Ensemble members weights
	end

	WeightedVoteCombiner(L::Int) = WeightedVoteCombiner(L,fill(1.0,L))

	"""
	Naive Bayes label combiner.

	# Fields
	 * `L::Int` is the expected ensemble size on whose output it operates
	 * `labels::Vector{S}` is the training labels vector
	 * `labelcount::Vector{Float64}` is the unique label count vector (used to calculate priors)
	 * `CM::Array{Float64,3}` is the confusion matrix array; one confusion matrix for each ensemble member, hence it is required that `size(CM,3)==L`
	"""
	struct NaiveBayesCombiner{S} <: LabelCombiner
		L::Int 						# Ensemble size
		labels::Vector{S}				# List of labels (should be in some order)
		labelcount::Vector{Float64}			# Counts corresponding to labels (in the same order)
		CM::Array{Float64,3}				# Confusion matrices for each classifier:
								#  - CM(k,s,i) indicates ensemble member 'i', true class 'k' and estimate 's';
								#  - k and s are indices in the sorted vector of unique labels.
	end

	"""
	BKS ("Behavious Knowledge Space") label combiner.

	# Fields
	 * `L::Int` is the expected ensemble size on whose output it operates
	 * `labels::Vector{S}` is the training labels vector
	 * `labelcount::Matrix{Int}` is the unique label count matrix
	 * `E::Matrix{S}` is matrix of ensemble decisions for the training data
	 * `default::S` is the default label
	"""
	struct BKSCombiner{S} <: LabelCombiner
		L::Int 						# Ensemble size
		labels::Vector{S}				# List of labels (should be in some order)
		labelcount::Matrix{Int}				# Class counts corresponding to label combinations (size: number of classes x unique label combinations)
		E::Matrix{S}				 	# E(i,j) is a matrix containing unique label combinations (size: number of classes  x unique label combinations)
		default::S					# The default label
	end

	"""
	Generalized mean combiner.

	# Fields
	 * `L::Int` is the expected ensemble size on whose output it operates
	 * `C::Int` is the number of classes
	 * `α::Float64` is the parameter of the generalized mean
	"""
	struct GeneralizedMeanCombiner <: ContinuousCombiner
		L::Int 						# Ensemble size
		C::Int						# Number of classes
		α::Float64					# Generalized mean parameter
	end

	"""
	Weighted generalized mean combiner.

	# Fields
	 * `L::Int` is the expected ensemble size on whose output it operates
	 * `C::Int` is the number of classes
	 * `α::Float64` is the parameter of the generalized mean
	 * `weights::Vector{Float64}` weight vector; its length should be equal to the field `L`
	"""
	struct WeightedMeanCombiner <: ContinuousCombiner
		L::Int 						# Ensemble size
		C::Int						# Number of classes
		α::Float64					# Generalized mean parameter
		weights::Vector{Float64}			# Ensemble members weights
	end

	"""
	Median combiner.

	# Fields
	 * `L::Int` is the expected ensemble size on whose output it operates
	 * `C::Int` is the number of classes
	"""
	struct MedianCombiner <: ContinuousCombiner
		L::Int 						# Ensemble size
		C::Int						# Number of classes
	end

	"""
	Median combiner.

	# Fields
	 * `L::Int` is the expected ensemble size on whose output it operates
	 * `C::Int` is the number of classes
	"""
	struct ProductCombiner <: ContinuousCombiner
		L::Int 						# Ensemble size
		C::Int						# Number of classes
	end

	"""
	Decision Template combiner.

	# Fields
	 * `L::Int` is the expected ensemble size on whose output it operates
	 * `C::Int` is the number of classes
	 * `DT::Array{Float64,3}` is the decision templates array
	 * `d::Distances.PreMetric` is the distance used to evaluate the similarity between the
	 current record decision profile and the templates
	"""
	struct DTCombiner{T<:Distances.PreMetric} <: ContinuousCombiner
		L::Int 						# Ensemble size
		C::Int						# Number of classes
		DT::Array{Float64,3}				# Decision templates
		d::T						# Distance
	end

	"""
		decision_profile_labels(x)

	Transforms the input `x` into a matrix; `x` can be a `Vector{Vector{T} where T}` in which case the result is a `Matrix{T}` where each element of `x` becomes a row in the output matrix.
	"""
	function decision_profile_labels(x::T where T<:AbstractVector{<:AbstractVector{S}}) where S
		# Check that all dimensions of inner vectors match
		@assert all(i == j for i in length.(x), j in length.(x)) "[Classifier Combiner] The dimensions of the vectors to be combined must match."

		n = length(x[1])	# Number of samples
		m = length(x)		# Size of ensemble
		out = Matrix{S}(m,n)	# Output

		for (i,v) in enumerate(x)
			out[i,:] = v
		end

		return out
	end

	# Desired form already, do nothing
	function decision_profile_labels(x::T where T<:AbstractMatrix{S}) where S
		x
	end

	"""
		decision_profile_continuous(x [,C])

	Transforms the input `x` into an `Array{T,3}`; `x` can be `Vector{Vector{T} where T}`, `Vector{Matrix{T} where T}` or `Matrix{T}`; `C` is the number of classes and can be either
	inferred from the sizes of `x` as:
	 * `1` if the input is a vector of vectors
	 * `size(x[1],1)` (the number of rows of the first matrix), if the input is a vector of matrices; ofcourse, the number of rows has to be the same for all matrices
	 * if the input is a `Matrix{T}`, `C` has to be explicitly specified as the ensemble outputs are already stacked together as rows of the matrix.

	So, the first dimension of the output is the number of classes, the second the number of samples and the third, the number of ensemble memebers.
	"""
	function decision_profile_continuous(x::T where T<:AbstractVector{<:AbstractVector{S}}, C::Int=1) where S # C is ignored.
		# Check that all dimensions of inner vectors match
		@assert all(i == j for i in length.(x), j in length.(x)) "[Classifier Combiner] The dimensions of vectors to be combined must match."

		n = length(x[1])	# Number of samples
		m = length(x)		# Size of ensemble
		out = Array{S,3}(1,n,m)	# Output

		@inbounds for (i,v) in enumerate(x)
			out[1,:,i] = v
		end

		return out
	end

	function decision_profile_continuous(x::T where T<:AbstractVector{<:AbstractMatrix{S}}, C::Int=size(x[1],1)) where S # C is ignored here as well
		# Check that all dimensions of inner vectors match
		@assert all(i == j for i in size.(x), j in size.(x)) "[Classifier Combiner] The dimensions of matrices to be combined must match."

		n = size(x[1],2)	# Number of samples
		m = length(x)		# Size of ensemble
		out = Array{S,3}(size(x[1],1),n,m) # Output

		@inbounds for (i,v) in enumerate(x)
			out[:,:,i] = v
		end

		return out
	end

	function decision_profile_continuous(x::T where T<:AbstractMatrix{S}, C::Int) where S
		# Check that all dimensions of inner vectors match
		@assert mod(size(x,1),C) == 0 "[Classifier Combiner] The first dimension of the input matrix should be a multiple of the number of classes."

		n = size(x,2)		# Number of samples
		m = div(size(x,1),C)	# Size of ensemble
		out = Array{S,3}(C,n,m)	# Output

		@inbounds for (i,ki) in enumerate(1:C:size(x,1))
			out[:,:,i] = x[ki:ki+C-1,:]
		end

		return out
	end

	"""
		combiner_train(predictions, labels, combiner_type)

	Trains a label combiner based on multiple ensemble outputs and returns a combiner object of `combiner_type`.

	# Arguments
	 * `predictions` are the estimated labels or continuous outputs; can be either a `Vector{Vector}` or `Matrix` for label outputs or, a `Vector{Matrix}` or `Matrix` for continuous outputs
	 * `labels` is a `Vector` containing the training labels
	 * `combiner_type` is the type of the resulting combiner and must be `<:ClassifierCombiner.AbstractCombiner`

	 It is assumed that the input size is consistent (e.g. if a `Vector{Vector}` the lengths of the inner vectors are equal), and that
	 the type of the elements of the input matches the element type of the training labels.
	"""
	# Function for training the vote combiners
	function combiner_train(predictions::T where T<:AbstractArray, labels::S where S<:AbstractVector, ::Type{NoCombiner})
		NoCombiner()
	end

	function combiner_train(predictions::T where T<:AbstractArray, labels::S where S<:AbstractVector, ::Type{VoteCombiner})
		_ = decision_profile_labels(predictions) # fake, just to check internal consistency of the input
		VoteCombiner(length(predictions), length(unique(labels)))
	end

	function combiner_train(predictions::T where T<:AbstractArray, labels::S where S<:AbstractVector, ::Type{WeightedVoteCombiner})
		ip = decision_profile_labels(predictions)
		m = size(ip,1)		# Size of ensemble
		n = size(ip,2)		# Number of samples
		w = zeros(m) 		# Weights
		c = length(unique(labels)) #number of classes

		@assert n==length(labels) "[Classifier Combiner] The size of the training labels must match the number of output samples."

		@inbounds @simd for i in 1:m
			@fastmath p = sum(isequal.(ip[i,:], labels))/n # individual ensemble member accuracy
			@fastmath w[i] = log((p+eps())/(1-p+eps()))
		end
		return WeightedVoteCombiner(m, c, w)
	end

	function combiner_train(predictions::T where T<:AbstractArray, labels::S where S<:AbstractVector{V}, ::Type{NaiveBayesCombiner}) where V
		ip = decision_profile_labels(predictions)
		m = size(ip,1)						# Size of ensemble
		n = size(ip,2)						# Number of samples
		ulabels::Vector{V} = sort(unique(labels))		# Sorted unique labels
		labelcount::Vector{Float64} = countapp(labels, ulabels)	# Label count
		C::Int = length(ulabels)				# Number of classes
		CM::Array{Float64,3} = zeros(Float64,C,C,m)		# Confusion matrices (for all ensemble members)

		@assert n==length(labels) "[Classifier Combiner] The size of the training labels must match the number of output samples."

		# Build confusion matrices
		@inbounds @simd for i in 1:m 		# for each ensemble memeber
			for s in 1:C			# for each "estimated" class
				for k in 1:C		# for each "true" class
					@fastmath CM[k,s,i] = 1.0/(labelcount[k]+1.0) * (sum( isequal.(labels, ulabels[k]) .& isequal.(ip[i,:], ulabels[s]) ) + 1.0/C)
				end
			end
		end

		return NaiveBayesCombiner(m, ulabels, labelcount, CM)
	end

	function combiner_train(predictions::T where T<:AbstractArray, labels::S where S<:AbstractVector{V}, ::Type{BKSCombiner}) where V
		ip = decision_profile_labels(predictions)
		m = size(ip,1)						# Size of ensemble
		n = size(ip,2)						# Number of samples
		ulabels::Vector{V} = sort(unique(labels))		# Sorted unique labels
		C::Int = length(ulabels)				# Number of classes

		@assert n==length(labels) "[Classifier Combiner] The size of the training labels must match the number of output samples."

		E = unique(ip,2)  					# unique decision profile i.e. only unique combinations of labels

		labelcount = zeros(Int, C,size(E,2))			# label counts for each class

		# Search for unique label output combinations and count how many fall in each class
		@inbounds for i in 1:size(ip,2) # i - samples
			for j in 1:size(E,2) 	# j - label combinations
				if view(E,:,j) == view(ip,:,i)
					labelcount[findin(ulabels,view(labels,i)),j].+=1
					continue # skip search for other combinations
				end
			end
		end

		return BKSCombiner(m, ulabels, labelcount, E, ulabels[indmax(sum(labelcount,2))])
	end

	function combiner_train(predictions::T where T<:AbstractArray, labels::S where S<:AbstractVector{V}, d::Distances.PreMetric, ::Type{DTCombiner}) where V
		ulabels::Vector{V} = sort(unique(labels))		# Sorted unique labels
		C::Int = length(ulabels)				# Number of classes

		ip = decision_profile_continuous(predictions, C)
		m = size(ip,3)						# Size of ensemble
		n = size(ip,2)						# Number of samples

		@assert n==length(labels) "[Classifier Combiner] The size of the training labels must match the number of output samples."

		# decision templates: 1'st dimension - mean class posteriors, 2'nd dimension - ensemble member, 3'rd dimension true class,
		DT = zeros(C,m,C)
		@inbounds for i in 1:C	# Calculate the mean estimated class probabilities over all ensemble members, for each class
			DT[:,:,i] = squeeze(mean(ip[:,labels.==ulabels[i],:],2),2)
		end

		return DTCombiner(m, C, DT, d)
	end

	"""
		combiner_exec(combiner, predictions)

	Combines the predictions from multiple ensemble members using a combiner returning the combined predictions

	# Arguments
	 * `combiner` is an combiner object of type `<:ClassifierCombiner.AbstractCombiner` used to combine `predictions`
	 * `predictions` are the estimated labels or continuous outputs; can be either a `Vector{Vector}` or `Matrix` for label outputs or, a `Vector{Matrix}` or `Matrix` for continuous outputs

	 It is assumed that the input size is consistent (e.g. if a `Vector{Vector}` the lengths of the inner vectors are equal), that
	 the type of the elements of the input matches the element type of the training labels and that the first dimension of `predictions` matched the learned ensemble size, `combiner.L`.
	"""
	# Functions for the execution of the combiners
	function combiner_exec(x::NoCombiner, predictions::T where T<:AbstractArray)
		predictions
	end

	function combiner_exec(x::VoteCombiner, predictions::T where T<:AbstractArray)
		p = decision_profile_labels(predictions)
		@assert size(p,1) == x.L "[Classifier combiner] Mismatch between expected and actual ensemble results size"

		out = Vector{eltype(p)}(size(p,2))	# Output
		pu = unique(p)				# Unique predicted labels
		@inbounds for i in 1:size(p,2)
			@fastmath out[i] = pu[findmax(countapp(view(p,:,i), pu))[2]]
		end

		return out
	end

	function combiner_exec(x::WeightedVoteCombiner, predictions::T where T<:AbstractArray)
		p = decision_profile_labels(predictions)
		m = size(p,1)				# Ensemble size
		n = size(p,2)				# Number of samples
		@assert m == x.L == length(x.weights) "[Classifier combiner] Mismatch between expected and actual ensemble results size or, weight vector size"

		out = Vector{eltype(p)}(n)		# Output
		pu = unique(p)				# Unique predicted labels
		@inbounds for i in 1:n
			@fastmath out[i] = pu[findmax(countappw(view(p,:,i), pu, x.weights, -Inf))[2]]
		end

		return out
	end

	function combiner_exec(x::S where S<:NaiveBayesCombiner{B}, predictions::T where T<:AbstractArray) where B
		p = decision_profile_labels(predictions)
		@assert size(p,1) == x.L "[Classifier combiner] Mismatch between expected and actual ensemble results size"

		n = size(p,2)				# Number of samples
		classes::Vector{B} = x.labels		# Ordered observed classes (in training)
		C = length(classes)			# Number of classes
		priors::Vector{Float64} = x.labelcount./sum(x.labelcount) # Vector of class priors
		out = Vector{B}(n)			# Output

		# Function that searches fast a value through a vector of unique values
		function _fastsearch_(v::T,uv::Vector{T})::Int where T
			@inbounds for (i,vi) in enumerate(uv)
				if isequal(vi,v)
					return i
				end
			end
			return 0
		end

		tmp = Vector{Float64}(C)		# Temporary class-wise priors
		@inbounds for j in 1:n			# for each sample,
			tmp = priors			# initialize posterior
			for k in 1:C			# for each "target" class
				for i in 1:size(x.CM,3)	# loop through confusion matrices
 					# Search for the estimated class (from each ensemble member) in the unique class vector;
					# use its position to retrieve the corresponding likelihood from the learned confusion matrices
					tmp[k] *= x.CM[k,_fastsearch_(p[i,j],classes)[1],i]
				end
			end
			out[j] = classes[findmax(tmp)[2]]
		end

		return out
	end

	function combiner_exec(x::S where S<:BKSCombiner{B}, predictions::T where T<:AbstractArray) where B
		p = decision_profile_labels(predictions)
		@assert size(p,1) == x.L "[Classifier combiner] Mismatch between expected and actual ensemble results size"

		n = size(p,2)				# Number of samples
		classes::Vector{B} = x.labels		# Ordered observed classes (in training)
		C = length(classes)			# Number of classes
		out = fill(x.default,n)			# Output

		@inbounds for i in 1:size(p,2)    	# i - samples
			for j in 1:size(x.E,2) 		# j - label combinations
				if view(x.E,:,j) == view(p,:,i)
					out[i] = x.labels[indmax(view(x.labelcount,:,j))]
					continue 	# skip search for other combinations
				end
			end
		end
		return out
	end

	function combiner_exec(x::GeneralizedMeanCombiner, predictions::T where T<:AbstractArray)
		p = decision_profile_continuous(predictions, x.C)
		@assert size(p,1) == x.C "[Classifier combiner] Mismatch between expected and input number of classes"
		@assert size(p,3) == x.L "[Classifier combiner] Mismatch between expected and input ensemble size"

		alpha::Float64=x.α
		return (squeeze(mean((p.+eps()).^alpha,3),3)).^(1/alpha)
	end

	function combiner_exec(x::WeightedMeanCombiner, predictions::T where T<:AbstractArray)
		p = decision_profile_continuous(predictions, x.C)
		@assert size(p,1) == x.C "[Classifier combiner] Mismatch between expected and input number of classes"
		@assert size(p,3) == x.L "[Classifier combiner] Mismatch between expected and input ensemble size"
		@assert size(p,3) == length(x.weights) "[Classifier combiner] Mismatch between ensemble size and weights length"
		alpha::Float64=x.α
		return mean(((p[:,:,i]*x.weights[i].+eps()).^alpha for i in 1:length(x.weights))).^(1.0/alpha)
	end

	function combiner_exec(x::MedianCombiner, predictions::T where T<:AbstractArray)
		p = decision_profile_continuous(predictions, x.C)
		@assert size(p,1) == x.C "[Classifier combiner] Mismatch between expected and input number of classes"
		@assert size(p,3) == x.L "[Classifier combiner] Mismatch between expected and input ensemble size"

		return Matrix{eltype(p)}(squeeze(median(p,3),3)) # median returns Float64, make sure original type is preserved
	end

	function combiner_exec(x::ProductCombiner, predictions::T where T<:AbstractArray)
		p = decision_profile_continuous(predictions, x.C)
		@assert size(p,1) == x.C "[Classifier combiner] Mismatch between expected and input number of classes"
		@assert size(p,3) == x.L "[Classifier combiner] Mismatch between expected and input ensemble size"

		return squeeze(prod(p,3),3)
	end

	function combiner_exec(x::DTCombiner, predictions::T where T<:AbstractArray)
		p = decision_profile_continuous(predictions, x.C)
		@assert size(p,1) == x.C "[Classifier combiner] Mismatch between expected and input number of classes"
		@assert size(p,3) == x.L "[Classifier combiner] Mismatch between expected and input ensemble size"
		n = size(p,2)
		out = Matrix{Float64}(x.C, n)

		for j in 1:n
			for i in 1:x.C
				out[i,j] = 1.0 .- 1.0/(x.C*x.L)*Distances.evaluate(x.d, x.DT[:,:,i], p[:,j,:])
			end
		end
		return out
	end
end

########################################
# FunctionCell and DataCell Interfaces #
########################################

## Label combiners

"""
	votecombiner(L::Int [,C::Int])

Majority vote label combiner. `L` is the size of the upstream ensemble. The number of classes `C` is most of the time not necessary are defaults to `-1`.
"""
votecombiner(L::Int, C::Int=-1) = FunctionCell(genericcombiner, Model(ClassifierCombiner.VoteCombiner(L,C), ModelProperties(L,1)), "Vote combiner")

"""
	wvotecombiner(L::Int)

Generates an untrained function cell that when piped data into, trains a weighted vote label combiner. `L` is the size of the upstream ensemble.
"""
wvotecombiner(L::Int) = FunctionCell(wvotecombiner, (L,), ModelProperties(L,1), "Weighted-Vote combiner" ) # untrained function cell

"""
	wvotecombiner(x,L)

Trains a weighted vote label combiner using the data `x`, considering an upstream ensemble size `L`.
"""
# Training
wvotecombiner(x::T where T<:CellDataL, L::Int) = wvotecombiner((getx!(x), gety(x)), L)
wvotecombiner(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractVector, L::Int) = wvotecombiner((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]), L)
wvotecombiner(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractVector, L::Int) = begin

	@assert nobs(x[1]) == nobs(x[2]) "[wvotecombiner] Expected $(nobs(x[1])) labels/values, got $(nobs(x[2]))."

	# Train combiner
	wvcombiner = ClassifierCombiner.combiner_train(getobs(x[1]), getobs(x[2]), ClassifierCombiner.WeightedVoteCombiner)

	# Check that the specified ensemble number is equal to the determined one
	@assert L==wvcombiner.L "[wvotecombiner] Expected ensemble dimension is $L and the one determined from input data is $(wvcombiner.L)"

	FunctionCell(genericcombiner, Model(wvcombiner, ModelProperties(L,1)), "Weighted-Vode combiner")
end

"""
	naivebayescombiner(L::Int)

Generates an untrained function cell that when piped data into, trains a Naive Bayes label combiner. `L` is the size of the upstream ensemble.
"""
naivebayescombiner(L::Int) = FunctionCell(naivebayescombiner, (L,), ModelProperties(L,1), "Naive-Bayes combiner" ) # untrained function cell

"""
	naivebayescombiner(x,L)

Trains a Naive Bayes label combiner using the data `x`, considering an upstream ensemble size `L`.
"""
# Training
naivebayescombiner(x::T where T<:CellDataL, L::Int) = naivebayescombiner((getx!(x), gety(x)), L)
naivebayescombiner(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractVector, L::Int) = naivebayescombiner((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]), L)
naivebayescombiner(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractVector, L::Int) = begin

	@assert nobs(x[1]) == nobs(x[2]) "[naivebayescombiner] Expected $(nobs(x[1])) labels/values, got $(nobs(x[2]))."

	# Train combiner
	nbcombiner = ClassifierCombiner.combiner_train(getobs(x[1]), getobs(x[2]), ClassifierCombiner.NaiveBayesCombiner)

	# Check that the specified ensemble number is equal to the determined one
	@assert L==nbcombiner.L "[naivebayescombiner] Expected ensemble dimension is $L and the one determined from input data is $(nbcombiner.L)"

	FunctionCell(genericcombiner, Model(nbcombiner, ModelProperties(L,1)), "Naive-Bayes combiner")
end

"""
	bkscombiner(L::Int)

Generates an untrained function cell that when piped data into, trains a BKS label combiner. `L` is the size of the upstream ensemble.
"""
bkscombiner(L::Int) = FunctionCell(bkscombiner, (L,), ModelProperties(L,1), "BKS combiner" ) # untrained function cell

"""
	bkscombiner(x,L)

Trains a BKS label combiner using the data `x`, considering an upstream ensemble size `L`.
"""
# Training
bkscombiner(x::T where T<:CellDataL, L::Int) = bkscombiner((getx!(x), gety(x)), L)
bkscombiner(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractVector, L::Int) = bkscombiner((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]), L)
bkscombiner(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractVector, L::Int) = begin

	@assert nobs(x[1]) == nobs(x[2]) "[bkscombiner] Expected $(nobs(x[1])) labels/values, got $(nobs(x[2]))."

	# Train combiner
	bkscombiner = ClassifierCombiner.combiner_train(getobs(x[1]), getobs(x[2]), ClassifierCombiner.BKSCombiner)

	# Check that the specified ensemble number is equal to the determined one
	@assert L==bkscombiner.L "[bkscombiner] Expected ensemble dimension is $L and the one determined from input data is $(bkscombiner.L)"

	FunctionCell(genericcombiner, Model(bkscombiner, ModelProperties(L,1)), "BKS combiner")
end

# Continuous output combiners

"""
	meancombiner(L::Int, C::Int [;α=1.0])

Trains a mean continuous output combiner. `L` is the size of the upstream ensemble, `C` is the expected number of classes and `α` is the parameter of the generalized mean formula.
"""
meancombiner(L::Int, C::Int;α::Float64=1.0) =
	FunctionCell(genericcombiner, Model(ClassifierCombiner.GeneralizedMeanCombiner(L,C,α), ModelProperties(L*C,C)),
	      kwtitle("Generalized mean combiner",((:α,α),)))

"""
	wmeancombiner(L::Int, C::Int, weights::Vector{Float64} [;α=1.0])

Trains a weighted mean continuous output combiner. `L` is the size of the upstream ensemble, `C` is the expected number of classes, `weights` are the
individual weights of the ensemble memebers and `α` is the parameter of the generalized mean formula.
"""
wmeancombiner(L::Int, C::Int, weights::Vector{Float64};α::Float64=1.0) =
	FunctionCell(genericcombiner, Model(ClassifierCombiner.WeightedMeanCombiner(L,C,α,weights), ModelProperties(L*C,C)),
	      kwtitle("Weighted mean combiner",((:α,α),)))

"""
	productcombiner(L::Int, C::Int)

Trains a product continuous output combiner. `L` is the size of the upstream ensemble, `C` is the expected number of classes.
"""
productcombiner(L::Int, C::Int) = FunctionCell(genericcombiner, Model(ClassifierCombiner.ProductCombiner(L,C), ModelProperties(L*C,C)), "Product combiner")

"""
	mediancombiner(L::Int, C::Int)

Trains a median continuous output combiner. `L` is the size of the upstream ensemble, `C` is the expected number of classes.
"""
mediancombiner(L::Int, C::Int) = FunctionCell(genericcombiner, Model(ClassifierCombiner.MedianCombiner(L,C), ModelProperties(L*C,C)), "Median combiner")

"""
	dtcombiner(L::Int, C::Int [;d=Distances.Euclidean()])

Generates an untrained function cell that when piped data into, trains a DT (decision templates) continuous combiner.
`L` is the size of the upstream ensemble, `C` is the expected number of classes and `d::Distances.PreMetric` is a
distance used to evaluate the similarity between the ensemble outputs for a given observation and the class templates.
"""
dtcombiner(L::Int, C::Int; d::Distances.PreMetric=Distances.Euclidean()) =
	FunctionCell(dtcombiner, (L,C), ModelProperties(L*C,C), "DT combiner";d=d) # untrained function cell

"""
	dtcombiner(x, L::Int, C::Int, [;d=Distances.Euclidean()])

Trains a DT continuous output combiner. `L` is the size of the upstream ensemble, `C` is the expected number of classes and
`d::Distances.PreMetric` is a distance used to evaluate the similarity between the ensemble outputs for a given observation
and the class templates.
"""
# Training
dtcombiner(x::T where T<:CellDataL, L::Int, C::Int; d::Distances.PreMetric=Distances.Euclidean()) =
	dtcombiner((getx!(x), gety(x)), L, C; d=d)
dtcombiner(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractVector, L::Int, C::Int; d::Distances.PreMetric=Distances.Euclidean()) =
	dtcombiner((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]), L, C; d=d)
dtcombiner(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractVector, L::Int, C::Int; d::Distances.PreMetric=Distances.Euclidean()) = begin

	@assert nobs(x[1]) == nobs(x[2]) "[dtcombiner] Expected $(nobs(x[1])) labels/values, got $(nobs(x[2]))."

	# Train combiner
	dtcombiner = ClassifierCombiner.combiner_train(getobs(x[1]), getobs(x[2]), d, ClassifierCombiner.DTCombiner)

	# Check that the specified ensemble number is equal to the determined one
	@assert L==dtcombiner.L "[dtcombiner] Expected ensemble dimension is $L and the one determined from input data is $(dtcombiner.L)"
	@assert C==dtcombiner.C "[dtcombiner] Expected $C classes and the number determined from input data is $(dtcombiner.C)"

	FunctionCell(genericcombiner, Model(dtcombiner, ModelProperties(L*C,C)), "DT combiner")
end

# Execution methods for all combiners ;)
genericcombiner(x::T where T<:CellData, model::Model{<:ClassifierCombiner.AbstractCombiner}) =
	datacell(genericcombiner(getx!(x), model), gety(x))
genericcombiner(x::T where T<:AbstractVector, model::Model{<:ClassifierCombiner.AbstractCombiner}) =
	genericcombiner(mat(x, LearnBase.ObsDim.Constant{2}()), model)
genericcombiner(x::T where T<:AbstractMatrix, model::Model{<:ClassifierCombiner.AbstractCombiner}) =
	mat(ClassifierCombiner.combiner_exec(model.data, getobs(x)), LearnBase.ObsDim.Constant{2}())


text:
# Simulate a 2D actomyosin network
# Alex Tam, 12/10/2020

# Simulate network
"Control function for actomyosin network simulations"
function actomyosin_network(parN, parA, parM, par, trial)
    # Specify domain width
    Lxx::Float64 = 2.5; Lxy::Float64 = 0; Lyx::Float64 = 0; Lyy::Float64 = 2.5; # Actual domain widths
    # Pre-allocate
    Force = [[0.0, 0.0, 0.0, 0.0] for idx in 1:parN.nT]; # [pN] Network force
    Curvature = [0.0 for idx in 1:parN.nT]; # Mean network curvature
    Index = [0.0 for idx in 1:parN.nT]; # Mean two-filament index
    Filament_Speed = [0.0 for idx in 1:parN.nT]; # [μm/s] Mean filament node speed
    Motor_Speed = [0.0 for idx in 1:parN.nT]; # [μm/s] Mean motor head speed
    Angle_ROC = [0.0 for idx in 1:parN.nT]; # [rad/s] Mean rate of change of angle
    Motor_Pos = [0.0 for idx in 1:parN.nT]; # [--] Mean myosin motor head position
    Motor_Angle = [0.0 for idx in 1:parN.nT]; # [rad] Mean angle between filaments at motor binding site
    Bins = Vector{Float64}(); Hist = Vector{Int}(); # [μm, --] Pre-allocate paired-distance histogram
    # Generate initial conditions
    state = State{Float64}(Vector{Vector{Vector}}(), Vector{Vector}()); # Initialise empty State struct
    mm = Vector{Myosin_Motor}(); # Pre-allocate empty myosin motors
    af, state = actin_ic(state, parN, parA, Lxx, Lyy); # Initialise actin filaments
    xl = intersection_search(state, af, mm); # Initialise cross-links
    mm, xl, state = myosin_ic(state, mm, parN, xl, Lxx, Lxy, Lyx, Lyy); # Initialise myosin motors
    state_old = state; # Store initial state to compute force
    # Time-stepping
    animation = @animate for i = 1:parN.nT
        # Write final DOF vector to file
        if i == parN.nT
            dof = build_dof(state);
            writedlm("dof-$par-trial-$trial-$i.csv", dof);
        end
        # Compute network characteristics
        Curvature[i] = curvature(af, state, Lxx, Lxy, Lyx, Lyy); # savefig("actomyosin_curvature-$par-trial-$trial.svg"); # Mean filament curvature
        Index[i] = two_filament_index(mm, state, Lxx, Lxy, Lyx, Lyy); # savefig("actomyosin_2f_index-$par-trial-$trial.svg"); # Mean two-filament index
        Bins, Hist = pcf(af, state, Lxx, Lyy); savefig("actomyosin_pcf-$par-trial-$trial.svg"); # Paired distances between filament nodes
        Filament_Speed[i] = filament_speed(parN, af, state, state_old, Lxx, Lxy, Lyx, Lyy); # savefig("actomyosin_filament_speed-$par-trial-$trial.svg"); # Filament node speeds
        Motor_Speed[i] = motor_speed(parN, mm, state, state_old, Lxx, Lxy, Lyx, Lyy); # savefig("actomyosin_motor_speed-$par-trial-$trial.svg"); # Motor head speeds
        Angle_ROC[i] = motor_angle_roc(parN, mm, state, state_old, Lxx, Lxy, Lyx, Lyy); # savefig("actomyosin_angle_roc-$par-trial-$trial.svg"); # Angle rate of change
        Motor_Pos[i] = motor_position(mm, state); # savefig("actomyosin_angle_roc-$par-trial-$trial.svg"); # Myosin motor head positions
        Motor_Angle[i] = motor_angle(mm, state, parN, Lxx, Lxy, Lyx, Lyy); # savefig("actomyosin_angle_roc-$par-trial-$trial.svg"); # Myosin motor angles
        Force[i] = network_force(state, state_old, af, xl, mm, parN, parA, parM, Lxx, Lxy, Lyx, Lyy); # Force
        # Draw network
        draw_network(state, af, xl, mm, parN, parA, Force[i], Lxx, Lxy, Lyx, Lyy);
        if i == 1
            savefig("actomyosin_ic-par-$par-trial-$trial.svg"); # Save image of initial condition
        end
        # Simulate one time step
        if i != parN.nT # Ensure correct looping sequence
            # Turnover and polymerisation
            af, mm, state = actin_turnover(state, af, mm, parN, parA, Lxx, Lyy);
            af = segment_translations(state, af); # Update filament translations
            xl = intersection_search(state, af, mm); # Update cross-links
            mm, xl, state = myosin_turnover(state, xl, mm, parN, parM, Lxx, Lxy, Lyx, Lyy); # Myosin binding/unbinding
            # Compute solution and store data
            state_old = state; # Store current state for energy functional
            new_dof = optimise_network(state_old, af, xl, mm, parN, parA, parM, Lxx, Lxy, Lyx, Lyy); # Compute network solution
            state = build_state(new_dof, af, mm); # Construct State from vector
        end
    end
    # Output .gif and image of the network
    gif(animation, "actomyosin-par-$par-trial-$trial.gif", fps = 10); savefig("actomyosin_end-par-$par-trial-$trial.svg");
    # Plot quantities versus time
    draw_stress(parN, Force, parN.nT, Lxx, Lyy); savefig("actomyosin_stress-par-$par-trial-$trial.svg"); # Stress components
    times, Bulk_Stress, Bulk_Stress_Int = draw_bulk_stress(parN, Force, parN.nT, Lxx, Lyy); savefig("actomyosin_bulk_stress-par-$par-trial-$trial.svg"); # Bulk stress
    @printf("Time-averaged bulk stress is %f pN/μm*s.\n", Bulk_Stress_Int[end]/times[end])
    # Compute time-integrated network statistics
    Curvature_Int, Index_Int, Filament_Speed_Int, Motor_Speed_Int, Angle_ROC_Int, Motor_Pos_Int, Motor_Angle_Int = integrated_statistics(parN, Curvature, Index, Filament_Speed, Motor_Speed, Angle_ROC, Motor_Pos, Motor_Angle)
    # Save mean stress and spatial measures time series data to files
    writedlm("times.csv", times);
    writedlm("stress-par-$par-trial-$trial.csv", Bulk_Stress);
    writedlm("curvature-par-$par-trial-$trial.csv", Curvature);
    writedlm("index-par-$par-trial-$trial.csv", Index);
    writedlm("filament_speed-par-$par-trial-$trial.csv", Filament_Speed);
    writedlm("motor_speed-par-$par-trial-$trial.csv", Motor_Speed);
    writedlm("motor_pos-par-$par-trial-$trial.csv", Motor_Pos);
    writedlm("motor_angle-par-$par-trial-$trial.csv", Motor_Angle);
    return state, af, mm, xl, Bulk_Stress_Int[end]/times[end], Curvature_Int[end]/times[end], Index_Int[end]/times[end], Filament_Speed_Int[end]/times[end], Motor_Speed_Int[end]/times[end], Angle_ROC_Int[end]/times[end], Motor_Pos_Int[end]/times[end], Motor_Angle_Int[end]/times[end], Bins, Hist
end


text:
"""
Check the correctness of the filtering implementations vs. basic readable math code
"""

using ODEFilters
using LinearAlgebra

@testset "PREDICT" begin
    # Setup
    d = 5
    m = rand(d)
    L_p = LowerTriangular(rand(d,d))
    P = L_p*L_p'

    A = rand(d,d)
    L_Q = LowerTriangular(rand(d,d))
    Q = L_Q*L_Q'

    # PREDICT
    m_p = A*m
    P_p = A*P*A' + Q

    x_curr = Gaussian(m,P)
    x_out = copy(x_curr)

    @testset "predict!" begin
        ODEFilters.predict!(x_out, x_curr, A, Q)
        @test m_p == x_out.μ
        @test P_p == x_out.Σ
    end

    @testset "predict" begin
        x_out = ODEFilters.predict(x_curr, A, Q)
        @test m_p == x_out.μ
        @test P_p == x_out.Σ
    end

    @testset "predict! with PSDMatrix" begin
        x_curr = Gaussian(m, PSDMatrix(L_p))
        x_out = copy(x_curr)
        ODEFilters.predict!(x_out, x_curr, A, PSDMatrix(Q))
        @test m_p == x_out.μ
        @test P_p ≈ x_out.Σ
    end

end

@testset "UPDATE" begin
    # Setup
    d = 5
    m_p = rand(d)
    L_P_p = LowerTriangular(rand(d,d))
    P_p = L_P_p * L_P_p'

    # Measure
    o = 3
    H = rand(o, d)
    L_R = 0.0*rand(o, o)
    R = L_R*L_R'

    z_data = zeros(o)
    z = H*m_p
    S = H*P_p*H' + R

    # UPDATE
    S_inv = inv(S)
    K = P_p * H' * S_inv
    m = m_p + K * (z_data .- z)
    P = P_p - K*S*K'

    x_pred = Gaussian(m_p, P_p)
    x_out = copy(x_pred)
    measurement = Gaussian(z, S)

    @testset "update!" begin
        ODEFilters.update!(x_out, x_pred, measurement, H, R)
        @test m == x_out.μ
        @test P ≈ x_out.Σ
    end

    @testset "update!" begin
        x_out = ODEFilters.update(x_pred, measurement, H, R)
        @test m == x_out.μ
        @test P ≈ x_out.Σ
    end

end

@testset "SMOOTH" begin
    # Setup
    d = 5
    m, m_s = rand(d), rand(d)
    L_P, L_P_s = LowerTriangular(rand(d,d)), LowerTriangular(rand(d,d))
    P, P_s = L_P*L_P', L_P_s*L_P_s'

    A = rand(d,d)
    L_Q = LowerTriangular(rand(d,d))
    Q = L_Q*L_Q'

    # PREDICT first
    m_p = A*m
    P_p = A*P*A' + Q

    # SMOOTH
    G = P * A' * inv(P_p)
    m_smoothed = m + G * (m_s - m_p)
    P_smoothed = P + G * (P_s - P_p) * G'

    x_curr = Gaussian(m, P)
    x_smoothed = Gaussian(m_s, P_s)

    ODEFilters.smooth(x_curr, x_smoothed, A, Q)

    @testset "smooth" begin
        x_out, _ = ODEFilters.smooth(x_curr, x_smoothed, A, Q)
        @test m_smoothed == x_out.μ
        @test P_smoothed ≈ x_out.Σ
    end

end


text:
# """
# 	$(SIGNATURES)

# Write sbatch file to local machine

# IN
# 	filePath
# 		path for the sbatch file
# 	juliaCmd
# 		Julia commands to be written after
# OUT
# 	filePath
# 		path that was actually written
# """
# function write_sbatch_file(filePath :: String, cmdPath :: String,
# 	sbatchCmdV :: Vector{String})

# 	if isempty(dirname(filePath))
# 		filePath = joinpath(sbatch_dir(compName = defaultLocal),  filePath);
# 	end

# 	lineV = sbatch_commands(outPath;
# 		timeStr = timeStr,  nNodes = nNodes,  nCores = nCores,  memPerCpu = memKB);

# 	open(filePath, "w") do io
# 		for line in lineV
# 			write(io, line);
# 		end
# 		write(io, "julia '$cmdPath' \n")
# 	end

# 	return filePath
# end

"""
	$(SIGNATURES)

Write the `#SBATCH` portion of an sbatch file and the "export number of threads" command.

# Arguments
- outPath :: String
	path where output of the run will be written
- timeStr :: String
	e.g. "3-00" for 3 days
- memPerCpu :: Int
	in KB
"""
function sbatch_commands(outPath :: String;
	timeStr :: String = "3-00",  nNodes :: Int = 1,  nCores :: Int = 1,
	memPerCpu :: Int = 16000,
	emailStr = "lhendri@email.unc.edu")

	lineV = [
		"#!/bin/bash \n",
		"#SBATCH -N $nNodes ",
		"#SBATCH -n $nCores ",
		"#SBATCH -t $timeStr ",
		"#SBATCH -p general ",
		"#SBATCH --mem-per-cpu $memPerCpu ",
		"#SBATCH -o '$outPath' ",
		"#SBATCH --mail-type=end ",
		"#SBATCH --mail-user=$emailStr ",
		" ",
		"export JULIA_NUM_THREADS=$nCores "
	];
	return lineV
end

# ----------