text:
"""
    iht(b, x::BEDFile, y, k, g)

If used with a `BEDFile` object `x`, then the temporary arrays `b0`, `Xb`, `Xb0`, and `sortidx` are all initialized as `SharedArray`s of the proper dimensions.
The additional optional arguments are:

- `pids`, a vector of process IDs. Defaults to `procs()`.
- `means`, a vector of SNP means. Defaults to `mean(T, x, shared=true, pids=procs()`.
- `invstds`, a vector of SNP precisions. Defaults to `invstd(x, means, shared=true, pids=procs()`.
"""
function iht{T <: Float}(
    b        :: DenseVector{T},
    x        :: BEDFile,
    y        :: DenseVector{T},
    k        :: Int,
    g        :: DenseVector{T};
    n        :: Int              = length(y),
    p        :: Int              = length(b),
    pids     :: DenseVector{Int} = procs(),
    means    :: DenseVector{T}   = mean(T, x, shared=true, pids=pids),
    invstds  :: DenseVector{T}   = invstd(x, means, shared=true, pids=pids),
    b0       :: DenseVector{T}   = SharedArray(T, p, init = S -> S[localindexes(S)] = b[localindexes(S)], pids=pids),
    Xb       :: DenseVector{T}   = xb(x,b,IDX,k, means=means, invstds=invstds, pids=pids),
    Xb0      :: DenseVector{T}   = SharedArray(T, p, init = S -> S[localindexes(S)] = Xb[localindexes(S)], pids=pids),
    xk       :: DenseMatrix{T}   = zeros(T,n,k),
    xgk      :: DenseVector{T}   = zeros(T,n),
    gk       :: DenseVector{T}   = zeros(T,k),
    IDX      :: BitArray{1}      = falses(p),
    IDX0     :: BitArray{1}      = copy(IDX),
    iter     :: Int              = 1,
    max_step :: Int              = 50,
)

    # which components of beta are nonzero?
    update_indices!(IDX, b, p=p)

    # if current vector is 0,
    # then take largest elements of d as nonzero components for b
    if sum(IDX) == 0
        a = select(g, k, by=abs, rev=true)
#        threshold!(IDX, g, abs(a), n=p)
        IDX[abs(g) .>= abs(a)-2*eps()] = true
    end

    # if support has not changed between iterations,
    # then xk and gk are the same as well
    # avoid extracting and computing them if they have not changed
    # one exception: we should always extract columns on first iteration
    if !isequal(IDX, IDX0) || iter < 2
        decompress_genotypes!(xk, x, IDX, means=means, invstds=invstds)
    end

    # store relevant components of gradient
    fill_perm!(sdata(gk), sdata(g), IDX, k=k, p=p)  # gk = g[IDX]

    # now compute subset of x*g
    BLAS.gemv!('N', one(T), sdata(xk), sdata(gk), zero(T), sdata(xgk))

    # warn if xgk only contains zeros
    all(xgk .== zero(T)) && warn("Entire active set has values equal to 0")

    # compute step size
    mu = sumabs2(sdata(gk)) / sumabs2(sdata(xgk))

    # notify problems with step size
    isfinite(mu) || throw(error("Step size is not finite, is active set all zero?"))
    mu <= eps(typeof(mu))  && warn("Step size $(mu) is below machine precision, algorithm may not converge correctly")

    # take gradient step
    BLAS.axpy!(p, mu, sdata(g), 1, sdata(b), 1)

    # preserve top k components of b
    project_k!(b, k)

    # which indices of new beta are nonzero?
    copy!(IDX0, IDX)
    update_indices!(IDX, b, p=p)

    # must correct for equal entries at kth pivot of b
    # this is a total hack! but matching magnitudes are very rare
    # should not drastically affect performance, esp. with big data
    # hack randomly permutes indices of duplicates and retains one
    if sum(IDX) > k
        a = select(b, k, by=abs, rev=true)          # compute kth pivot
        duples = find(x -> abs(x) .== abs(a), b)    # find duplicates
        c = randperm(length(duples))                # shuffle
        d = duples[c[2:end]]                        # permute, clipping top
        b[d] = zero(T)                             # zero out duplicates
        IDX[d] = false                              # set corresponding indices to false
    end

    # update xb
    xb!(Xb,x,b,IDX,k, means=means, invstds=invstds, pids=pids)

    # calculate omega
    omega_top = sqeuclidean(sdata(b),sdata(b0))
    omega_bot = sqeuclidean(sdata(Xb),sdata(Xb0))

    # backtrack until mu sits below omega and support stabilizes
    mu_step = 0
    while mu*omega_bot > 0.99*omega_top && sum(IDX) != 0 && sum(IDX $ IDX0) != 0 && mu_step < max_step

        # stephalving
        mu /= 2

        # warn if mu falls below machine epsilon
        mu <= eps(typeof(mu)) && warn("Step size equals zero, algorithm may not converge correctly")

        # recompute gradient step
        copy!(b,b0)
        BLAS.axpy!(p, mu, sdata(g), 1, sdata(b), 1)

        # recompute projection onto top k components of b
        project_k!(b, k)

        # which indices of new beta are nonzero?
        update_indices!(IDX, b, p=p)

        # must correct for equal entries at kth pivot of b
        # this is a total hack! but matching magnitudes are very rare
        # should not drastically affect performance, esp. with big data
        # hack randomly permutes indices of duplicates and retains one
        if sum(IDX) > k
            a = select(b, k, by=abs, rev=true)          # compute kth pivot
            duples = find(x -> abs(x) .== abs(a), b)    # find duplicates
            c = randperm(length(duples))                # shuffle
            d = duples[c[2:end]]                        # permute, clipping top
            b[d] = zero(T)                             # zero out duplicates
            IDX[d] = false                              # set corresponding indices to false
        end

        # recompute xb
        xb!(Xb,x,b,IDX,k, means=means, invstds=invstds, pids=pids)

        # calculate omega
        omega_top = sqeuclidean(sdata(b),sdata(b0))
        omega_bot = sqeuclidean(sdata(Xb),sdata(Xb0))

        # increment the counter
        mu_step += 1
    end

    return mu, mu_step
end

"""
    L0_reg(x::BEDFile, y, k)

If used with a `BEDFile` object `x`, then the temporary floating point arrays are all initialized as `SharedArray`s of the proper dimensions.
The additional optional arguments are:

- `pids`, a vector of process IDs. Defaults to `procs()`.
- `means`, a vector of SNP means. Defaults to `mean(T, x, shared=true, pids=procs()`.
- `invstds`, a vector of SNP precisions. Defaults to `invstd(x, means, shared=true, pids=procs()`.
"""
function L0_reg{T <: Float}(
    X        :: BEDFile,
    Y        :: DenseVector{T},
    k        :: Int;
    n        :: Int              = length(Y),
    p        :: Int              = size(X,2),
    pids     :: DenseVector{Int} = procs(),
    Xk       :: DenseMatrix{T}   = SharedArray(T, (n,k), init = S -> S[localindexes(S)] = zero(T), pids=pids),
    means    :: DenseVector{T}   = mean(T,X, shared=true, pids=pids),
    invstds  :: DenseVector{T}   = invstd(X,means, shared=true, pids=pids),
    b        :: DenseVector{T}   = SharedArray(T, p, init = S -> S[localindexes(S)] = zero(T), pids=pids),
    b0       :: DenseVector{T}   = SharedArray(T, p, init = S -> S[localindexes(S)] = zero(T), pids=pids),
    df       :: DenseVector{T}   = SharedArray(T, p, init = S -> S[localindexes(S)] = zero(T), pids=pids),
    r        :: DenseVector{T}   = SharedArray(T, n, init = S -> S[localindexes(S)] = zero(T), pids=pids),
    Xb       :: DenseVector{T}   = SharedArray(T, n, init = S -> S[localindexes(S)] = zero(T), pids=pids),
    Xb0      :: DenseVector{T}   = SharedArray(T, n, init = S -> S[localindexes(S)] = zero(T), pids=pids),
    tempn    :: DenseVector{T}   = SharedArray(T, n, init = S -> S[localindexes(S)] = zero(T), pids=pids),
    gk       :: DenseVector{T}   = SharedArray(T, k, init = S -> S[localindexes(S)] = zero(T), pids=pids),
    support  :: BitArray{1}      = falses(p),
    support0 :: BitArray{1}      = falses(p),
    tol      :: Float            = convert(T, 1e-4),
    max_iter :: Int              = 100,
    max_step :: Int              = 50,
    quiet    :: Bool             = true
)

    # start timer
    tic()

    # first handle errors
    k        >= 0      || throw(ArgumentError("Value of k must be nonnegative!\n"))
    max_iter >= 0      || throw(ArgumentError("Value of max_iter must be nonnegative!\n"))
    max_step >= 0      || throw(ArgumentError("Value of max_step must be nonnegative!\n"))
    tol      >  eps(T) || throw(ArgumentError("Value of global tol must exceed machine precision!\n"))

    # initialize return values
    mm_iter   = 0                 # number of iterations of L0_reg
    mm_time   = zero(T)           # compute time *within* L0_reg
    next_loss = oftype(tol,Inf)   # loss function value

    # initialize floats
    current_obj = oftype(tol,Inf) # tracks previous objective function value
    the_norm    = zero(T)         # norm(b - b0)
    scaled_norm = zero(T)         # the_norm / (norm(b0) + 1)
    mu          = zero(T)         # Landweber step size, 0 < tau < 2/rho_max^2

    # initialize integers
    i       = 0                   # used for iterations in loops
    mu_step = 0                   # counts number of backtracking steps for mu

    # initialize booleans
    converged = false             # scaled_norm < tol?

    # update Xb, r, and gradient
    if sum(support) == 0
        fill!(Xb,zero(T))
        copy!(r,sdata(Y))
    else
        xb!(Xb,X,b,support,k, means=means, invstds=invstds, pids=pids)
        difference!(r, Y, Xb)
    end
    xty!(df, X, r, means=means, invstds=invstds, p=p, pids=pids)

    # formatted output to monitor algorithm progress
    if !quiet
         println("\nBegin MM algorithm\n")
         println("Iter\tHalves\tMu\t\tNorm\t\tObjective")
         println("0\t0\tInf\t\tInf\t\tInf")
    end

    # main loop
    for mm_iter = 1:max_iter

        # notify and break if maximum iterations are reached.
        if mm_iter >= max_iter

            if !quiet
                print_with_color(:red, "MM algorithm has hit maximum iterations $(max_iter)!\n")
                print_with_color(:red, "Current Objective: $(current_obj)\n")
            end

            # send elements below tol to zero
            threshold!(b, tol, n=p)

            # stop timer
            mm_time = toq()

            # these are output variables for function
            # wrap them into a Dict and return
            output = Dict{String, Any}("time" => mm_time, "loss" => next_loss, "iter" => mm_iter, "beta" => b)

            return output
        end

        # save values from previous iterate
        copy!(b0,b)             # b0 = b
        copy!(Xb0,Xb)           # Xb0 = Xb
        current_obj = next_loss

        # now perform IHT step
#        (mu, mu_step) = iht(b,X,Y,k,df, n=n, p=p, max_step=max_step, IDX=support, IDX0=support0, b0=b0, Xb=Xb, Xb0=Xb0, xgk=tempn, xk=Xk, sortidx=indices, gk=idx, means=means, invstds=invstds, iter=mm_iter, pids=pids)
        (mu, mu_step) = iht(b,X,Y,k,df, n=n, p=p, max_step=max_step, IDX=support, IDX0=support0, b0=b0, Xb=Xb, Xb0=Xb0, xgk=tempn, xk=Xk, gk=gk, means=means, invstds=invstds, iter=mm_iter, pids=pids)

        # the IHT kernel gives us an updated x*b
        # use it to recompute residuals and gradient
        difference!(r,Y,Xb)
        xty!(df, X, r, means=means, invstds=invstds, pids=pids)

        # update loss, objective, and gradient
        next_loss = sumabs2(sdata(r)) / 2

        # guard against numerical instabilities
        # ensure that objective is finite
        # if not, throw error
        isnan(next_loss) && throw(error("Objective function is NaN, aborting..."))
        isinf(next_loss) && throw(error("Objective function is Inf32, aborting..."))

        # track convergence
        the_norm    = chebyshev(b,b0)
        scaled_norm = the_norm / ( norm(b0,Inf) + 1)
        converged   = scaled_norm < tol

        # output algorithm progress
        quiet || @printf("%d\t%d\t%3.7f\t%3.7f\t%3.7f\n", mm_iter, mu_step, mu, the_norm, next_loss)

        # check for convergence
        # if converged and in feasible set, then algorithm converged before maximum iteration
        # perform final computations and output return variables
        if converged

            # send elements below tol to zero
            threshold!(b, tol, n=p)

            # stop time
            mm_time = toq()

            if !quiet
                println("\nMM algorithm has converged successfully.")
                println("MM Results:\nIterations: $(mm_iter)")
                println("Final Loss: $(next_loss)")
                println("Total Compute Time: $(mm_time)")
            end

            # these are output variables for function
            # wrap them into a Dict and return
            output = Dict{String, Any}("time" => mm_time, "loss" => next_loss, "iter" => mm_iter, "beta" => b)

            return output
        end

        # algorithm is unconverged at this point.
        # if algorithm is in feasible set, then rho should not be changing
        # check descent property in that case
        # if rho is not changing but objective increases, then abort
        if next_loss > current_obj + tol
            if !quiet
                print_with_color(:red, "\nMM algorithm fails to descend!\n")
                print_with_color(:red, "MM Iteration: $(mm_iter)\n")
                print_with_color(:red, "Current Objective: $(current_obj)\n")
                print_with_color(:red, "Next Objective: $(next_loss)\n")
                print_with_color(:red, "Difference in objectives: $(abs(next_loss - current_obj))\n")
            end
            throw(ErrorException("Descent failure!"))
#           output = Dict{String, Any}("time" => -1.0, "loss" => -1.0, "iter" => -1, "beta" => fill!(b,Inf32))

            return output
        end
    end # end main loop
end # end function

"""
    iht_path(x::BEDFile, y, path)

If used with a `BEDFile` object `x`, then the temporary arrays are all initialized as `SharedArray`s of the proper dimensions.
The additional optional arguments are:

- `pids`, a vector of process IDs. Defaults to `procs()`.
- `means`, a vector of SNP means. Defaults to `mean(T, x, shared=true, pids=procs()`.
- `invstds`, a vector of SNP precisions. Defaults to `invstd(x, means, shared=true, pids=procs()`.
"""
function iht_path{T <: Float}(
    x        :: BEDFile,
    y        :: DenseVector{T},
    path     :: DenseVector{Int};
    pids     :: DenseVector{Int} = procs(),
    means    :: DenseVector{T}   = mean(T,x, shared=true, pids=pids),
    invstds  :: DenseVector{T}   = invstd(x,means, shared=true, pids=pids),
    tol      :: Float            = convert(T, 1e-4),
    max_iter :: Int              = 100,
    max_step :: Int              = 50,
    quiet    :: Bool             = true
)

    # size of problem?
    n = length(y)
    p = size(x,2)

    # how many models will we compute?
    num_models = length(path)

    # preallocate SharedArrays for intermediate steps of algorithm calculations
    b       = SharedArray(T, p, init = S -> S[localindexes(S)] = zero(T), pids=pids)        # for || Y - XB ||_2^2
    b0      = SharedArray(T, p, init = S -> S[localindexes(S)] = zero(T), pids=pids)        # previous iterate beta0
    df      = SharedArray(T, p, init = S -> S[localindexes(S)] = zero(T), pids=pids)        # (negative) gradient
    r       = SharedArray(T, n, init = S -> S[localindexes(S)] = zero(T), pids=pids)        # for || Y - XB ||_2^2
    Xb      = SharedArray(T, n, init = S -> S[localindexes(S)] = zero(T), pids=pids)        # X*beta
    Xb0     = SharedArray(T, n, init = S -> S[localindexes(S)] = zero(T), pids=pids)        # X*beta0
    tempn   = SharedArray(T, n, init = S -> S[localindexes(S)] = zero(T), pids=pids)        # temporary array of n floats

    # index vector for b has more complicated initialization
#    indices = SharedArray(Int, p, init = S -> S[localindexes(S)] = localindexes(S), pids=pids)

    # allocate the BitArrays for indexing in IHT
    # also preallocate matrix to store betas
    support    = falses(p)                # indicates nonzero components of beta
    support0   = copy(support)            # store previous nonzero indicators
    betas      = spzeros(T,p,num_models)  # a matrix to store calculated models

    # compute the path
    @inbounds for i = 1:num_models

        # model size?
        q = path[i]

        # monitor progress
        quiet || print_with_color(:blue, "Computing model size $q.\n\n")

        # these arrays change in size from iteration to iteration
        # we must allocate them for every new model size
        Xk     = zeros(T,n,q)     # store q columns of X
        gk     = zeros(T,q)       # another temporary array of q floats

        # store projection of beta onto largest k nonzeroes in magnitude
        project_k!(b, q)

        # now compute current model
#        output = L0_reg(x,y,q, n=n, p=p, b=b, tol=tol, max_iter=max_iter, max_step=max_step, quiet=quiet, Xk=Xk, r=r, Xb=Xb, Xb=Xb0, b0=b0, df=df, idx=idx, tempn=tempn, indices=indices, support=support, support0=support0, means=means, invstds=invstds, pids=pids)
        output = L0_reg(x,y,q, n=n, p=p, b=b, tol=tol, max_iter=max_iter, max_step=max_step, quiet=quiet, Xk=Xk, r=r, Xb=Xb, Xb=Xb0, b0=b0, df=df, gk=gk, tempn=tempn, support=support, support0=support0, means=means, invstds=invstds, pids=pids)

        # extract and save model
        copy!(sdata(b), output["beta"])

        # ensure that we correctly index the nonzeroes in b
        update_indices!(support, b, p=p)
        fill!(support0, false)

        # put model into sparse matrix of betas
        betas[:,i] = sparsevec(sdata(b))
    end

    # return a sparsified copy of the models
    return betas
end

"""
    one_fold(x::BEDFile, y, path, folds, fold)

If used with a `BEDFile` object `x`, then the additional optional arguments are:

- `pids`, a vector of process IDs. Defaults to `procs()`.
- `means`, a vector of SNP means. Defaults to `mean(T, x, shared=true, pids=procs()`.
- `invstds`, a vector of SNP precisions. Defaults to `invstd(x, means, shared=true, pids=procs()`.
"""
function one_fold{T <: Float}(
    x        :: BEDFile,
    y        :: DenseVector{T},
    path     :: DenseVector{Int},
    folds    :: DenseVector{Int},
    fold     :: Int;
    pids     :: DenseVector{Int} = procs(),
    means    :: DenseVector{T}   = mean(T,x, shared=true, pids=pids),
    invstds  :: DenseVector{T}   = invstd(x,means, shared=true, pids=pids),
    max_iter :: Int              = 100,
    max_step :: Int              = 50,
    quiet    :: Bool             = true
)

    # make vector of indices for folds
    test_idx = folds .== fold
    test_size = sum(test_idx)

    # preallocate vector for output
    myerrors = zeros(T, test_size)

    # train_idx is the vector that indexes the TRAINING set
    train_idx = !test_idx

    # allocate the arrays for the training set
    x_train = x[train_idx,:]
    y_train = y[train_idx]
    Xb      = SharedArray(T, test_size, init = S -> S[localindexes(S)] = zero(T), pids=pids)
    b       = SharedArray(T, size(x,2), init = S -> S[localindexes(S)] = zero(T), pids=pids)
    r       = SharedArray(T, test_size, init = S -> S[localindexes(S)] = zero(T), pids=pids)

    # compute the regularization path on the training set
    betas = iht_path(x_train,y_train,path, max_iter=max_iter, quiet=quiet, max_step=max_step, means=means, invstds=invstds, pids=pids)

    # compute the mean out-of-sample error for the TEST set
    # do this for every computed model in regularization path
    for i = 1:size(betas,2)

        # pull ith model in dense vector format
        b2 = full(vec(betas[:,i]))

        # copy it into SharedArray b
        copy!(b,b2)

        # indices stores Boolean indexes of nonzeroes in b
        update_indices!(indices, b, p=p)

        # compute estimated response Xb with $(path[i]) nonzeroes
        xb!(Xb,x_test,b,indices,path[i],test_idx, means=means, invstds=invstds, pids=pids)

        # compute residuals
        difference!(r,y,Xb)

        # compute out-of-sample error as squared residual averaged over size of test set
        myerrors[i] = sumabs2(r) / test_size / 2
    end

    return myerrors
end

"""
    cv_iht(x::BEDFile, y, path, q)

If used with a `BEDFile` object `x`, then the additional optional arguments are:

- `pids`, a vector of process IDs. Defaults to `procs()`.
- `means`, a vector of SNP means. Defaults to `mean(T, x, shared=true, pids=procs()`.
- `invstds`, a vector of SNP precisions. Defaults to `invstd(x, means, shared=true, pids=procs()`.
"""
function cv_iht{T <: Float}(
    x             :: BEDFile,
    y             :: DenseVector{T},
    path          :: DenseVector{Int},
    q             :: Int;
    folds         :: DenseVector{Int} = cv_get_folds(sdata(y),q),
    pids          :: DenseVector{Int} = procs(),
    means         :: DenseVector{T}   = mean(T,x, shared=true, pids=pids),
    invstds       :: DenseVector{T}   = invstd(x,means, shared=true, pids=pids),
    tol           :: Float            = convert(T, 1e-4),
    n             :: Int              = length(y),
    p             :: Int              = size(x,2),
    max_iter      :: Int              = 100,
    max_step      :: Int              = 50,
    quiet         :: Bool             = true,
    compute_model :: Bool             = false
)

    # how many elements are in the path?
    num_models = length(path)

    # preallocate vectors used in xval
    errors  = zeros(T, num_models)    # vector to save mean squared errors

    # want to compute a path for each fold
    # the folds are computed asynchronously
    # the @sync macro ensures that we wait for all of them to finish before proceeding
    @sync for i = 1:q
        # one_fold returns a vector of out-of-sample errors (MSE for linear regression)
        # @fetch(one_fold(...)) spawns one_fold() on a CPU and returns the MSE
        errors[i] = @fetch(one_fold(x, y, path, folds, i, max_iter=max_iter, max_step=max_step, quiet=quiet, means=means, invstds=invstds, pids=pids))
    end

    # average the mses
    errors ./= q

    # what is the best model size?
    k = convert(Int, floor(mean(path[errors .== minimum(errors)])))

    # print results
    if !quiet
        println("\n\nCrossvalidation Results:")
        println("k\tMSE")
        @inbounds for i = 1:num_models
            println(path[i], "\t", errors[i])
        end
        println("\nThe lowest MSE is achieved at k = ", k)
    end

    # recompute ideal model
    if compute_model

        # initialize model
        b = SharedArray(T, p, init = S -> S[localindexes(S)] = zero(T), pids=pids)

        # first use L0_reg to extract model
        output = L0_reg(x,y,k, max_iter=max_iter, max_step=max_step, quiet=quiet, tol=tol, b=b, pids=pids)

        # which components of beta are nonzero?
        inferred_model = output["beta"] .!= zero(T)
        bidx = find( x -> x .!= zero(T), b)

        # allocate the submatrix of x corresponding to the inferred model
        x_inferred = zeros(T,n,sum(inferred_model))
        decompress_genotypes!(x_inferred,x)

        # now estimate b with the ordinary least squares estimator b = inv(x'x)x'y
        Xty = BLAS.gemv('T', one(T), x_inferred, y)
        xtx = BLAS.gemm('T', 'N', one(T), x_inferred, x_inferred)
        b   = xtx \ Xty
        return errors, b, bidx
    end
    return errors
end


text:
module TestNormalForm
include("preamble.jl")

using Bifurcations.Continuations: as, ContinuationSolution, sweeps_as_vectors
using Bifurcations.Examples.Pitchfork
using Bifurcations.Examples.Transcritical

@testset "$normal" for normal in [Pitchfork, Transcritical]
    solver = init(normal.prob; h0=0.3)
    solve!(solver)
    sol = as(solver.sol, ContinuationSolution)
    rtol = solver.opts.rtol
    atol = solver.opts.atol
    for (us, ps) in zip(sweeps_as_vectors(sol, 1),
                        sweeps_as_vectors(sol, 0))
        analytic = hcat(normal.closest_analytic.(us, ps)...)'
        numeric = hcat(us, ps)
        @test numeric ≈ analytic  rtol=2rtol atol=atol
        # TODO: Don't double rtol; I think rtol=rtol is not working
        # because the way closest_analytic chooses the "closest" point
        # is actually not the closest.
    end

    # One of the first two sweeps has to detect a simple bifurcation
    num_bifurcations = [length(sweep.simple_bifurcation)
                        for sweep in sol.sweeps]
    @test Set(num_bifurcations[1:2]) == Set([0, 1])
    _, i = findmax(num_bifurcations)
    sweep = sol.sweeps[i]
    b = sweep.simple_bifurcation[1]
    @test b.u0[end] <= 0 <= b.u1[end]
end

end  # module

ext: jl


text:
@static if Sys.isunix()

const PROT_READ     = Cint(1)
const PROT_WRITE    = Cint(2)
const MAP_SHARED    = Cint(1)
const MAP_PRIVATE   = Cint(2)
const MAP_ANONYMOUS = Cint(Sys.isbsd() ? 0x1000 : 0x20)

const MPOL_DEFAULT     = Cint(0)
const MPOL_PREFERRED   = Cint(1)
const MPOL_BIND        = Cint(2)
const MPOL_INTERLEAVE  = Cint(3)
const MPOL_LOCAL       = Cint(4)
const MPOL_MF_STRICT   = Cint(1)
const MPOL_MF_MOVE     = Cint(2)
const MPOL_MF_MOVE_ALL = Cint(4)

const MPOL_F_NODE      = Cint(1)
const MPOL_F_ADDR      = Cint(2)

function get_mempolicy(ptr::Ptr)
  mode = Ref(~Culong(0))
  nodemask = Ref{Culong}(1 << nodeid)
  ret = ccall((:get_mempolicy, libnuma), Cint, (Ptr{Cvoid}, Ptr{Culong}, Culong, Culong, Culong),
              mode, nodemask, sizeof(Culong)*8, ptr, MPOL_F_ADDR)
  systemerror("get_mempolicy failed", ret == -1)

  mode[], nodemask[]
end

function mcurrent_node(ptr::Ptr{T}) where {T}
  page_size = page_size_bytes()
  ptr = reinterpret(Ptr{T}, reinterpret(UInt, ptr) & ~(page_size-1))
  pages = [ptr]

  status = Ref{Cint}(-1)
  ret = ccall((:move_pages, libnuma), Clong, (Cint, Culong, Ptr{Ptr{Cvoid}}, Ptr{Cint}, Ptr{Cint}, Cint),
              0, 1, pages, C_NULL, status, 0)
  systemerror("move_pages failed", ret == -1)
  status
end

function mmap_mbind(::Type{Array{T,N}}, dims::NTuple{N, Integer}, nodeid::Int64) where {T,N}
    size = prod(dims) * sizeof(T)
    size >= 0 || throw(ArgumentError("size should be >= 0, but $size"))

    prot = PROT_READ | PROT_WRITE
    flags = MAP_PRIVATE | MAP_ANONYMOUS
    ptr = ccall(:jl_mmap, Ptr{Cvoid}, (Ptr{Cvoid}, Csize_t, Cint, Cint, Cint, Int64),
            C_NULL, size, prot, flags, -1, 0)
    systemerror("memory mapping failed", reinterpret(Int, ptr) == -1)

    try
        nodemask = Ref{Culong}(1 << nodeid)
        ret = ccall((:mbind, libnuma), Clong, (Ptr{Cvoid}, Culong, Cint, Ptr{Culong}, Culong, Cuint),
                ptr, size, MPOL_BIND, nodemask, sizeof(Culong)*8, MPOL_MF_MOVE)
        systemerror("mbind failed", ret == -1)
    catch
        systemerror("munmap",  ccall(:munmap, Cint, (Ptr{Cvoid}, Int), ptr, size) != 0)
        rethrow()
    end

    A = unsafe_wrap(Array, convert(Ptr{T}, ptr), dims)
    finalizer(A) do x
        @static if Sys.isunix()
            systemerror("munmap",  ccall(:munmap, Cint, (Ptr{Cvoid}, Int), ptr, size) != 0)
        end
    end

    A
end

const MBytes = 1024*1024
const KBytes = 1024

page_size_bytes() = Int(ccall(:jl_getpagesize, Clong, ()))
function huge_page_size_in_bytes()
    for line in eachline("/proc/meminfo")
        if startswith(line, "Hugepagesize")
            _, count, unit = split(line)
            size = parse(Int64, count)
            return if unit == "kB"
                size * KBytes
            else
                error("unknown unit $unit")
            end
        end
    end
end

function parse_node(descr::AbstractString)
    node, pages = split(descr, '=')
    parse(Int64, node[2:end]), parse(Int64, pages)
end

function numastat(pid::Int32=getpid())
    categories = ["huge", "heap", "stack", "N"]

    default_page_size = page_size_bytes()
    huge_page_size = huge_page_size_in_bytes()

    NT = NamedTuple{(:huge, :heap, :stack, :private, :total), NTuple{5, Float64}}

    numa_maps = "/proc/$pid/numa_maps"
    open(numa_maps) do io
        table = Dict{Int64, NT}()

        for line in eachline(io)
            category = 4
            parts = split(line)

            for p in parts
                if category == 4
                    idx = findfirst(cat -> startswith(p, cat), categories)
                    idx !== nothing && (category = idx)
                end

                # per-node page quantity
                if p[1] == 'N'
                    node, pages = parse_node(p)

                    multiplier = category == 1 ? huge_page_size : default_page_size
                    mbytes = pages * multiplier / MBytes

                    vals = get(table, node_num, NT((0.0, 0.0, 0.0, 0.0, 0.0)))
                    table[node_num] = NT((vals[1] + (category == 1 ? mbytes : 0.),
                                         vals[2] + (category == 2 ? mbytes : 0.),
                                         vals[3] + (category == 3 ? mbytes : 0.),
                                         vals[4] + (category == 4 ? mbytes : 0.),
                                         vals[5] + mbytes))
                end
            end
        end

        table
    end
end

function find_numa_mapping(ptr::UInt64, pid::Int32=getpid())
    default_page_size = page_size_bytes()
    huge_page_size = huge_page_size_in_bytes()

    numa_maps = "/proc/$pid/numa_maps"
    line = open(numa_maps) do io
        prev = nothing
        for line in eachline(io)
            parts = split(line; limit=2)
            address = parse(UInt64, parts[1], base=16)
            address > ptr && break
            prev = line
        end
        return prev
    end

    if line === nothing
        error("mapping for $ptr not found")
    end

    parts = split(line)
    offset = ptr - parse(UInt64, parts[1], base=16)
    nodepages = map(parse_node, filter(startswith("N"), parts))

    multiplier = if "huge" in parts
        huge_page_size
    else
        default_page_size
    end

    pages = sum(last, nodepages)
    bytes = pages * multiplier
    @assert offset < bytes

    Dict(((n, p*multiplier/MBytes) for (n,p) in nodepages))
end

find_numa_mapping(ptr::Ptr) = find_numa_mapping(convert(UInt64, ptr))
find_numa_mapping(a::AbstractArray) = find_numa_mapping(pointer(a))

function numaprofile(f::Function; gcfirst::Bool=true, warmup::Int64=0)
    node_stores = event"node-stores"
    node_loads = event"node-loads"
    node_store_misses = event"node-store-misses"
    node_load_misses = event"node-load-misses"
    profile(f, Event[node_stores, node_loads, node_store_misses, node_load_misses], gcfirst=gcfirst, warmup=warmup)
end

"""
    numaprofile(ex, args...)

Convience macro for numa profiling an expression.
"""
macro numaprofile(ex, args...)
    params = collect(args)
    for ex in params
        if isa(ex, Expr) && ex.head == :(=)
            ex.head = :kw
        end
    end
    quote
        numaprofile(() -> $(esc(ex)), $(params...))
    end
end

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# The epochs used for date rounding are based ISO 8601's "year zero" notation
const DATEEPOCH = value(Date(0))
const DATETIMEEPOCH = value(DateTime(0))

# According to ISO 8601, the first day of the first week of year 0000 is 0000-01-03
const WEEKEPOCH = value(Date(0, 1, 3))

"""
    epochdays2date(days) -> Date

Takes the number of days since the rounding epoch (`0000-01-01T00:00:00`) and returns the
corresponding `Date`.
"""
epochdays2date(i) = Date(UTD(DATEEPOCH + Int64(i)))

"""
    epochms2datetime(milliseconds) -> DateTime

Takes the number of milliseconds since the rounding epoch (`0000-01-01T00:00:00`) and
returns the corresponding `DateTime`.
"""
epochms2datetime(i) = DateTime(UTM(DATETIMEEPOCH + Int64(i)))

"""
    date2epochdays(dt::Date) -> Int64

Takes the given `Date` and returns the number of days since the rounding epoch
(`0000-01-01T00:00:00`) as an `Int64`.
"""
date2epochdays(dt::Date) = value(dt) - DATEEPOCH

"""
    datetime2epochms(dt::DateTime) -> Int64

Takes the given `DateTime` and returns the number of milliseconds since the rounding epoch
(`0000-01-01T00:00:00`) as an `Int64`.
"""
datetime2epochms(dt::DateTime) = value(dt) - DATETIMEEPOCH

function Base.floor(dt::Date, p::Year)
    value(p) < 1 && throw(DomainError())
    years = year(dt)
    return Date(years - mod(years, value(p)))
end

function Base.floor(dt::Date, p::Month)
    value(p) < 1 && throw(DomainError())
    y, m = yearmonth(dt)
    months_since_epoch = y * 12 + m - 1
    month_offset = months_since_epoch - mod(months_since_epoch, value(p))
    target_month = mod(month_offset, 12) + 1
    target_year = div(month_offset, 12) - (month_offset < 0 && target_month != 1)
    return Date(target_year, target_month)
end

function Base.floor(dt::Date, p::Week)
    value(p) < 1 && throw(DomainError())
    days = value(dt) - WEEKEPOCH
    days = days - mod(days, value(Day(p)))
    return Date(UTD(WEEKEPOCH + Int64(days)))
end

function Base.floor(dt::Date, p::Day)
    value(p) < 1 && throw(DomainError())
    days = date2epochdays(dt)
    return epochdays2date(days - mod(days, value(p)))
end

Base.floor(dt::DateTime, p::DatePeriod) = DateTime(Base.floor(Date(dt), p))

function Base.floor(dt::DateTime, p::TimePeriod)
    value(p) < 1 && throw(DomainError())
    milliseconds = datetime2epochms(dt)
    return epochms2datetime(milliseconds - mod(milliseconds, value(Millisecond(p))))
end

"""
    floor(dt::TimeType, p::Period) -> TimeType

Returns the nearest `Date` or `DateTime` less than or equal to `dt` at resolution `p`.

For convenience, `p` may be a type instead of a value: `floor(dt, Dates.Hour)` is a shortcut
for `floor(dt, Dates.Hour(1))`.

```jldoctest
julia> floor(Date(1985, 8, 16), Dates.Month)
1985-08-01

julia> floor(DateTime(2013, 2, 13, 0, 31, 20), Dates.Minute(15))
2013-02-13T00:30:00

julia> floor(DateTime(2016, 8, 6, 12, 0, 0), Dates.Day)
2016-08-06T00:00:00
```
"""
Base.floor(::Dates.TimeType, ::Dates.Period)

"""
    ceil(dt::TimeType, p::Period) -> TimeType

Returns the nearest `Date` or `DateTime` greater than or equal to `dt` at resolution `p`.

For convenience, `p` may be a type instead of a value: `ceil(dt, Dates.Hour)` is a shortcut
for `ceil(dt, Dates.Hour(1))`.

```jldoctest
julia> ceil(Date(1985, 8, 16), Dates.Month)
1985-09-01

julia> ceil(DateTime(2013, 2, 13, 0, 31, 20), Dates.Minute(15))
2013-02-13T00:45:00

julia> ceil(DateTime(2016, 8, 6, 12, 0, 0), Dates.Day)
2016-08-07T00:00:00
```
"""
function Base.ceil(dt::TimeType, p::Period)
    f = floor(dt, p)
    return (dt == f) ? f : f + p
end

"""
    floorceil(dt::TimeType, p::Period) -> (TimeType, TimeType)

Simultaneously return the `floor` and `ceil` of a `Date` or `DateTime` at resolution `p`.
More efficient than calling both `floor` and `ceil` individually.
"""
function floorceil(dt::TimeType, p::Period)
    f = floor(dt, p)
    return f, (dt == f) ? f : f + p
end

"""
    round(dt::TimeType, p::Period, [r::RoundingMode]) -> TimeType

Returns the `Date` or `DateTime` nearest to `dt` at resolution `p`. By default
(`RoundNearestTiesUp`), ties (e.g., rounding 9:30 to the nearest hour) will be rounded up.

For convenience, `p` may be a type instead of a value: `round(dt, Dates.Hour)` is a shortcut
for `round(dt, Dates.Hour(1))`.

```jldoctest
julia> round(Date(1985, 8, 16), Dates.Month)
1985-08-01

julia> round(DateTime(2013, 2, 13, 0, 31, 20), Dates.Minute(15))
2013-02-13T00:30:00

julia> round(DateTime(2016, 8, 6, 12, 0, 0), Dates.Day)
2016-08-07T00:00:00
```

Valid rounding modes for `round(::TimeType, ::Period, ::RoundingMode)` are
`RoundNearestTiesUp` (default), `RoundDown` (`floor`), and `RoundUp` (`ceil`).
"""
function Base.round(dt::TimeType, p::Period, r::RoundingMode{:NearestTiesUp})
    f, c = floorceil(dt, p)
    return (dt - f) < (c - dt) ? f : c
end

Base.round(dt::TimeType, p::Period, r::RoundingMode{:Down}) = Base.floor(dt, p)
Base.round(dt::TimeType, p::Period, r::RoundingMode{:Up}) = Base.ceil(dt, p)

# No implementation of other `RoundingMode`s: rounding to nearest "even" is skipped because
# "even" is not defined for Period; rounding toward/away from zero is skipped because ISO
# 8601's year 0000 is not really "zero".
Base.round(::TimeType, ::Period, ::RoundingMode) = throw(DomainError())

# Default to RoundNearestTiesUp.
Base.round(dt::TimeType, p::Period) = Base.round(dt, p, RoundNearestTiesUp)

# Make rounding functions callable using Period types in addition to values.
Base.floor(dt::TimeType, p::Type{<:Period}) = Base.floor(dt, p(1))
Base.ceil(dt::TimeType, p::Type{<:Period}) = Base.ceil(dt, p(1))

function Base.round(dt::TimeType, p::Type{<:Period}, r::RoundingMode=RoundNearestTiesUp)
    return Base.round(dt, p(1), r)
end

max_stars_repo_name: odow/julia

max_issues_repo_name: odow/julia

max_forks_repo_name: odow/julia


text:
#!/bin/bash
# -*- mode: julia -*-
#=
exec julia --project=@. --compile=min -O0 --startup-file=no -e "include(popfirst!(ARGS))" "${BASH_SOURCE[0]}" "$@"
=#

using Raw2CSV

# Basically:
#
# `gather_csv_from_folder`: takes one folder and a list of fields to extract,
# extract all fields from all files inside the folder and return an String
# with a CSV.
#
# Alternative:
# `gather_csv_from_folders`: the same as above but take a list of folders
#  (i.e., consider all files inside all folders).
#
# `key_equals` takes two arguments: the name of the field to extract (this
# is both the name in the output file and the name of the column in the CSV),
# and the default value (if the field is not found). If there is no default
# value, then you need to pass `NoDefault{TypeOfTheField}()` so the code knows
# which is the type of the field (to better print the CSV). Valid types for
# fields are: `Int` (for integer data), `Float64` (for floating point data),
# and `String` (for categories and everything else that may be parsed later).
#
# If you wanna to have a column name different from the name of the field
# inside the files then, instead of `key_equals` just use:
#
# ("column_name", key_equals_extractor("field_name_in_log", <default>))
#

csv = gather_csv_from_folder("./data", [
	# The code unfortunately does not automatically save the name of the log
	# file in a column (so you know exactly from where each line came), you
	# need to save the name in the log and extract it.
	key_equals("output_filename", NoDefault{String}()),
	# Here we extract a field that only exists in some files and say that,
	# if the field does not appear in the file, the default is 0 (the code
	# knows the field has integers because 0 is an integer).
	key_equals("only_in_1_and_2", 0),
	# If the above was instead:
	#
	# key_equals("only_in_1_and_2", NoDefault{Int}()),
	#
	# then you would have an error message because one of the files
	# does not have this field. If it was:
	#
	# key_equals("only_in_1_and_2", 0.0),
	#
	# then the numbers would be saved as floating point values in the CSV.
	("custom_column_name", key_equals_extractor("in_all", NoDefault{String}())),
	# float_data field exists in all files, but there is no problem in
	# having a default if you do not want to trigger errors in case it
	# could not be present.
	key_equals("float_data", 0.0)
])
print(csv)


text:
"""
    eom_var1coupled(x, p, n) -> SVector{3}

Equations of motion a system consisting of two mutually coupled first order
autoregressive processes. For `c = 0`, only Y -> X. For `c > 0`, there is
bidirectional coupling. Each subcomponent is dynamically affected by
an independent Gaussian random processes (with standard deviations `σx` and
`σy`) at each time step.

"""
function eom_var1coupled(x, p, n)
    c, σx, σy = (p...,)
    x, y = (x...,)

    ϵx = rand(Normal(0, σx))
    ϵy = rand(Normal(0, σy))

    dx = 0.5*x + 0.2*y + ϵx
    dy = c*x + 0.7*y + ϵy
    return SVector{2}(dx, dy)
end

function var1coupled(uᵢ, c, σx, σy)
    p = [c, σx, σy]
    logistic_system = DiscreteDynamicalSystem(eom_var1coupled, uᵢ, p)

    return logistic_system
end

"""
    var1coupled(;uᵢ = rand(2), c = 0, σx = sqrt(0.1), σy = sqrt(0.1))

Initialise a system consisting of two mutually coupled first order
autoregressive processes. For `c = 0`, only Y -> X. For `c > 0`, there is
bidirectional coupling. Each subcomponent is dynamically affected by
an independent Gaussian random processes (with standard deviations `σx` and
`σy`) at each time step.
"""
var1coupled(;uᵢ = rand(2), c = 0, σx = sqrt(0.1), σy = sqrt(0.1)) =
var1coupled(uᵢ, c, σx, σy)


text:
# Copyright 2017-2019 Eric S. Tellez
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using StatsBase
export accuracy, precision_recall, precision, recall, f1, scores

"""
It computes the recall between the gold dataset and the list of predictions `predict`

It applies the desired weighting scheme for binary and multiclass problems
- `:macro` performs a uniform weigth to each class
- `:weigthed` the weight of each class is proportional to its population in gold
- `:micro` returns the global recall, without distinguishing among classes
"""
function recall(gold, predict; weight=:macro)::Float64
    P = precision_recall(gold, predict)
    if weight == :macro
        mean(x -> x.recall, values(P.per_class))
    elseif weight == :weighted
        mean(x -> x.recall * x.population / length(gold), values(P.per_class))
    elseif :micro
        P.recall
    else
        throw(Exception("Unknown weighting method $weight"))
    end
end

"""
It computes the precision between the gold dataset and the list of predictions `predict`

It applies the desired weighting scheme for binary and multiclass problems
- `:macro` performs a uniform weigth to each class
- `:weigthed` the weight of each class is proportional to its population in gold
- `:micro` returns the global precision, without distinguishing among classes
"""
function precision(gold, predict; weight=:macro)::Float64
    P = precision_recall(gold, predict)
    if weight == :macro
        mean(x -> x.precision, values(P.per_class))
    elseif weight == :weighted
        mean(x -> x.precision * x.population / length(gold), values(P.per_class))
    elseif weight == :micro
        P.precision
    else
        throw(Exception("Unknown weighting method $weight"))
    end
end

"""
It computes the F1 score between the gold dataset and the list of predictions `predict`

It applies the desired weighting scheme for binary and multiclass problems
- `:macro` performs a uniform weigth to each class
- `:weigthed` the weight of each class is proportional to its population in gold
- `:micro` returns the global F1, without distinguishing among classes
"""
function f1(gold, predict; weight=:macro)::Float64
    P = precision_recall(gold, predict)
    function _f1(x)::Float64
        d = x.precision + x.recall
        if d == 0
            0.0
        else
            2 * x.precision * x.recall / d
        end
    end
    if weight == :macro
        mean(x -> _f1(x), values(P.per_class))
    elseif weight == :weighted
        mean(x -> _f1(x) * x.population / length(gold), values(P.per_class))
    elseif weight == :micro
        _f1(P)
    else
        throw(Exception("Unknown weighting method $weight"))
    end
end

"""
Computes precision, recall, and f1 scores, for global and per-class granularity
"""
function scores(gold, predicted)
    P = precision_recall(gold, predicted)
    m = Dict(
        :micro_f1 => 2 * P.precision * P.recall / (P.precision + P.recall),
		:precision => P.precision,
		:recall => P.recall,
		:class_f1 => Dict(),
		:class_precision => Dict(),
		:class_recall => Dict()
    )

    for (k, v) in P.per_class
        m[:class_f1][k] = 2 * v.precision * v.recall / (v.precision + v.recall)
		m[:class_precision][k] = v.precision
		m[:class_recall][k] = v.recall
    end

    m[:macro_recall] = mean(values(m[:class_recall]))
    m[:macro_f1] = mean(values(m[:class_f1]))
    m[:accuracy] = accuracy(gold, predicted)
    m
end

"""
It computes the global and per-class precision and recall values between the gold standard
and the predicted set
"""
function precision_recall(gold, predicted)
    labels = unique(gold)
    M = Dict{typeof(labels[1]), NamedTuple}()
    tp_ = 0
    tn_ = 0
    fn_ = 0
    fp_ = 0

    for label in labels
        lgold = label .== gold
        lpred = label .== predicted

        tp = 0
        tn = 0
        fn = 0
        fp = 0
        for i in 1:length(lgold)
            if lgold[i] == lpred[i]
                if lgold[i]
                    tp += 1
                else
                    tn += 1
                end
            else
                if lgold[i]
                    fn += 1
                else
                    fp += 1
                end
            end
        end

        tp_ += tp
        tn_ += tn
        fn_ += fn
        fp_ += fp
        precision = tp / (tp + fp)

        if precision <= eps(typeof(precision))
            @info "precision is zero for label '$label'; #classes=$(length(labels)) "
        end
        M[label] = (precision=precision, recall=tp / (tp + fn), population=sum(lgold) |> Int)
    end

    (precision=tp_ / (tp_ + fp_), recall=tp_ / (tp_ + fn_), per_class=M)
end

"""
It computes the accuracy score between the gold and the predicted sets
"""
function accuracy(gold, predicted)
    #  mean(gold .== predicted)
    c = 0
    for i in 1:length(gold)
        c += (gold[i] == predicted[i])
    end

    c / length(gold)
end

######### Regression ########

export pearson, spearman, isqerror
"""
Pearson correlation score
"""
function pearson(X::AbstractVector{F}, Y::AbstractVector{F}) where {F <: AbstractFloat}
    X̄ = mean(X)
    Ȳ = mean(Y)
    n = length(X)
    sumXY = 0.0
    sumX2 = 0.0
    sumY2 = 0.0
    for i in 1:n
        x, y = X[i], Y[i]
        sumXY += x * y
        sumX2 += x * x
        sumY2 += y * y
    end
    num = sumXY - n * X̄ * Ȳ
    den = sqrt(sumX2 - n * X̄^2) * sqrt(sumY2 - n * Ȳ^2)
    num / den
end

"""
Spearman rank correleation score
"""
function spearman(X::AbstractVector{F}, Y::AbstractVector{F}) where {F <: AbstractFloat}
    n = length(X)
    x = invperm(sortperm(X))
    y = invperm(sortperm(Y))
    d = x - y
    1 - 6 * sum(d.^2) / (n * (n^2 - 1))
end

"""
Negative squared error (to be used for maximizing algorithms)
"""
function isqerror(X::AbstractVector{F}, Y::AbstractVector{F}) where {F <: AbstractFloat}
    n = length(X)
    d = 0.0

    @inbounds for i in 1:n
        d += (X[i] - Y[i])^2
    end

    -d
end


text:
using DataFrames, MixedModels

sleepstudy =  MixedModels.dataset(:sleepstudy);

using VegaLite

# this works nicely, but only if I don't also do a transformation
sleepstudy |> @vlplot(
    #transform=[{regression=:reaction, on=:days}],
    wrap=:subj,
    columns=3,
    x=:days,y=:reaction,
    mark=:point,
)

using Gadfly

## attempt 1

import Base.Iterators: partition

function partition_facets(facets, n)
    parts = Vector.(collect(partition(facets, n)))

    len = length(first(parts))
    empty = Gadfly.context()

    while length(last(parts)) < length(first(parts))
        push!(last(parts), empty)
    end
    parts
end

facets = map(collect(groupby(sleepstudy, :subj))) do df
    pp = plot(df,
        x=:days,
        y=:reaction,
        Guide.title(first(df.subj)),
        Geom.point,
        layer(Stat.smooth(method=:lm), Geom.line)
    )
    pp
end

ncols = 4

rows = map(partition_facets(facets, ncols)) do rr
    hstack(rr);
end

p = vstack(rows...)

##

## attempt 2
ncols = 4
groups = DataFrame(subj = unique(sleepstudy.subj))
groups.col = rem.(0:nrow(groups)-1, ncols)
groups.row = div.(0:nrow(groups)-1, ncols)

dat = leftjoin(sleepstudy, groups, on=:subj)

plot(dat,
    x=:days,
    y=:reaction,
    xgroup=:row,
    ygroup=:col,
    Geom.subplot_grid(Geom.point,
                     layer(Stat.smooth(method=:lm), Geom.line),
    ),
    Guide.xlabel("Days"),
    Guide.ylabel("Reaction time (ms)"),
    Scale.xgroup(),
    Scale.ygroup(),
)

##

## this doesn't work
using AlgebraOfGraphics, AbstractPlotting, CairoMakie
vars = AlgebraOfGraphics.style(:days, :reaction)
grp = AlgebraOfGraphics.style(layout_x = :subj)
geom = spec(Scatter) + AlgebraOfGraphics.linear
data(sleepstudy) * vars * grp * geom |> draw

max_stars_repo_name: palday/small-multiples

max_issues_repo_name: palday/small-multiples

max_forks_repo_name: palday/small-multiples


text:
#=
    Interpolate using weighted essentially non-oscillatory (WENO) techniques.
=#

export WENO4Interpolation

"""
    WENOInterpolationType

Abstract class for all types of WENO interpolation.
"""
abstract type WENOInterpolationType <: InterpolationType end

"""
    WENO4Interpolation

Fourth-order WENO interpolation. Based on Janett et al. (2019),
"A novel fourth-order WENO interpolation technique. A possible new tool designed for radiative transfer",
https://ui.adsabs.harvard.edu/abs/2019A&A...624A.104J
"""
struct WENO4Interpolation{T} <: WENOInterpolationType where T <: Number
    nodes::Vector{T}
    α2::Vector{T}  # coefficients for the α2 and α3 weights
    α3::Vector{T}  # coefficients for the α2 and α3 weights
    q2::Array{T, 2}  # coefficients for the quadratic Lagrange polynomials q2
    q3::Array{T, 2}  # coefficients for the quadratic Lagrange polynomials q3
end

function Base.:(==)(o1::WENO4Interpolation, o2::WENO4Interpolation)
    o1.nodes == o2.nodes &&
    o1.α2 == o2.α2 &&
    o1.α3 == o2.α3 &&
    o1.q2 == o2.q2 &&
    o1.q3 == o2.q3
end

lbounds(itp::WENO4Interpolation) = (first(itp.nodes),)
ubounds(itp::WENO4Interpolation) = (last(itp.nodes),)

function interpolate(nodes::AbstractVector{<:Number}, A::AbstractVector{<:Number})
    n = length(nodes)
    α2 = Vector{eltype(nodes)}(undef, n)
    α3 = Vector{eltype(nodes)}(undef, n)
    q2 = Array{eltype(nodes)}(undef, 3, n)
    q3 = Array{eltype(nodes)}(undef, 3, n)

    for i ∈ 1:n
        xim, xi, xip, xipp = get_stencil_nodes(i, nodes)
        yim, yi, yip, yipp = get_stencil_values(i, A)
        q2[:, i], q3[:, i] = weno4_q(xim, xi, xip, xipp, yim, yi, yip, yipp)
        α2[i], α3[i] = weno4_α(xim, xi, xip, xipp, yim, yi, yip, yipp)
    end
    WENO4Interpolation(nodes, α2, α3, q2, q3)
end

function (itp::WENO4Interpolation)(x::Number)
    @boundscheck (checkbounds(Bool, itp, x) || Base.throw_boundserror(itp, (x,)))
    n = length(itp.nodes)
    i = binary_search(itp.nodes, x)
    xim, xi, xip, xipp = get_stencil_nodes(i, itp.nodes)
    Δim = x - xim
    Δi = x - xi
    Δip = x - xip
    Δipp = x - xipp
    q2 = itp.q2[1, i] * Δi * Δip + itp.q2[2, i] * Δim * Δip +  itp.q2[3, i] * Δim * Δi
    q3 = itp.q3[1, i] * Δip * Δipp + itp.q3[2, i] * Δi * Δipp +  itp.q3[3, i] * Δi * Δip
    if i == 1
        result = q3
    elseif i == n - 1
        result = q2
    else
        α2 = Δipp * itp.α2[i]
        α3 = Δip * itp.α3[i]
        ω2 = α2 / (α2 + α3)
        ω3 = α3 / (α2 + α3)
        result = ω2 * q2 + ω3 * q3
    end
    return result
end

function get_stencil_nodes(index::Int, nodes::AbstractVector{<:Number})
    n = length(nodes)
    if index == n
        index -= 1
    end
    xi = nodes[index]
    xip = nodes[index+1]
    if index == 1
        xim = zero(eltype(nodes))
        xipp = nodes[index+2]
    elseif index == n - 1
        xim = nodes[index-1]
        xipp = zero(eltype(nodes))
    else
        xim = nodes[index-1]
        xipp = nodes[index+2]
    end
    return (xim, xi, xip, xipp)
end

function get_stencil_values(index::Int, A::AbstractVector{<:Number})
    n = length(A)
    if index == n
        index -= 1
    end
    yi = A[index]
    yip = A[index+1]
    if index == 1
        yim = zero(eltype(A))
        yipp = A[index+2]
    elseif index == n - 1
        yim = A[index-1]
        yipp = zero(eltype(A))
    else
        yim= A[index-1]
        yipp = A[index+2]
    end
    return (yim, yi, yip, yipp)
end

"""
    binary_search(a::AbstractVector{<:Number}, x)

Return the index of the last value in `a` less than or equal to `x`.
`a` is assumed to be sorted.
"""
function binary_search(a::AbstractVector{<:Number}, x)
    n = length(a)
    ileft = 2
    i = n - 1
    while ileft <= i
        mid = (ileft + i) >> 1
        if @inbounds x < a[mid]
            i = mid - 1
        else
            ileft = mid + 1
        end
    end
    return i
end

"""
Compute the smoothness indicators β2 and β3 for 4th order WENO interpolation.
Arguments are the x and y values for the 4-point stencil.
"""
function weno4_α(xim, xi, xip, xipp, yim, yi, yip, yipp)
    ε = 1.0f-6
    him = xi - xim
    hi = xip - xi
    hip = xipp - xip
    H = him + hi + hip
    # Derivatives
    yyim = -((2 * him + hi) * H + him * (him + hi)) / (him * (him + hi) * H) * yim
    yyim += ((him + hi) * H) / (him * hi * (hi + hip)) * yi
    yyim -= (him * H) / ((him + hi) * hi * hip) * yip
    yyim += (him * (him + hi)) / ((hi + hip) * hip * H) * yipp
    yyi = -(hi * (hi + hip)) / (him * (him + hi) * H) * yim
    yyi += (hi * (hi + hip) - him * (2 * hi + hip)) / (him * hi * (hi + hip)) * yi
    yyi += (him * (hi + hip)) / ((him + hi) * hi * hip) * yip
    yyi -= (him * hi) / ((hi + hip) * hip * H) * yipp
    yyip = (hi * hip) / (him * (him + hi) * H) * yim
    yyip -= (hip * (him + hi)) / (him * hi * (hi + hip)) * yi
    yyip += ((him + 2 * hi) * hip - (him + hi) * hi) / ((him + hi) * hi * hip) * yip
    yyip += ((him + hi) * hi) / ((hi + hip) * hip * H) * yipp
    yyipp = -((hi + hip) * hip) / (him * (him + hi) * H) * yim
    yyipp += (hip * H) / (him * hi * (hi + hip)) * yi
    yyipp -= ((hi + hip) * H) / ((him + hi) * hi * hip) * yip
    yyipp += ((2 * hip + hi) * H + hip * (hi + hip)) / ((hi + hip) * hip * H) * yipp
    # Smoothness indicators
    β2 = (hi + hip)^2 * (abs(yyip - yyi) / hi - abs(yyi - yyim) / him)^2
    β3 = (him + hi)^2 * (abs(yyipp - yyip) / hip - abs(yyip - yyi) / hi)^2
    α2 = - one(typeof(xim)) / ((xipp - xim) * (ε + β2))
    α3 = one(typeof(xim)) / ((xipp - xim) * (ε + β3))
    return (α2, α3)
end

"""
Compute the q2 and q3 Lagrange interpolation factors for 4th order WENO interpolation.
Arguments are the x and y values for the 4-point stencil.
"""
function weno4_q(xim, xi, xip, xipp, yim, yi, yip, yipp)
    him = xi - xim
    hi = xip - xi
    hip = xipp - xip
    q2 = Vector{eltype(xim)}(undef, 3)
    q3 = Vector{eltype(xim)}(undef, 3)
    q2[1] = yim / (him * (him + hi))
    q2[2] = -yi / (him * hi)
    q2[3] = yip / ((him + hi) * hi)
    q3[1] = yi / (hi * (hi + hip))
    q3[2] = -yip / (hi * hip)
    q3[3] = yipp / ((hi + hip) * hip)
    return (q2, q3)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Notcurses_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Notcurses")
JLLWrappers.@generate_main_file("Notcurses", UUID("56927b6a-812d-5b37-a4c6-b46ab6359194"))
end  # module Notcurses_jll


text:
using GridWorld
gameboard = zeros(Int8, 3, 3)
gameboard[3,1] = 1
gameboard[2,2] = 2
gameboard[1,3] = 1
function ttt(row,col)
    if gameboard[row,col] == 0
        return nothing
    end
    if gameboard[row,col] == 1
        return Label("X")
    end
    Label("O"; color="red")
end
g = Grid(3, 3, ttt, margin_top=20, margin_bottom=20, margin_left=20, margin_right=20)
threecol = Grid([g for _ in 1:5])
onerow = Grid([g for _ in 1:5], 1, :)
onecol = Grid([g for _ in 1:5], :, 1)
twocol = Grid([g, onecol]; valign=:top)


text:
push!(LOAD_PATH, "../src/")
using Documenter, P2P

makedocs(sitename="P2P.ai",
    pages=[
        "Introduction" => "index.md"
        "Define Problem" => "define.md"
        "How did we get here" => "path.md"
        "what do we know" => "process.md"
        "get Facts straight" => "facts.md"
        "graph" => "graph.md"
        "data insights" => "query.md"
        "analytics" => "analytics.md"
        "AI predictions" => "ml.md"
        "API" => "api.md"
    ],
)


text:
# The Patchable backend

import Patchwork
import Patchwork.Elem

# Names that can be easily imported
export Patchable

type Patchable <: Backend
    width::Float64
    height::Float64

    jsheader::Vector{String}
    jsmodules::Vector{(String, String)}
    clip_paths::Dict{ClipPrimitive, String}
    vector_properties::Vector
    function Patchable(width, height, absolute_elems=Elem[])
        width = size_measure(width)
        height = size_measure(height)
        if !isabsolute(width) || !isabsolute(height)
            error("SVG image size must be specified in absolute units.")
        end
        new(width.abs,
            height.abs,
            String[],
            (String, String)[],
            Dict{ClipPrimitive, String}(),
            Any[])
    end
end

iswithjs(::Patchable) = false
iswithousjs(::Patchable) = true

vector_properties(img::Patchable) = if !isempty(img.vector_properties)
    img.vector_properties[end]
end

clip_path_id(img, path) =
   if (haskey(img.clip_paths, path))
       return img.clip_paths[path]
   else
       id = string("clippath-", length(img.clip_paths) + 1)
       img.clip_paths[path] = id
       return id
   end

# Declare that we need the recursive drawing procedure
function draw(img::Patchable, root::Context)
    root = Elem(:svg, :svg,
         draw_recursive(img, root),
         width=string(img.width, "mm"),
         height=string(img.height, "mm"),
         viewBox=string("0 0 ", img.width, ' ', img.height),
         stroke=svg_fmt_color(default_stroke_color),
         fill=svg_fmt_color(default_fill_color))

    root &= @compat Dict("stroke-width" => svg_fmt_float(default_line_width.abs),
                 "font-size" => svg_fmt_float(default_font_size.abs))

    if !isempty(img.clip_paths)
        defs = Elem(:svg, :defs)
        for (path, id) in img.clip_paths
            defs <<= Elem(:svg, :clipPath,
                          Elem(:svg, :path, d=svg_fmt_path(path.points)),
                          id=id)
        end
        root = root << defs
    end
    root
end

root_box(img::Patchable) =
    AbsoluteBoundingBox(0.0, 0.0, img.width, img.height)

init_context(::Patchable, ::Context) = Elem(:svg, :g)

typealias SVGPart Union(Elem, Dict, Nothing)

addto(::Patchable, acc::Nothing, child::Nothing) = nothing
addto(::Patchable, acc::Nothing, child::SVGPart) = child
addto(::Patchable, acc::SVGPart, child::Nothing) = acc
addto(::Patchable, acc::Elem{:svg, :g}, child::Elem) = acc << child

addto(p::Patchable, acc::Nothing, child::Array) = addto(p, Elem(:svg, :g), child)
function addto(p::Patchable, acc::Elem{:svg, :g}, child::Array)
    for i=1:length(child)
        @inbounds acc = addto(p, acc, child[i])
    end
    acc
end

addto(::Patchable, acc::Elem, child::Elem) = Elem(:svg, :g, acc, child)
addto(::Patchable, acc::Elem, child::Dict) = acc & child
addto(::Patchable, acc::Nothing, child::Dict) = Elem(:svg, :g) & child

function push_property_frame(img::Patchable, vector_props)
    push!(img.vector_properties, vector_props)
end

function properties_at_index(img, prop_vecs, i)
    props = Dict()
    for (proptype, property) in prop_vecs
        if i > length(property.primitives)
            error("Vector of properties and vector of forms have different length")
        end
        draw!(img, property.primitives[i], props)
    end
    props
end

function pop_property_frame(img::Patchable)
    pop!(img.vector_properties)
end

# Form Drawing
# ------------

function draw(img::Patchable, form::Form)
    acc = Array(Any, length(form.primitives))
    properties = vector_properties(img)

    for i in 1:length(form.primitives)
        elem = draw(img, form.primitives[i])
        if properties !== nothing && !isempty(properties)
            props = properties_at_index(img, properties, i)
            elem = addto(img, elem, props)
        end
        acc[i] = elem
    end
    acc
end

function draw(img::Patchable, prim::BitmapPrimitive)
    warn("Patchable backend does not yet support bitmap primitives")
    nothing
end

draw(img::Patchable, prim::CirclePrimitive) =
    Elem(:svg, :circle,
         cx=prim.center.x.abs,
         cy=prim.center.y.abs,
         r=prim.radius.abs)

draw(img::Patchable, prim::CurvePrimitive) =
    Elem(:svg, :path,
         fill="none",
         path=string("M"
                     , prim.anchor0.x.abs, ','
                     , prim.anchor0.y.abs, "C"
                     , prim.ctrl0.x.abs, ','
                     , prim.ctrl0.y.abs, ' '
                     , prim.ctrl1.x.abs, ','
                     , prim.ctrl1.y.abs, ' '
                     , prim.anchor1.x.abs, ','
                     , prim.anchor1.y.abs))

function draw(img::Patchable, prim::EllipsePrimitive)
    cx = prim.center.x.abs
    cy = prim.center.y.abs
    rx = sqrt((prim.x_point.x.abs - cx)^2 +
              (prim.x_point.y.abs - cy)^2)
    ry = sqrt((prim.y_point.x.abs - cx)^2 +
              (prim.y_point.y.abs - cy)^2)
    theta = rad2deg(atan2(prim.x_point.y.abs - cy,
                          prim.x_point.x.abs - cx))

    if !all(isfinite([cx, cy, rx, ry, theta]))
        return nothing
    end
    el = Elem(:svg, :ellipse,
              cx=cx, cy=cy, rx=rx, ry=ry)
    if abs(theta) > 1e-4
        el = el & @compat Dict(:transform => string("rotate(",
                                        svg_fmt_float(theta), ' ',
                                        svg_fmt_float(cx), ' ',
                                        svg_fmt_float(cy), ')'))
    end
    el
end

function svg_fmt_path(points::Vector{Point}, bridge_gaps::Bool=false)
    io=IOBuffer()
    print_svg_path(io, points, bridge_gaps)
    takebuf_string(io)
end

function draw(img::Patchable, prim::LinePrimitive)
     n = length(prim.points)
     if n <= 1; return; end

     Elem(:svg, :path, fill="none", d=svg_fmt_path(prim.points, true))
end

function svg_fmt_path_ops(ops)
    io = IOBuffer()
    for op in ops
        svg_print_path_op(io, op)
    end
    takebuf_string(io)
end

draw(img::Patchable, prim::PathPrimitive) =
    Elem(:svg, :path, d=svg_fmt_path_ops(prim.ops))

function draw(img::Patchable, prim::PolygonPrimitive)
     n = length(prim.points)
     if n <= 1; return; end

     Elem(:svg, :path, d=svg_fmt_path(prim.points, true) * " z")
end

function draw(img::Patchable, prim::RectanglePrimitive)
    width = max(prim.width.abs, 0.01)
    height = max(prim.height.abs, 0.01)

    Elem(:svg, :rect,
         x=prim.corner.x.abs,
         y=prim.corner.y.abs,
         width=width,
         height=height)
end

function draw(img::Patchable, prim::TextPrimitive)
    el = pango_to_elems(prim.value) & @compat Dict(
            :x=>prim.position.x.abs,
            :y=>prim.position.y.abs)
    if is(prim.halign, hcenter)
        el &= @compat Dict("text-anchor" => "middle")
    elseif is(prim.halign, hright)
        el &= @compat Dict("text-anchor" => "end")
    end
    if is(prim.valign, vcenter)
        el &= @compat Dict(:dy=>"0.35em")
    elseif is(prim.halign, vtop)
        el &= @compat Dict(:dy=>"0.6mm")
    end

    if abs(prim.rot.theta) > 1e-4
        el &= @compat Dict(:transform => string("rotate(",
                                     rad2deg(prim.rot.theta), ' ',
                                     svg_fmt_float(prim.rot.offset.x.abs), ' ',
                                     svg_fmt_float(prim.rot.offset.y.abs), ')'))
     end
     el

end

# Property Primitives
# -------------------
function draw(img::Patchable, prop::Property)
    dict = Dict()
    for prim in prop.primitives
        draw!(img, prim, dict)
    end
    dict
end

function draw!(img::Patchable, prim::PropertyPrimitive, dict)
    item = draw(img, prim)
    if !is(item, nothing)
        k, v = item

        dict[k] = v
    end
end

function draw(img::Patchable, prim::ClipPrimitive)
    id = clip_path_id(img, prim)
    "clip-path", "url(#$id)"
end

draw(img::Patchable, prim::FillOpacityPrimitive) =
    :opacity, prim.value

function draw!(img::Patchable, prim::FillPrimitive, props)
    if isa(prim.color, AlphaColorValue)
        props[:fill] = svg_fmt_color(prim.color.c)
        props["fill-opacity"] = prim.color.alpha
    elseif isa(prim.color, ColorValue)
        props[:fill] = svg_fmt_color(prim.color)
    else
        props[:fill] = "none"
    end
end

draw(img::Patchable, prim::FontPrimitive) =
    "font-family", escape_string(prim.family)

draw(img::Patchable, prim::FontSizePrimitive) =
    "font-size", prim.value.abs

function draw(img::Patchable, prim::JSCallPrimitive)
    nothing
end

function draw(img::Patchable, prim::JSIncludePrimitive)
    push!(img.jsheader, prim.value)
    if prim.jsmodule != nothing
        push!(img.jsmodules, prim.jsmodule)
    end
    nothing
end

draw(img::Patchable, prim::LineWidthPrimitive) =
    "stroke-width", prim.value.abs

draw(img::Patchable, prim::SVGAttributePrimitive) = nothing
    #prim.attribute, escape_string(prim.value)

draw(img::Patchable, prim::SVGClassPrimitive) =
    :class, escape_string(prim.value)

draw(img::Patchable, prim::SVGIDPrimitive) =
    :id, escape_string(prim.value)

function draw(img::Patchable, prim::StrokeDashPrimitive)
    if isempty(prim.value)
        "stroke-dasharray", :none
    else
        "stroke-dasharray", join([v.abs for v in prim.value], ',')
    end
end

draw(img::Patchable, prim::StrokeLineCapPrimitive) =
    "stroke-linecap", svg_fmt_linecap(prim.value)

draw(img::Patchable, prim::StrokeLineJoinPrimitive) =
    "stroke-linejoin", svg_fmt_linejoin(prim.value)

draw(img::Patchable, prim::StrokeOpacityPrimitive) =
    "stroke-opacity", prim.value

function draw!(img::Patchable, prim::StrokePrimitive, dict)
    if isa(prim.color, AlphaColorValue)
        dict[:stroke] = svg_fmt_color(prim.color.c)
        dict["stroke-opacity"] = property.color.value
    else
        dict[:stroke] = svg_fmt_color(prim.color)
    end
end

draw(img::Patchable, prim::VisiblePrimitive) =
    :visibility, prim.value ? "visible" : "hidden"

# Pango markup to Patchwork
function pango_to_elems(text::String)
    pat = r"<(/?)\s*([^>]*)\s*>"
    input = text
    output = Elem[Elem(:svg, :text)] # Stack
    lastpos = 1

    baseline_shift = 0.0
    open_tag = false

    for mat in eachmatch(pat, text)
        txt = input[lastpos:mat.offset-1]
        if length(txt) > 0
            output[end] <<= Patchwork.Text(txt)
        end

        closing_tag = mat.captures[1] == "/"

        if open_tag && !closing_tag
            el = pop!(output)
            output[end] <<= el
        end

        if mat.captures[2] == "sup"
            if mat.captures[1] == "/"
                el = pop!(output)
                output[end] <<= el
            else
                el = Elem(:svg, :tspan,
                         style=(@compat Dict("dominant-baseline" => :inherit)),
                         dy="-0.6em") & @compat Dict("font-size" => "83%")
                push!(output, el)
                baseline_shift = -0.6 * 0.83
            end
        elseif mat.captures[2] == "sub"
            if mat.captures[1] == "/"
                el = pop!(output)
                output[end] <<= el
            else
                el = Elem(:svg, :tspan,
                          style=(@compat Dict("dominant-baseline" => :inherit)),
                          dy="0.6em") & @compat Dict("font-size" => "83%")
                push!(output, el)
                baseline_shift = 0.6 * 0.83
            end
        elseif mat.captures[2] == "i"
            if mat.captures[1] == "/"
                el = pop!(output)
                output[end] <<= el
            else
                el = Elem(:svg, :tspan,
                          style = @compat Dict("dominant-baseline" => :inherit)) &
                                    @compat Dict("font-style"=>"italic")
                push!(output, el)
            end
        elseif mat.captures[2] == "b"
            if mat.captures[1] == "/"
                el = pop!(output)
                output[end] <<= el
            else
                el = Elem(:svg, :tspan,
                          style=@compat Dict("dominant-baseline" => :inherit)) &
                                    @compat Dict("font-style"=>"bold")
                push!(output, el)
            end
        end

        if closing_tag && baseline_shift != 0.0
            el = Elem(:svg, :tspan,
                      dy="$(-baseline_shift)em")
            push!(output, el)
            baseline_shift = 0.0
            open_tag = true
        end

        lastpos = mat.offset + length(mat.match)
    end
    txt = input[lastpos:end]
    if length(txt) > 0
        output[end] <<= Patchwork.Text(string(txt))
    end
    if open_tag
        el = pop!(output)
        output[end] <<= el
    end
    output[1]
end

# writemime for signals
if isinstalled("Reactive")

    import Base: writemime
    import Reactive: Signal, lift

    if isdefined(Main, :IJulia)
        import IJulia: metadata
        metadata{T <: ComposeNode}(::Signal{T}) = Dict()
    end

    function writemime{T <: ComposeNode}(io::IO, m::MIME"text/html", ctx::Signal{T})
        writemime(io, m, lift(c -> draw(
            Patchable(
                Compose.default_graphic_width,
                Compose.default_graphic_height
            ), c), ctx))
    end
end


text:
function toplanetext(data::String)
   return data
end

function tojson(data::String)
   return JSON.parse(data)
end

function toxml(data::String)
   return parse_xml(content)
end

formats = Dict(
    "application/json" =>  tojson
  #, "application/xml"  =>  toxml
  , "*/*"              =>  toplanetext
  , "text/plain"       =>  toplanetext
)


text:
module Gillespie

using Distributions
using DataFrames
using QuadGK
using Roots
using StaticArrays

export
    ssa,
    gillespie,
#     jensen,
    truejump,
    ssa_data,
    pfsample,
    SSAArgs,
    SSAStats,
    SSAResult

include("SSA.jl")

"
This function performs stochastic simulation, currently using either the Doob-Gillespie algorithm or Jensen's method. It takes the following arguments:

- **x0** : a `Vector` of `Int64`, representing the initial states of the system.
- **F** : a `Function` or a callable type, which itself takes two arguments; x, a `Vector` of `Int64` representing the states, and parms, a `Vector` of `Float64` representing the parameters of the system. In the case of time-varying rates (for algorithms `:jensen` and `:tjm`), there should be a third argument, the time of the system.
- **nu** : a `Matrix` of `Int64`, representing the transitions of the system, organised by row.
- **parms** : a `Vector` of `Float64` representing the parameters of the system.
- **maxnsteps** : the maximum number of jumps (`Int64`)

There are several named arguments:

- **algo**: the algorithm to use (`Symbol`, either `:gillespie` (default), ':jensen', or ':tjm').
- **max_rate**: the maximum rate (`Float64`, for Jensen's method only).
- **thin**: (`Bool`) whether to thin jumps for Jensens method (default: `true`).

"
function ssa(x0::AbstractVector{Int64},F::Base.Callable,nu::AbstractMatrix{Int64},parms::AbstractVector{Float64},maxnsteps::Int64; algo=:gillespie, max_rate::Float64=0.0, thin::Bool=true)
  @assert algo in [:gillespie,:jensen,:tjm] "Available algorithms are :gillespie, :jensen, and :tjm"
  if algo == :gillespie
    return gillespie(x0,F,nu,parms,maxnsteps)
  end
#   if algo == :jensen
#     return jensen(x0,F,nu,parms,tf,max_rate,thin)
#   end
  if algo == :tjm
    return tjm(x0,F,nu,parms,maxnsteps)
  end
end

end # module


text:
abstract type AbstractNFFTPlan{T,D,DIM} end

##################
# define interface
##################
# in-place NFFT
@mustimplement nfft!(p::AbstractNFFTPlan{D,DIM,T}, f::AbstractArray, fHat::AbstractArray) where {D,DIM,T}
# in-place adjoint NFFT
@mustimplement nfft_adjoint!(p::AbstractNFFTPlan{D,DIM,T}, fHat::AbstractArray, f::AbstractArray) where {D,DIM,T}
# size of the input array for the NFFT
@mustimplement size(p::AbstractNFFTPlan{D,DIM,T}) where {D,DIM,T}
# size of the output array for the NFFT
@mustimplement numFourierSamples(p::AbstractNFFTPlan{D,DIM,T}) where {D,DIM,T}

#################
# derived methods
#################
"""
        nfft(p, f) -> fHat

For a **non**-directional `D` dimensional plan `p` this calculates the NFFT of a `D` dimensional array `f` of size `N`.
`fHat` is a vector of length `M`.
(`M` and `N` are defined in the plan creation)

For a **directional** `D` dimensional plan `p` both `f` and `fHat` are `D`
dimensional arrays, and the dimension specified in the plan creation is
affected.
"""
function nfft(p::AbstractNFFTPlan{D,0,T}, f::AbstractArray{U,D}, args...) where {D,T,U}
    fHat = similar(f,Complex{T}, p.M)
    nfft!(p, f, fHat, args...)
    return fHat
end

function nfft(p::AbstractNFFTPlan{D,DIM,T}, f::AbstractArray{U,D}, args...) where {D,DIM,T,U}
    sz = [p.N...]
    sz[DIM] = p.M
    fHat = similar(f, Complex{T}, Tuple(sz))
    nfft!(p, f, fHat, args...)
    return fHat
end

"""
        nfft_adjoint(p, fHat) -> f

For a **non**-directional `D` dimensional plan `p` this calculates the adjoint NFFT of a length `M` vector `fHat`
`f` is a `D` dimensional array of size `N`.
(`M` and `N` are defined in the plan creation)

For a **directional** `D` dimensional plan `p` both `f` and `fHat` are `D`
dimensional arrays, and the dimension specified in the plan creation is
affected.
"""
function nfft_adjoint(p::AbstractNFFTPlan{D,DIM,T}, fHat::AbstractArray{U}, args...) where {D,DIM,T,U}
    f = similar(fHat, Complex{T}, p.N)
    nfft_adjoint!(p, fHat, f, args...)
    return f
end

@generated function consistencyCheck(p::AbstractNFFTPlan{D,DIM,T}, f::AbstractArray{U,D},
                                     fHat::AbstractArray{Y}) where {D,DIM,T,U,Y}
    quote
        M = numFourierSamples(p)
        N = size(p)
        if $DIM == 0
            fHat_test = (M == length(fHat))
        elseif $DIM > 0
            fHat_test = @nall $D d -> ( d == $DIM ? size(fHat,d) == M : size(fHat,d) == N[d] )
        end

        if N != size(f) || !fHat_test
            throw(DimensionMismatch("Data is not consistent with NFFTPlan"))
        end
    end
end


text:
module Interface

using Test, Random, MLJBase,
      Tables, CategoricalArrays,
      OrderedCollections, DataFrames,
      TypedTables, MLJModelInterface

using ..Models

include("data_utils.jl")

include("model_api.jl")

@test trait(rms) == :measure

end

true


text:
export delimited
"""
    delimited(infile::String; delim::Union{Char,String,Regex} = "auto", digits::Int64 = 3, silent::Bool = false)
Load a delimited-type file into memory as a PopData object.

### Arguments
- `infile` : path to file

### Keyword Arguments
- `delim` : delimiter characters. By default uses auto-parsing of `CSV.File`
- `digits` : number of digits denoting each allele (default: `3`)
- `silent`   : whether to print file information during import (default: `true`)
- `allow_monomorphic::Bool` : whether to keep monomorphic loci in the dataset (default: `false`)

### File formatting:
- The first row is column names (names don't matter)
- The columns must be in this order:
    1. sample name
    2. population name
    3. longitude
    4. latitude
    5. locus_1 genotypes
    6. locus_2 genotypes
    7. etc...

#### Missing data
##### Genotypes
Missing genotypes can be formatted as all-zeros `000000`, left empty, or negative-nine `-9`

##### Location data
If location data is missing for a sample (which is ok!), make sure the value is
blank, otherwise there will be transcription errors! (look at line 3 in the example below)

## Example
```
lizardsCA = delimited("CA_lizards.csv", digits = 3);
```

### Formatting example:
```
name,population,long,lat,Locus1,Locus2,Locus3
sierra_01,mountain,11.11,-22.22,001001,-9,001001
sierra_02,mountain,11.12,-22.21,001001,001001,001002
snbarb_01,coast,,,001001,000000,001002
snbarb_02,coast,11.14,-22.24,001001,001001,001001
snbarb_03,coast,11.15,,001002,001001,001001
```
"""
function delimited(
    infile::String;
    delim::Union{Char,String,Regex} = "auto",
    digits::Int = 3,
    silent::Bool = false,
    allow_monomorphic::Bool = false
    )
    isfile(infile) || throw(ArgumentError("$infile not found."))
    dlm = delim == "auto" ? nothing : delim
    file_parse = CSV.read(infile, DataFrame, delim = dlm, missingstring = ["-9", ""])
    locinames = names(file_parse)[5:end]
    #meta = select(file_parse, 1:4)
    coords = select(file_parse, 3,4)
    # force strings for this field
    #meta.population = string.(meta.population)
    #rename!(meta, [:name, :population, :longitude, :latitude])
    select!(file_parse, Not(3:4))
    geno_parse = DataFrames.stack(file_parse, DataFrames.Not(1:2))
    rename!(geno_parse, [:name, :population, :locus, :genotype])

    select!(geno_parse,
        :name => (i -> PooledArray(i, compress = true)) => :name,
        :population => (i -> PooledArray(string.(i), compress = true)) => :population,
        :locus => (i -> PooledArray(i, compress = true)) => :locus,
        :genotype
    )

    try
        geno_parse.genotype = PooledArray(phase.(geno_parse.genotype, Int8, digits), compress = true)
    catch
        geno_parse.genotype = phase.(geno_parse.genotype, Int16, digits)
    end

    if !silent
        @info "\n $(truncatepath(abspath(infile)))\n data: loci = $(length(locinames)), samples = $(length(meta[!, 1])), populations = $(length(unique(meta[!,2])))"
        println()
    end
    popinfo = PopDataInfo(geno_parse)
    if any(.!ismissing.(coords[:,1]))
        popinfo.sampleinfo.longitude = coords[:,1]
    end
    if any(.!ismissing.(coords[:,2]))
        popinfo.sampleinfo.latitude = coords[:,2]
    end
    pd_out = PopData(popinfo, geno_parse)
    !allow_monomorphic && dropmonomorphic!(pd_out, silent = silent)
    return pd_out
end

"""
    delimited(data::PopData; filename::String, delim::String = ",", digits::Integer = 3, format::String = "wide", miss::Int = 0)
Write PopData to a text-delimited file.
### Keyword Arguments
- `filename`: a `String` of the output filename
- `digits` : an `Integer` indicating how many digits to format each allele as (e.g. `(1, 2)` => `001002` for `digits = 3`)
- `format` : a `String` indicating whether to output in`"wide"` or `"long"` (aka `"tidy"`) format
  - `wide` : the standard format CSV for importing into PopGen.jl
  - `long` : the `loci` table with `longitude` and `latitude` columns added
- `delim` : the `String` delimiter to use for writing the file
- `miss` : an `Integer` for how you would like missing values written
    - `0` : As a genotype represented as a number of zeroes equal to `digits × ploidy` like `000000` (default)
    - `-9` : As a single value `-9`

### Example
```julia
cats = @nancycats;
fewer_cats = omit(cats, name = samplenames(cats)[1:10]);
delimited(fewer_cats, filename = "filtered_nancycats.gen", digits = 3, format = "wide", delim = " ")
```
"""
function delimited(data::PopData; filename::String, delim::String = ",", digits::Integer = 3, format::String = "wide", miss::Int = 0)
    # create empty String column
    unphased_df = insertcols!(copy(data.genodata), :string_geno => fill("", length(data.genodata.name)))
    # grouping to facilitate pulling out ploidy values for coding missing genotypes
    grp_df = groupby(unphased_df, :name)
    for sample in grp_df
        samplename = sample.name[1]
        sample_ploidy = convert(Int, data.sampleinfo.ploidy[data.sampleinfo.name .== samplename][1])
        sample.string_geno .= unphase.(sample.genotype, digits = digits, ploidy = sample_ploidy, miss = miss)
    end
    # pop out original genotype column
    select!(unphased_df, 1, 2, 3, 5)
    if format == "wide"
        wide_df = DataFrames.unstack(unphased_df, :locus, :string_geno)
        insertcols!(wide_df, 3, :longitude => data.metadata.longitude, :latitude => data.metadata.latitude)
        CSV.write(filename, wide_df, delim = delim) ;
        return
    else
        out_df = sort(
            transform(
                unphased_df, 1:2,
                3 => (i -> Vector{Union{Float32, Missing}}(undef, length(i))) => :longitude,
                4 => (i -> Vector{Union{Float32, Missing}}(undef, length(i))) => :latitude, 3:4
                )
                    )
        for i in 1:length(data.sampleinfo.name)
            out_df[out_df.name .== data.sampleinfo.name[i], :longitude] .= data.metadata.longitude[i]
            out_df[out_df.name .== data.sampleinfo.name[i], :latitude] .= data.metadata.latitude[i]
        end
        CSV.write(filename, out_df, delim = delim) ;
        return
    end
end


text:
using DNABarcode
using Base.Test
using Bio.Seq

function hamming_distance(x, y)
    @assert length(x) == length(y)
    dist = 0
    for i in 1:endof(x)
        dist += x[i] != y[i]
    end
    return dist
end

function randdna()
    r = rand()
    if r < 0.20
        return DNA_A
    elseif r < 0.40
        return DNA_C
    elseif r < 0.60
        return DNA_G
    elseif r < 0.80
        return DNA_T
    else
        return DNA_N
    end
end

@testset "HammingCode" begin
    for k in 2:6, mindist in 2:k
        barcodes = genbarcodes(HammingCode{k}, mindist)
        @test !isempty(barcodes)
        ok = true
        for i in 1:length(barcodes)
            for j in i+1:length(barcodes)
                if hamming_distance(barcodes[i], barcodes[j]) < mindist
                    ok = false
                end
            end
        end
        @test ok

        for i in rand(1:length(barcodes), 10)
            barcode = barcodes[i]
            seq = DNASequence(barcode) * dna"ACGTAT"
            @test demultiplex(barcodes, seq) == i
            for _ in 1:fld(mindist - 1, 2)
                seq[rand(1:k)] = randdna()
            end
            @test demultiplex(barcodes, seq) == i
        end
    end
end


text:
using Test
using CausalDiscovery
using Random
using MLStyle

aexpr = au"""(program
(= GRID_SIZE 16)

(object Particle (Cell 0 0 "blue"))

(: particles (List Particle))
(= particles
(initnext (list)
(updateObj (prev particles) (--> obj (Particle (uniformChoice (adjPositions (.. obj origin))))))))

(on clicked (= particles (addObj (prev particles) (Particle (Position (.. click x) (.. click y))))))
)"""

aexpr2 = au"""(program
(= GRID_SIZE 16)

(object Particle (Cell 0 0 "blue"))

(: particles (List Particle))
(= particles
(initnext (list)
(updateObj (prev particles) (--> obj (Particle (uniformChoice (adjPositions (.. obj origin))))))))

(on clicked (= particles (addObj (prev particles) (Particle (Position (.. click x) (.. click y))))))
)"""

aexpr3 = au"""(program
(= GRID_SIZE 16)

(object Particle (Cell 0 0 "red"))

(: particles (List Particle))
(= particles
(initnext (list)
(updateObj (prev particles) (--> obj (Particle (uniformChoice (adjPositions (.. obj origin))))))))

(on clicked (= particles (addObj (prev particles) (Particle (Position (.. click x) (.. click y))))))
)"""

aumod = eval(compiletojulia(aexpr))
state = aumod.init(nothing, nothing, nothing, nothing, nothing, MersenneTwister(0))
state = aumod.next(state, nothing, nothing, nothing, nothing, nothing)
state = aumod.next(state, aumod.Click(5,5), nothing, nothing, nothing, nothing)

# p = eval(compiletojulia(aexpr))
aumod2 = eval(compiletojulia(aexpr2))

e = [(aum -> [aum.Cell(aum.Position(5, 5), "blue", 0.8)],
[[nothing, nothing, nothing, nothing, nothing, MersenneTwister(0)],
[nothing, nothing, nothing, nothing, nothing],
[aumod2.Click(5,5), nothing, nothing, nothing, nothing]]),
(aum -> [aum.Cell(aum.Position(6, 6), "blue", 0.8), aum.Cell(aum.Position(8, 8), "red", 0.8), aum.Cell(aum.Position(5, 5), "blue", 0.8)],
[[nothing, nothing, nothing, nothing, nothing, MersenneTwister(0)],
[nothing, nothing, nothing, nothing, nothing],
[aumod2.Click(5,4), nothing, nothing, nothing, nothing]])]

@inferred get_fitness(aumod2, e, aexpr2)
fitness_p1 = get_fitness(aumod2, e, aexpr2)

e = [(aum -> [aum.Cell(aum.Position(5, 5), "blue", 0.8)],
[[nothing, nothing, nothing, nothing, nothing, MersenneTwister(0)],
[nothing, nothing, nothing, nothing, nothing],
[aumod2.Click(5,5), nothing, nothing, nothing, nothing]])]

@inferred get_fitness(aumod2, e, aexpr2)
fitness_p2 = get_fitness(aumod2, e, aexpr2)

@test fitness_p1 > fitness_p2

aumod3 = eval(compiletojulia(aexpr3))

e = [(aum -> [aum.Cell(aum.Position(5, 5), "blue", 0.8)],
[[nothing, nothing, nothing, nothing, nothing, MersenneTwister(0)],
[nothing, nothing, nothing, nothing, nothing],
[aumod3.Click(5,5), nothing, nothing, nothing, nothing]])]

@inferred get_fitness(aumod3, e, aexpr3)
fitness_p3 = get_fitness(aumod3, e, aexpr3)

@test fitness_p3 > fitness_p1
@test fitness_p3 > fitness_p2


text:
#==========================================================================
  Regularized Linear Discriminant Analysis Solver
==========================================================================#

immutable ModelLDA{T<:BlasReal}
    order::Union{Type{Val{:row}},Type{Val{:col}}}
    is_cda::Bool
    W::AbstractMatrix{T}  # Whitening matrix
    M::Matrix{T}          # Matrix of class means (one per row)
    priors::Vector{T}     # Vector of class priors
    gamma::Nullable{T}
end

function show(io::IO, model::ModelLDA)
    println(io, (model.is_cda ? "Canonical" : "Linear") * " Discriminant Model")
    println(io, "\nRegularization Parameters:")
    print(io, "    γ = ")
    isnull(model.gamma) ? print(io, "N/A") : showcompact(get(model.gamma))
    println(io, "\n\nClass Priors:")
    for i in eachindex(model.priors)
        print(io, "    Class ", i, " Probability: ")
        showcompact(io, model.priors[i]*100)
        println(io, "%")
    end
    println(io, "\nClass Means:")
    for j in eachindex(model.priors)
        print(io, "    Class ", j, " Mean: [")
        print(io, join([sprint(showcompact, v) for v in subvector(model.order, model.M, j)], ", "))
        println(io, "]")
    end
    print(io, "\n")
end

# Extracts Fisher components
function components!{T<:BlasReal}(
         ::Type{Val{:row}},
        W::AbstractMatrix{T},
        M::Matrix{T},
        w::Vector{T}
    )
    H = broadcast!(-, M, M, w'M)
    UDVᵀ = svdfact!(H * W)
    d = min(size(M,2),size(M,1)-1)
    Vᵀ = sub(UDVᵀ[:Vt], 1:d, :) # (UDVᵀ[:Vt])[1:d,:]
    W * transpose(Vᵀ)
end

function components!{T<:BlasReal}(
         ::Type{Val{:col}},
        W::AbstractMatrix{T},
        M::Matrix{T},
        w::Vector{T}
    )
    H = broadcast!(-, M, M, M*w)
    UDVᵀ = svdfact!(transpose(W * H))
    d = min(size(M,1), size(M,2)-1)
    Vᵀ = sub(UDVᵀ[:Vt], 1:d, :) #(UDVᵀ[:Vt])[1:d,:]
    Vᵀ * W
end

for (scheme, dim_obs) in ((:(:row), 1), (:(:col), 2))
    isrowmajor = dim_obs == 1
    dim_param = isrowmajor ? 2 : 1

    _ij, _ji = isrowmajor ? (:i, :j) : (:j, :i)  # Swapped variables for row and column ordering
    _nk, _kn = isrowmajor ? (:n, :k) : (:k, :n)
    _np, _pn = isrowmajor ? (:n, :p) : (:p, :n)
    W, H     = isrowmajor ? (:W, :H) : (:H, :W)

    @eval begin
        # X in uncentered data matrix
        # M is matrix of class means (one per row)
        # y is one-based vector of class IDs
        # γ is nullable regularization parameter in [0,1]
        function lda!{T<:BlasReal,U}(
                 ::Type{Val{$scheme}},
                X::Matrix{T},
                M::Matrix{T},
                y::RefVector{U},
                γ::Nullable{T}
            )
            p = size(X, $dim_param)
            p >= 1 || throw(DimensionMismatch("Must have at least one parameter"))
            n = size(X, $dim_obs)
            n >= 2 || throw(DimensionMismatch("Must have at least two observations"))
            H = centerclasses!(Val{$scheme}, X, M, y)
            isnull(γ) ? whitendata_qr!(Val{$scheme}, H) : whitendata_svd!(Val{$scheme}, H, get(γ))
        end

        # same rules as lda! for common arguments
        # priors is a vector of class weights
        function cda!{T<:BlasReal,U}(
                 ::Type{Val{$scheme}},
                X::Matrix{T},
                M::Matrix{T},
                y::RefVector{U},
                γ::Nullable{T},
                priors::Vector{T}
            )
            k = convert(Int64, y.k)
            if k != length(priors)
                throw(DimensionMismatch("Argument priors length does not match class count"))
            end
            W = lda!(Val{$scheme}, X, M, y, γ)
            size(X, $dim_param) < k ? W : components!(Val{$scheme}, W, M, priors)
        end

        function discriminants_lda{T<:BlasReal}(
                 ::Type{Val{$scheme}},
                W::AbstractMatrix{T},
                M::Matrix{T},
                priors::Vector{T},
                Z::Matrix{T}
            )
            n = size(Z, $dim_obs)
            p = size(Z, $dim_param)
            d = size(W, $dim_param)
            k = length(priors)
            D   = Array(T, $_nk, $_kn)  # discriminant function values
            H   = Array(T, $_np, $_pn)  # temporary array to prevent re-allocation k times
            hᵀh = Array(T, n)           # diag(H'H)
            for j = 1:k
                broadcast!(-, H, Z, subvector(Val{$scheme}, M, j))
                dotvectors!(Val{$scheme}, $H * $W, hᵀh)
                for i = 1:n
                    D[$_ij, $_ji] = -hᵀh[i]/2 + log(priors[j])
                end
            end
            D
        end
    end
end

doc"`lda(X, y; order, M, gamma, priors)` Fits a regularized linear discriminant model to the data
in `X` based on class identifier `y`."
function lda{T<:BlasReal,U<:Integer}(
        X::Matrix{T},
        y::AbstractVector{U};
        order::Union{Type{Val{:row}},Type{Val{:col}}} = Val{:row},
        M::Matrix{T} = classmeans(order, X, RefVector(y)),
        gamma::Union{T,Nullable{T}} = Nullable{T}(),
        priors::Vector{T} = ones(T,maximum(y))/maximum(y)
    )
    all(priors .> 0) || error("Argument priors must have positive values only")
    isapprox(sum(priors), one(T)) || error("Argument priors must sum to 1")
    yref = isa(y,RefVector) ? y : RefVector(y)
    length(priors) == yref.k || throw(DimensionMismatch("Prior length must match class count"))
    γ = isa(gamma, Nullable) ? gamma : Nullable(gamma)
    W = lda!(order, copy(X), M, yref, γ)
    ModelLDA{T}(order, false, W, M, priors, γ)
end

doc"`cda(X, y; order, M, gamma, priors)` Fits a regularized canonical discriminant model to the
data in `X` based on class identifier `y`."
function cda{T<:BlasReal,U<:Integer}(
        X::Matrix{T},
        y::AbstractVector{U};
        order::Union{Type{Val{:row}},Type{Val{:col}}} = Val{:row},
        M::Matrix{T} = classmeans(order, X, RefVector(y)),
        gamma::Union{T,Nullable{T}} = Nullable{T}(),
        priors::Vector{T} = ones(T, maximum(y))/maximum(y)
    )
    all(priors .> 0) || error("Argument priors must have positive values only")
    isapprox(sum(priors), one(T)) || error("Argument priors must sum to 1")
    yref = isa(y,RefVector) ? y : RefVector(y)
    length(priors) == yref.k || throw(DimensionMismatch("Prior length must match class count"))
    γ = isa(gamma, Nullable) ? gamma : Nullable(gamma)
    W = cda!(order, copy(X), copy(M), yref, γ, priors)
    ModelLDA{T}(order, true, W, M, priors, γ)
end

function discriminants{T<:BlasReal}(mod::ModelLDA{T}, Z::Matrix{T})
    discriminants_lda(mod.order, mod.W, mod.M, mod.priors, Z)
end

function classify{T<:BlasReal}(mod::ModelLDA{T}, Z::Matrix{T})
    vec(mapslices(indmax, discriminants(mod, Z), mod.order == Val{:row} ? 2 : 1))
end


text:
export IMatrix

"""
    IMatrix{N, Tv}()
    IMatrix{N}() where N = IMatrix{N, Int64}()
    IMatrix(A::AbstractMatrix{T}) where T -> IMatrix

IMatrix matrix, with size N as label, use `Int64` as its default type, both `*` and `kron` are optimized.
"""
struct IMatrix{N,Tv} <: AbstractMatrix{Tv} end
IMatrix{N}() where {N} = IMatrix{N,Bool}()
IMatrix(N::Int) = IMatrix{N}()

size(A::IMatrix{N}, i::Int) where {N} = N
size(A::IMatrix{N}) where {N} = (N, N)
getindex(A::IMatrix{N,T}, i::Integer, j::Integer) where {N,T} = T(i == j)

####### sparse matrix ######
nnz(M::IMatrix{N}) where {N} = N
nonzeros(M::IMatrix{N,T}) where {N,T} = ones(T, N)
findnz(M::IMatrix{N,T}) where {N,T} = (collect(1:N), collect(1:N), ones(T, N))
ishermitian(D::IMatrix) = true
isdense(::IMatrix) = false

similar(::IMatrix{N,Tv}, ::Type{T}) where {N,Tv,T} = IMatrix{N,T}()
copyto!(A::IMatrix{N}, B::IMatrix{N}) where {N} = A


text:
"""
    EpsilonGreedySelector <: AbstractActionSelector
    EpsilonGreedySelector(ϵ)

The best action is selected for a proportion ``1 - \\epsilon``
and a random action (with uniform probability) is selected for a proportion ``\\epsilon``.
ϵ can also be a decay. See the following examples.

## Example

```julia
julia> selector = EpsilonGreedySelector(0.1)
EpsilonGreedySelector(0.1)

julia> countmap(selector([1,2,1,1]) for _ in 1:1000)
Dict{Any,Int64} with 4 entries:
  4 => 37
  2 => 915
  3 => 22
  1 => 26
```

julia> ϵ = exp_decay(init=1.0, λ=1.0, decay_step=500, clip=0.1)
(::getfield(Ju, Symbol("#f#34")){Float64,Float64,Int64,Float64}) (generic function with 1 method)

julia> selector = EpsilonGreedySelector(ϵ)
EpsilonGreedySelector{getfield(Ju, Symbol("#f#34")){Float64,Float64,Int64,Float64}}(getfield(Ju, Symbol("#f#34")){Float64,Float64,Int64,Float64}(1.0, 1.0, 500, 0.1, Core.Box(-1)))

julia> countmap(selector([1,2,1,1]) for _ in 1:1000)
Dict{Any,Int64} with 4 entries:
  4 => 101
  2 => 677
  3 => 106
  1 => 116
"""
mutable struct EpsilonGreedySelector{T<:Union{Float64, <:Function}} <: AbstractActionSelector
    ϵ::T
end

_epsilon_select(values, ϵ) = rand() > ϵ ? sample(findallmax(values)[2]) : rand(1:length(values))

"""
    (p::EpsilonGreedySelector)(values::AbstractArray{T, 1}) where T

!!! note
    If multiple values with the same maximum value are found.
    Then a random one will be returned!

    `NaN` will be filtered unless all the values are `NaN`.
    In that case, a random one will be returned.
"""
(p::EpsilonGreedySelector{Float64})(values) = _epsilon_select(values, p.ϵ)

(p::EpsilonGreedySelector{<:Function})(values) = _epsilon_select(values, p.ϵ())

# TODO: add `dims` argument to support higher dimension?


text:
using BlackBoxOptim: pareto_dominates, pareto_dominates_hat

describe("Pareto dominance") do
  test("simple 1-objective cases where there is dominance") do
    @check pareto_dominates([0.0], [1.0])
    @check pareto_dominates([-100.0], [100.0])
  end

  test("simple 1-objective cases where there is NO dominance") do
    @check !pareto_dominates([-2.5], [-2.5])
    @check !pareto_dominates([12.32], [12.32])
  end

  test("simple 2-objective cases where there is dominance") do
    @check pareto_dominates([0, 1], [1, 2])
    @check pareto_dominates([-10, -5], [0, 0])
    @check pareto_dominates([0, 1], [0, 2])
    @check pareto_dominates([-1, 2], [0, 2])
  end

  test("simple 2-objective cases where there is NO dominance") do
    @check !pareto_dominates([0, 1], [0, 1])
    @check !pareto_dominates([-1, 1], [-1, 1])

    @check !pareto_dominates([0, 1], [1, 0])
  end

  rand_vector(size = rand(1:107)) = rand(1:103) * randn(size)
  rand_indices(n, numindices) = shuffle(collect(1:n))[1:numindices]

  # Generate a delta-vector for a vector, i.e. a vector of changes in a number of nonzero positions.
  delta_vector_for(u, numNonzero) = begin
    size = length(u)
    delta = abs(rand_vector(size))
    delta[rand_indices(size, size - min(numNonzero, size))] = 0.0
    delta
  end

  @repeat test("when the vectors are equal, i.e. NO dominance") do
    u = rand_vector()
    @check !pareto_dominates(u, u)
    @check pareto_dominates_hat(u, u) == 0
  end

  @repeat test("when the vectors have one larger and one smaller value, i.e. NO dominance") do
    u = rand_vector(rand(2:107))
    v = copy(u)
    idxs = rand_indices(length(u), 2)
    v[idxs[1]] += rand()
    v[idxs[2]] -= rand()
    @check !pareto_dominates(u, v)
    @check pareto_dominates_hat(u, v) == 0
  end

  @repeat test("when there is one dominating") do
    u = rand_vector()
    v = copy(u) .+ abs(delta_vector_for(u, rand(1:length(u))))
    @check pareto_dominates(u, v)
    @check !pareto_dominates(v, u)
    @check pareto_dominates_hat(u, v) == -1
    @check pareto_dominates_hat(v, u) ==  1
  end

  @repeat test("for single float fitness values") do
    fval1 = randn()
    f1 = makefitness(fval1)
    fval2 = fval1 + rand(0.0001:113.0)
    f2 = makefitness(fval2)
    @check pareto_dominates(f1, f2)
    @check !pareto_dominates(f2, f1)
    @check !pareto_dominates(f1, f1)
    @check !pareto_dominates(f2, f2)

    fval3 = fval1 - rand(0.00001:113.0)
    f3 = makefitness(fval3)
    @check pareto_dominates(f3, f1)
    @check pareto_dominates(f3, f2)
    @check !pareto_dominates(f1, f3)
    @check !pareto_dominates(f2, f3)
    @check !pareto_dominates(f3, f3)
  end
end


text:
function solution(str::String)::Array{String}
    # unicode safe
    result = Array{String}(undef,(1+length(str))÷2)
    i = 0
    for c in str
        i += 1
        if isodd(i)
            result[i÷2+1] = String([c])
        else
            result[i÷2] *= String([c])
        end
    end
    if isodd(i)
        result[end] *= "_"
    end
    result
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas