text:
module Kruskal

using LinearAlgebra
using Statistics
using Optim

export kruskalmds

# naming similar to [Kruskal, Nonmetric Multidimensional Scaling: A Numerical Method]
# X is a vector with values [x₁,y₁,...,x₂,y₂,...,xₙ,yₙ,...], i.e. a vectorized pxN matrix

# gradient-free S, for debugging purposes
# function stressreference{T}(X::AbstractVecOrMat{T},p,D)
#     N = div(length(X),p)

#     num = den = zero(T)

#     # for i>j
#     ind = 1 # index in vectorization of upper triangular part
#     for j=2:N
#         for i=1:j-1
#             i2 = p*(i-1) # 0-based index of xᵢ in X
#             j2 = p*(j-1) # 0-based index of xⱼ in X

#             t = zero(T)
#             for dim=1:p
#                 t += (X[i2+dim]-X[j2+dim])^2
#             end
#             num += (sqrt(t)-D[ind])^2
#             den += t

#             ind += 1
#         end
#     end
#     sqrt(num/den)
# end

mutable struct SData{T}
    lastX::Vector{T} # Vectorization of pxN data point matrix. Needed to figure out whether we moved or not since last call to _update
    p::Int # Number of dimensions in low-dimensional representation.
    N::Int # Number of data points
    targetD::Vector{T} # Vectorization of upper triangular part of target distance matrix.
    num::T
    den::T
    S::T
    numGrad::Vector{T} # gradient of numerator
    denGrad::Vector{T} # gradient of denominator
end
SData(p,N,targetD::Vector{T}) where {T} = SData(NaN*zeros(T,p*N),p,N,targetD,zero(T),zero(T),zero(T),Vector{T}(undef,p*N),Vector{T}(undef,p*N))

function update!(sData::SData{T},X) where {T}
    sData.lastX == X && return # nothing to do

    targetD = sData.targetD
    N = sData.N # number of data points
    p = sData.p
    numGrad = sData.numGrad
    denGrad = sData.denGrad
    num = den = zero(T)
    numGrad[:] .= 0
    denGrad[:] .= 0

    # for i<j, point indices
    ind = 1 # index in vectorization of upper triangular part
    for j=2:N
        for i=1:j-1
            i2 = p*(i-1) # 0-based index of xᵢ in X
            j2 = p*(j-1) # 0-based index of xⱼ in X

            t = zero(T)
            for dim=1:p
                t += (X[i2+dim]-X[j2+dim])^2
            end
            st = sqrt(t)
            num += (st-targetD[ind])^2
            den += t

            u = 2*(st-targetD[ind])/st # common computations for all dimensions in numGrad
            for dim=1:p
                v = (X[i2+dim]-X[j2+dim])
                #w = u*v
                w = abs(v)<1e-12 ? 0.0 : u*v # handle edge case by letting 0*Inf := 0
                numGrad[i2+dim] += w
                numGrad[j2+dim] -= w
                denGrad[i2+dim] += 2v
                denGrad[j2+dim] -= 2v
            end

            ind += 1
        end
    end

    sData.num = num
    sData.den = den
    sData.S = sqrt(num/den)
    nothing
end

# gradient-aware S (i.e. sharing data with gradient computation)
function stress(X,sData::SData)
    update!(sData,X)
    sData.S
end

function stressgradient(X,storage,sData::SData)
    update!(sData,X)

    targetD = sData.targetD
    num = sData.num
    den = sData.den
    S = sData.S
    numGrad = sData.numGrad
    denGrad = sData.denGrad

    C = 1.0/(2*S) # constant factor
    for i=1:length(X)
        storage[i] = C*(numGrad[i]*den - num*denGrad[i]) / den^2
    end

    nothing
end

function updatelocation(D::AbstractVector{T}, X0::AbstractMatrix{T}, tol) where {T}
    p = size(X0,1)
    N = size(X0,2)
    X0 = X0[:]

    # # Simple, gradient-free version
    # res = optimize(x->stressreference(x,p,D), X0, BFGS())
    # @assert Optim.converged(res)
    # X = reshape(Optim.minimizer(res),p,N)

    sData = SData(p,N,D) # init storage for reused computations
    # res = optimize(x->stress(x,sData), (storage,x)->stressgradient(x,storage,sData), X0, BFGS())
    res = optimize(x->stress(x,sData), (storage,x)->stressgradient(x,storage,sData), X0, BFGS(), Optim.Options(x_tol=tol, f_tol=tol))
    #res = optimize(x->stress(x,sData), X0, BFGS()) # doesn't use gradient, for comparison
    Optim.converged(res) || error("Kruskal MDS: Location optimization did not coverge.")
    X = reshape(Optim.minimizer(res),p,N)

    # Since the optimum is preserved when X is uniformly scaled, normalize scaling to keep numerical stability
    X = X .- mean(X;dims=2) # center X (doesn't change stress)
    X = X / sqrt(mean(sum(X.^2;dims=1);dims=2)[1]) # scale so that RMS distance to origin is equal to 1

    X, Optim.minimum(res)
end

# convert distance matrix to pdist-style vector
function distancevector(D::AbstractMatrix)
    # vectorize upper triangular part of D
    N = size(D,1)
    dist = zeros(div(N*(N-1),2))
    ind=1
    for j=2:N
        for i=1:j-1
            dist[ind] = D[i,j]
            ind += 1
        end
    end
    dist
end

# convert pdist-style distance vector to matrix
function distancematrix!(D::AbstractMatrix, dist::AbstractVector)
    N = size(D,1)

    for i=1:N
        D[i,i] = 0.0
    end

    ind=1
    for j=2:N
        for i=1:j-1
            D[i,j] = D[j,i] = dist[ind]
            ind += 1
        end
    end
    D
end
distancematrix(dist::AbstractVector{T}) where {T} = (M=length(dist); N=round(Int,1/2+sqrt(1/4+2M)); distancematrix!(Matrix{T}(undef,N,N),dist))

function isotonicregression(Y::AbstractVector{T}, tieGroups::Vector{UnitRange}=[]) where {T}
    X = Tuple{T,Int}[] # tuple with values and integer weight

    # reorder values in each tieGroup to put them in increasing order
    isempty(tieGroups) || (Y=copy(Y)) # avoid modifying the original vector
    tiePerms = map!(rng->sortperm(Y[rng]), Vector{Vector{Int}}(undef,length(tieGroups)), tieGroups)
    for (g,p) in zip(tieGroups,tiePerms)
        Y[g] = Y[g[p]]
    end

    for y in Y
        yw = 1

        while !isempty(X) && y<=X[end][1]
            # merge with previous point
            x,xw = pop!(X)
            y,yw = (x*xw+y*yw)/(xw+yw), xw+yw
        end
        push!(X, (y,yw)) # store point
    end

    # unwrap to full vector by duplicating elements by their weight
    Z = Vector{T}(undef,length(Y))
    zi = 1
    for x in X
        Z[zi:zi+x[2]-1] .= x[1]
        zi += x[2]
    end

    # inverse tieGroup reordering
    for (g,p) in zip(tieGroups,tiePerms)
        Z[g[p]] = Z[g]
    end

    Z
end

function pdist(X)
    # create pdist-style distance vector
    p = size(X,1)
    N = size(X,2)
    dist = similar(X, div(N*(N-1),2))
    ind=1
    for j=2:N
        for i=1:j-1
            dist[ind] = sqrt(sum((X[:,i]-X[:,j]).^2)) # D[i,j]
            ind += 1
        end
    end
    dist
end

function duplicateranges(v)
    @assert issorted(v)
    count = 0
    out = UnitRange[]
    for i=1:length(v)
        if i<length(v) && v[i]==v[i+1]
            count+=1
        else
            count>0 && push!(out, i-count:i)
            count = 0
        end
    end
    out
end

function kruskalmds(D::AbstractMatrix{T}, p::Integer, X0=convert(Matrix{T},randn(size(D,1),p)); maxIter=1000, tol=1e-6) where {T<:AbstractFloat}
    @assert p>0
    @assert D≈D' "Distance matrix must be symmetric"
    N = size(D,1)
    @assert size(X0)==(N,p)
    X0 = X0' # work with transpose internally.

    dist = distancevector(D)
    perm = sortperm(dist)
    iperm = invperm(perm)
    tieGroups = duplicateranges(dist[perm])

    X = X0 .- mean(X0;dims=2) # center X (doesn't change stress)
    X = X / sqrt(mean(sum(X.^2;dims=1);dims=2)[1]) # scale so that RMS distance to origin is equal to 1 (doesn't change stress, if the monotone function is updated accordingly)

    S = Inf
    converged = false

    # println("i: 0, S=", Sreference(X0,p,dist))

    for i=1:maxIter # nbr iterations
        prevS = S
        modelDist = pdist(X)[perm] # compute distance and permute
        permDist = isotonicregression(modelDist, tieGroups) # nonnegativity constraint will be enforced automatically since all values in modelDist are nonnegative
        dist = permDist[iperm]

        XOld = X
        X,S = updatelocation(dist,X,tol)
        stepSize = norm(X-XOld)

        # println("i: ", i, ", S=", S, ", stepSize=", stepSize, ", prevS-S=", prevS-S, ", tol*S=", tol*S)
        stepSize<tol && (converged=true; break)
        prevS-S<tol*S && (converged=true; break)
        prevS = S
    end

    # post-processing
    # the stress is invariant under rotation, so we can normalize rotation by principal axes
    F = svd(X)
    X = F.V*Diagonal(F.S)

    X, S, converged, distancematrix(dist)
end
kruskalmds(D::AbstractMatrix, p::Integer, X0=randn(size(D,1),p); kwargs...) = kruskalmds(convert(Matrix{Float64},D),p,X0;kwargs...)

end


text:
# Test general queries
@testset "General Queries" begin
    # Setup the infinite model
    optimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                         eval_objective_value=false)
    m = InfiniteModel(optimizer)
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_variable(m, inf(par))
    @hold_variable(m, g)
    # setup the mock solved transcription model
    tm = transcription_model(m)
    @variable(tm, x)
    @variable(tm, y)
    @objective(tm, Min, x^2)
    @constraint(tm, c1, 2 * x * y <= 1)
    tm.ext[:TransData].hold_to_var[g] = x
    tm.ext[:TransData].infinite_to_vars[inf] = [x, y]
    modelstring = """
    variables: x, y
    minobjective: 1*x*x
    c1: 2*x*y <= 1.0
    """
    model = MOIU.Model{Float64}()
    MOIU.loadfromstring!(model, modelstring)
    MOIU.test_models_equal(JuMP.backend(tm).model_cache, model, ["x","y"], ["c1"])
    JuMP.optimize!(tm)

    mockoptimizer = JuMP.backend(tm).optimizer.model
    MOI.set(mockoptimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    MOI.set(mockoptimizer, MOI.RawStatusString(), "solver specific string")
    MOI.set(mockoptimizer, MOI.ResultCount(), 2)
    MOI.set(mockoptimizer, MOI.PrimalStatus(1), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.DualStatus(1), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.PrimalStatus(2), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.DualStatus(2), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.SolveTime(), 0.42)
    # test termination_status
    @testset "JuMP.termination_status" begin
        @test termination_status(m) == MOI.OPTIMAL
    end
    # test primal_status
    @testset "JuMP.primal_status" begin
        @test primal_status(m) == MOI.FEASIBLE_POINT
        @test primal_status(m, result = 2) == MOI.FEASIBLE_POINT
    end
    # test dual_status
    @testset "JuMP.dual_status" begin
        @test dual_status(m) == MOI.FEASIBLE_POINT
        @test dual_status(m, result = 2) == MOI.FEASIBLE_POINT
    end
    # test solve_time
    @testset "JuMP.solve_time" begin
        @test solve_time(m) == 0.42
    end
    # test raw_status
    @testset "JuMP.raw_status" begin
        @test raw_status(m) == "solver specific string"
    end
end

# Test objective queries
@testset "Objective Queries" begin
    # Setup the infinite model
    optimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                         eval_objective_value=false)
    m = InfiniteModel(optimizer)
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_variable(m, inf(par))
    @hold_variable(m, g)
    # setup the mock solved transcription model
    tm = transcription_model(m)
    @variable(tm, x)
    @variable(tm, y)
    @objective(tm, Min, x^2)
    @constraint(tm, c1, 2 * x * y <= 1)
    tm.ext[:TransData].hold_to_var[g] = x
    tm.ext[:TransData].infinite_to_vars[inf] = [x, y]
    modelstring = """
    variables: x, y
    minobjective: 1*x*x
    c1: 2*x*y <= 1.0
    """
    model = MOIU.Model{Float64}()
    MOIU.loadfromstring!(model, modelstring)
    MOIU.test_models_equal(JuMP.backend(tm).model_cache, model, ["x","y"], ["c1"])
    JuMP.optimize!(tm)

    mockoptimizer = JuMP.backend(tm).optimizer.model
    MOI.set(mockoptimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    MOI.set(mockoptimizer, MOI.ObjectiveValue(), -1.0)
    MOI.set(mockoptimizer, MOI.DualObjectiveValue(), -2.0)
    MOI.set(mockoptimizer, MOI.ObjectiveBound(), 2.0)
    MOI.set(mockoptimizer, MOI.ResultCount(), 1)
    MOI.set(mockoptimizer, MOI.PrimalStatus(), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.DualStatus(), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.ConstraintDual(), JuMP.optimizer_index(c1), -2.0)
    @testset "JuMP.objective_bound" begin
        @test objective_bound(m) == 2.0
    end
    # test objective_bound
    @testset "JuMP.objective_value" begin
        @test objective_value(m) == -1.
    end
    # test dual_objective_bound
    @testset "JuMP.dual_objective_value" begin
        @test dual_objective_value(m) == -2.
    end
end

# Test variable queries
@testset "Variable Queries" begin
    # Setup the infinite model
    optimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                         eval_objective_value=false)
    m = InfiniteModel(optimizer)
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_variable(m, inf(par))
    @hold_variable(m, g)
    # setup the mock solved transcription model
    tm = transcription_model(m)
    @variable(tm, x)
    @variable(tm, y)
    @objective(tm, Min, x^2)
    @constraint(tm, c1, x <= 0)
    tm.ext[:TransData].hold_to_var[g] = x
    tm.ext[:TransData].infinite_to_vars[inf] = [x, y]
    modelstring = """
    variables: x, y
    minobjective: 1*x*x
    c1: x <= 0.0
    """
    model = MOIU.Model{Float64}()
    MOIU.loadfromstring!(model, modelstring)
    # MOIU.test_models_equal(JuMP.backend(tm).model_cache, model, ["x","y"], ["c1"])
    JuMP.optimize!(tm)

    mockoptimizer = JuMP.backend(tm).optimizer.model
    MOI.set(mockoptimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    MOI.set(mockoptimizer, MOI.ResultCount(), 1)
    MOI.set(mockoptimizer, MOI.PrimalStatus(1), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(x), 1.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(y), 0.0)
    # test has_values
    @testset "JuMP.has_values" begin
        @test has_values(m)
    end
    # test map_value
    @testset "map_value" begin
        @test map_value(inf, Val(:TransData), 1) == [1., 0.]
        @test map_value(g, Val(:TransData), 1) == 1.
    end
    # test value
    @testset "JuMP.value" begin
        @test value(inf) == [1., 0.]
        @test value(g) == 1.
    end
    # test map_optimizer_index
    @testset "map_optimizer_index" begin
        @test isa(map_optimizer_index(g, Val(:TransData)), MOI.VariableIndex)
        @test isa(map_optimizer_index(inf, Val(:TransData)), Vector{MOI.VariableIndex})
    end
    # test optimizer_index
    @testset "JuMP.optimizer_index" begin
        @test isa(optimizer_index(g), MOI.VariableIndex)
        @test isa(optimizer_index(inf), Vector{MOI.VariableIndex})
    end
    # test dual
    @testset "JuMP.dual" begin
        @test_throws ErrorException dual(g)
    end
end

# Test constraint queries
@testset "Constraint Queries" begin
    # Setup the infinite model
    optimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                         eval_objective_value=false)
    m = InfiniteModel(optimizer)
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_variable(m, inf(par))
    @hold_variable(m, g)
    @constraint(m, c1, g <= 0)
    @constraint(m, c2, inf >= 0)
    # setup the mock solved transcription model
    tm = transcription_model(m)
    @variable(tm, x)
    @variable(tm, y)
    @variable(tm, z)
    @objective(tm, Min, x^2)
    @constraint(tm, c3, x <= 0)
    @constraint(tm, c4, y >= 0)
    @constraint(tm, c5, z >= 0)
    tm.ext[:TransData].hold_to_var[g] = x
    tm.ext[:TransData].infinite_to_vars[inf] = [y, z]
    tm.ext[:TransData].finite_to_constr[c1] = c3
    tm.ext[:TransData].infinite_to_constrs[c2] = [c4, c5]
    modelstring = """
    variables: x, y, z
    minobjective: 1*x*x
    c3: x <= 0.0
    c4: y >= 0.0
    c5: z >= 0.0
    """
    model = MOIU.Model{Float64}()
    MOIU.loadfromstring!(model, modelstring)
    # MOIU.test_models_equal(JuMP.backend(tm).model_cache, model, ["x","y","z"],
    #                        ["c3", "c4", "c5"])
    JuMP.optimize!(tm)

    mockoptimizer = JuMP.backend(tm).optimizer.model
    MOI.set(mockoptimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    MOI.set(mockoptimizer, MOI.ResultCount(), 1)
    MOI.set(mockoptimizer, MOI.PrimalStatus(), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.DualStatus(), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.VariablePrimal(), JuMP.optimizer_index(x), -1.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(), JuMP.optimizer_index(y), 0.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(), JuMP.optimizer_index(z), 1.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(), JuMP.optimizer_index(c3), -1.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(), JuMP.optimizer_index(c4), 0.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(), JuMP.optimizer_index(c5), 1.0)
    # test map_value
    @testset "map_value" begin
        @test map_value(c1, Val(:TransData), 1) == -1.
        @test map_value(c2, Val(:TransData), 1) == [0., 1.]
    end
    # test value
    @testset "JuMP.value" begin
        @test value(c1) == -1.
        @test value(c2) == [0., 1.]
    end
    # test map_optimizer_index
    @testset "map_optimizer_index" begin
        @test isa(map_optimizer_index(c1, Val(:TransData)), MOI.ConstraintIndex)
        @test isa(map_optimizer_index(c2, Val(:TransData)), Vector{<:MOI.ConstraintIndex})
    end
    # test optimizer_index
    @testset "JuMP.optimizer_index" begin
        @test isa(optimizer_index(c1), MOI.ConstraintIndex)
        @test isa(optimizer_index(c2), Vector{<:MOI.ConstraintIndex})
    end
    # test has_values
    @testset "JuMP.has_duals" begin
        @test has_duals(m)
    end
    # test map_dual
    @testset "map_dual" begin
        @test map_dual(c1, Val(:TransData), 1) == -1.
        @test map_dual(c2, Val(:TransData), 1) == [0., 1.]
    end
    # test dual
    @testset "JuMP.dual" begin
        @test dual(c1) == -1.
        @test dual(c2) == [0., 1.]
    end
    # test map_shadow_price
    @testset "map_shadow_price" begin
        @test map_shadow_price(c1, Val(:TransData)) == -1.
        @test map_shadow_price(c2, Val(:TransData)) == [-0., -1.]
    end
    # test shadow_price
    @testset "JuMP.shadow_price" begin
        @test shadow_price(c1) == -1.
        @test shadow_price(c2) == [-0., -1.]
    end
end

# Test LP sensitivity methods
@testset "LP Sensitivities" begin
    # Setup the infinite model
    optimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                         eval_objective_value=false)
    m = InfiniteModel(optimizer)
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_variable(m, inf(par))
    @hold_variable(m, g)
    @constraint(m, c1, g <= 0)
    @constraint(m, c2, inf >= 0)
    # setup the mock solved transcription model
    tm = transcription_model(m)
    @variable(tm, x)
    @variable(tm, y)
    @variable(tm, z)
    @objective(tm, Max, 2x)
    @constraint(tm, c3, x <= 0)
    @constraint(tm, c4, y >= 0)
    @constraint(tm, c5, z >= 0)
    tm.ext[:TransData].hold_to_var[g] = x
    tm.ext[:TransData].infinite_to_vars[inf] = [y, z]
    tm.ext[:TransData].finite_to_constr[c1] = c3
    tm.ext[:TransData].infinite_to_constrs[c2] = [c4, c5]
    modelstring = """
    variables: x, y, z
    maxobjective: 2*x
    c3: x <= 0.0
    c4: y >= 0.0
    c5: z >= 0.0
    """
    model = MOIU.Model{Float64}();
    MOIU.loadfromstring!(model, modelstring);
    # MOIU.test_models_equal(JuMP.backend(tm).model_cache, model, ["x","y","z"],
    #                        ["c3", "c4", "c5"])
    JuMP.optimize!(tm)

    mockoptimizer = JuMP.backend(tm).optimizer.model;
    MOI.set(mockoptimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    MOI.set(mockoptimizer, MOI.ResultCount(), 1)
    MOI.set(mockoptimizer, MOI.PrimalStatus(), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.DualStatus(), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.ObjectiveValue(), 0.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(), JuMP.optimizer_index(x), 0.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(), JuMP.optimizer_index(y), 0.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(), JuMP.optimizer_index(z), 0.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(), JuMP.optimizer_index(c3), -2.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(), JuMP.optimizer_index(c4), 0.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(), JuMP.optimizer_index(c5), 0.0)
    MOI.set(mockoptimizer, MOI.ConstraintBasisStatus(), JuMP.optimizer_index(c3), MOI.NONBASIC)
    MOI.set(mockoptimizer, MOI.ConstraintBasisStatus(), JuMP.optimizer_index(c4), MOI.NONBASIC)
    MOI.set(mockoptimizer, MOI.ConstraintBasisStatus(), JuMP.optimizer_index(c5), MOI.NONBASIC)
    # test map_lp_rhs_perturbation_range
    @testset "map_lp_rhs_perturbation_range" begin
        @test map_lp_rhs_perturbation_range(c1, Val(:TransData), 1e-8) == (-Inf, Inf)
        @test map_lp_rhs_perturbation_range(c2, Val(:TransData), 1e-8) == [(-Inf, Inf),
                                                                           (-Inf, Inf)]
    end
    # test lp_rhs_perturbation_range
    @testset "JuMP.lp_rhs_perturbation_range" begin
        @test lp_rhs_perturbation_range(c1) == (-Inf, Inf)
        @test lp_rhs_perturbation_range(c2) == [(-Inf, Inf), (-Inf, Inf)]
    end
    # test map_lp_objective_perturbation_range
    @testset "map_lp_objective_perturbation_rangee" begin
        @test map_lp_objective_perturbation_range(g, Val(:TransData), 1e-8) == (-2.0, Inf)
        @test map_lp_objective_perturbation_range(inf, Val(:TransData), 1e-8) == [(-Inf, 0.0),
                                                                                (-Inf, 0.0)]
    end
    # test lp_objective_perturbation_range
    @testset "JuMP.lp_objective_perturbation_range" begin
        @test lp_objective_perturbation_range(g) == (-2.0, Inf)
        @test lp_objective_perturbation_range(inf) == [(-Inf, 0.0), (-Inf, 0.0)]
    end
end

ext: jl


text:
## same for all soundingss
nsamples = 100001
nchainspersounding = 5
ppn = 48
@assert mod(ppn, nchainspersounding+1) == 0
zfixed   = [-1e5]
ρfixed   = [1e12]
zstart = 0.0
showgeomplot = false
extendfrac, dz = 1.06, 1.5
nlayers = 50
ρbg = 10.
ntimesperdecade = 10
nfreqsperdecade = 5
useML = false
## make transD options
nmin, nmax = 2, 40
K = transD_GP.GP.OrstUhn()
demean = false
sampledc = true
fbounds = [-1 3]
sdpos = 0.05
sdprop = 0.05
sddc = 0.008
λ, δ = [2], 0.1
save_freq = 100
Tmax = 2.50
nchainsatone = 1
## plot n random soundings and a background response
using Random
nplot = 2
nplot = min(nplot, length(soundings))
opt = transD_GP.SkyTEM1DInversion.makeoperatorandoptions(
                            rseed = 2,
                            soundings,
                            zfixed = zfixed,
                            ρfixed = ρfixed,
                            zstart = zstart,
                            extendfrac = extendfrac,
                            dz = dz,
                            ρbg = ρbg,
                            nlayers = nlayers,
                            ntimesperdecade = ntimesperdecade,
                            nfreqsperdecade = nfreqsperdecade,
                            useML = useML,
                            nmin = nmin,
                            nmax = nmax,
                            K = K,
                            demean = demean,
                            sdpos = sdpos,
                            sdprop = sdprop,
                            sddc = sddc,
                            sampledc = sampledc,
                            fbounds = fbounds,
                            save_freq = save_freq,
                            showgeomplot = showgeomplot,
                            λ = λ,
                            δ = δ
)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
using CSVFiles
using IteratorInterfaceExtensions
using TableTraits
using FileIO
using DataValues
using Test

@testset "CSVFiles" begin

@testset "basic" begin
    array = collect(load(joinpath(@__DIR__, "data.csv")))
    @test length(array) == 3
    @test array == [(Name="John",Age=34.,Children=2),(Name="Sally",Age=54.,Children=1),(Name="Jim",Age=23.,Children=0)]

    output_filename = tempname() * ".csv"

    try
        array |> save(output_filename)

        array2 = collect(load(output_filename))

        @test array == array2
    finally
        GC.gc()
        rm(output_filename)
    end
end

@testset "traits" begin
    csvf = load(joinpath(@__DIR__, "data.csv"))

    @test IteratorInterfaceExtensions.isiterable(csvf) == true
    @test TableTraits.isiterabletable(csvf) == true
    @test TableTraits.supports_get_columns_copy_using_missing(csvf) == true
end

@testset "missing values" begin
    array3 = [(a=DataValue(3),b="df\"e"),(a=DataValue{Int}(),b="something")]

    @testset "default" begin
        output_filename2 = tempname() * ".csv"

        try
            array3 |> save(output_filename2)
        finally
            rm(output_filename2)
        end
    end

    @testset "alternate" begin
        output_filename2 = tempname() * ".csv"

        try
            array3 |> save(output_filename2, nastring="")
        finally
            rm(output_filename2)
        end
    end
end

@testset "Column interface" begin
    csvf2 = load(joinpath(@__DIR__, "data.csv"))
    @test TableTraits.supports_get_columns_copy_using_missing(csvf2) == true
    data = TableTraits.get_columns_copy_using_missing(csvf2)
    @test data == (Name=["John", "Sally", "Jim"], Age=[34.,54.,23.], Children=[2,1,0])
end

@testset "Less Basic" begin
    array = [(Name="John",Age=34.,Children=2),(Name="Sally",Age=54.,Children=1),(Name="Jim",Age=23.,Children=0)]

    @test_broken false # TODO Reenable download test once FileIO is fixed
    # @testset "remote loading" begin
    #     rem_array = collect(load("https://raw.githubusercontent.com/queryverse/CSVFiles.jl/v0.2.0/test/data.csv"))
    #     @test length(rem_array) == 3
    #     @test rem_array == array
    # end

    @testset "can round trip TSV" begin
        output_filename3 = tempname() * ".tsv"

        try
            array |> save(output_filename3)

            array4 = collect(load(output_filename3))
            @test length(array4) == 3
            @test array4 == array
        finally
            GC.gc()
            rm(output_filename3)
        end
    end

    @testset "no quote" begin
        output_filename4 = tempname() * ".csv"

        try
            array |> save(output_filename4, quotechar=nothing)

        finally
            GC.gc()
            rm(output_filename4)
        end
    end
end

@testset "Streams" begin
    data = [(Name="John",Age=34.,Children=2),(Name="Sally",Age=54.,Children=1),(Name="Jim",Age=23.,Children=0)]

    @testset "CSV"  begin
        stream = IOBuffer()
        mark(stream)
        fileiostream = FileIO.Stream{FileIO.format"CSV"}(stream)
        save(fileiostream, data)
        reset(stream)
        mark(stream)
        csvstream = load(fileiostream)
        reloaded_data = collect(csvstream)
        @test IteratorInterfaceExtensions.isiterable(csvstream)
        @test TableTraits.isiterabletable(csvstream)
        @test TableTraits.supports_get_columns_copy_using_missing(csvstream)
        @test reloaded_data == data

        reset(stream)
        csvstream = load(fileiostream)
        reloaded_data2 = TableTraits.get_columns_copy_using_missing(csvstream)
        @test reloaded_data2 == (Name=["John", "Sally", "Jim"], Age=[34., 54., 23.], Children=[2, 1, 0])
    end

    @testset "TSV" begin
        stream = IOBuffer()
        mark(stream)
        fileiostream = FileIO.Stream{FileIO.format"TSV"}(stream)
        save(fileiostream, data)
        reset(stream)
        mark(stream)
        csvstream = load(fileiostream)
        reloaded_data = collect(csvstream)
        @test IteratorInterfaceExtensions.isiterable(csvstream)
        @test TableTraits.isiterabletable(csvstream)
        @test TableTraits.supports_get_columns_copy_using_missing(csvstream)
        @test reloaded_data == data

        reset(stream)
        csvstream = load(fileiostream)
        reloaded_data2 = TableTraits.get_columns_copy_using_missing(csvstream)
        @test reloaded_data2 == (Name=["John", "Sally", "Jim"], Age=[34., 54., 23.], Children=[2, 1, 0])
    end
end

@testset "Compression" begin
    data = [(Name="John",Age=34.,Children=2),(Name="Sally",Age=54.,Children=1),(Name="Jim",Age=23.,Children=0)]

    @testset "CSV" begin
        output_filename = "output.csv.gz"
        try
            save(File{format"CSV"}(output_filename), data)
            reloaded_data = collect(load(File{format"CSV"}(output_filename)))
            @test reloaded_data == data
        finally
            rm(output_filename)
        end
    end

    @testset "TSV" begin
        output_filename = "output.tsv.gz"
        try
            save(File{format"TSV"}(output_filename), data)
            reloaded_data = collect(load(File{format"TSV"}(output_filename)))
            @test reloaded_data == data
        finally
            rm(output_filename)
        end
    end
end

@testset "show" begin
    x = load(joinpath(@__DIR__, "data.csv"))

    @test sprint(show, x) == """
    3x3 CSV file
    Name  │ Age  │ Children
    ──────┼──────┼─────────
    John  │ 34.0 │ 2
    Sally │ 54.0 │ 1
    Jim   │ 23.0 │ 0       """

    @test sprint((stream,data)->show(stream, "text/html", data), x) ==
        "

Name	Age	Children

"John"	34.0	2

"Sally"	54.0	1

"Jim"	23.0	0

"

    @test sprint((stream,data)->show(stream, "application/vnd.dataresource+json", data), x) ==
        "{\"schema\":{\"fields\":[{\"name\":\"Name\",\"type\":\"string\"},{\"name\":\"Age\",\"type\":\"number\"},{\"name\":\"Children\",\"type\":\"integer\"}]},\"data\":[{\"Name\":\"John\",\"Age\":34.0,\"Children\":2},{\"Name\":\"Sally\",\"Age\":54.0,\"Children\":1},{\"Name\":\"Jim\",\"Age\":23.0,\"Children\":0}]}"

    @test showable("text/html", x) == true
    @test showable("application/vnd.dataresource+json", x) == true

    open("data.csv", "r") do f
        x2 = load(Stream{format"CSV"}(f))

        @test sprint(show, x2) == """
    3x3 CSV file
    Name  │ Age  │ Children
    ──────┼──────┼─────────
    John  │ 34.0 │ 2
    Sally │ 54.0 │ 1
    Jim   │ 23.0 │ 0       """

        @test sprint((stream,data)->show(stream, "text/html", data), x2) ==
            "

Name	Age	Children

"John"	34.0	2

"Sally"	54.0	1

"Jim"	23.0	0

"

        @test sprint((stream,data)->show(stream, "application/vnd.dataresource+json", data), x2) ==
            "{\"schema\":{\"fields\":[{\"name\":\"Name\",\"type\":\"string\"},{\"name\":\"Age\",\"type\":\"number\"},{\"name\":\"Children\",\"type\":\"integer\"}]},\"data\":[{\"Name\":\"John\",\"Age\":34.0,\"Children\":2},{\"Name\":\"Sally\",\"Age\":54.0,\"Children\":1},{\"Name\":\"Jim\",\"Age\":23.0,\"Children\":0}]}"

        @test showable("text/html", x2) == true
        @test showable("application/vnd.dataresource+json", x2) == true
    end

end

@testset "savestreaming" begin
    using DataFrames

    df = DataFrame(A = 1:2:1000, B = repeat(1:10, inner=50), C = 1:500)
    df1 = df[1:5, :]
    df2 = df[6:10, :]

    # Test both csv and tsv formats
    for ext in ("csv", "tsv")
        fname = "output.$ext"
        s = savestreaming(fname, df1)
        write(s, df2)
        write(s, df2)   # add this slice twice
        close(s)

        new_df = DataFrame(load(fname))
        @test new_df[1:5,:]   == df1
        @test new_df[6:10,:]  == df2
        @test new_df[11:15,:] == df2

        rm(fname)
    end
end

end # Outer-most testset


text:
# Bit-parallel element nucleotide complementation
@inline function _complement_bitpar(a::A, head::UInt64, tail...) where {A<:NucleicAcidAlphabet{2}}
    return (~head, _complement_bitpar(a, tail...)...)
end

@inline _complement_bitpar(a::A) where {A<:NucleicAcidAlphabet{2}} = ()

@inline function pushfirst(x::Kmer{A,K,N}, nt) where {A,K,N}
    ntbits = UInt64(BioSequences.encode(A(), nt)) << (62 - (64N - 2K))
    #ntbits = UInt64(@inbounds BioSequences.twobitnucs[reinterpret(UInt8, nt) + 0x01]) << (62 - (64N - 2K))
    return Kmer{A,K,N}(_rightshift_carry(2, ntbits, x.data...))
end

@inline function pushlast(x::Kmer{A,K,N}, nt) where {A,K,N}
    ntbits = UInt64(BioSequences.encode(A(), nt))
    #ntbits = UInt64(@inbounds BioSequences.twobitnucs[reinterpret(UInt8, nt) + 0x01])
    _, newbits = _leftshift_carry(2, ntbits, x.data...)
    return Kmer{A,K,N}(newbits)
end

###
### Transformation methods
###

"""
    complement(seq::T) where {T<:Kmer}
Return the complement of a short sequence type `x`.
"""
@inline function BioSequences.complement(seq::T) where {T<:Kmer}
    return T(_complement_bitpar(Alphabet(seq), seq.data...))
end

"""
    reverse(seq::Kmer{A,K,N}) where {A,K,N}
Return the reverse of short sequence type variable `seq`.
"""
@inline function Base.reverse(seq::Kmer{A,K,N}) where {A,K,N}
#    rdata = _reverse(identity, BioSequences.BitsPerSymbol(seq), seq.data...)
    rdata = _reverse(BioSequences.BitsPerSymbol(seq), seq.data...)
    return Kmer{A,K,N}(rightshift_carry(rdata, 64N - 2K))
end

"""
    reverse_complement(x::Kmer)

Return the reverse complement of `x`.
"""
@inline function BioSequences.reverse_complement(x::Kmer{A,K,N}) where {A,K,N}
    return complement(reverse(x))
end

#=
@inline function reverse_complement2(seq::Kmer{A,K,N}) where {A,K,N}
    f = x -> complement_bitpar(x, A())
    rdata = _reverse(f, BioSequences.BitsPerSymbol(seq), seq.data...)
    return Kmer{A,K,N}(rightshift_carry(rdata, 64N - 2K))
end
=#

"""
    canonical(x::Kmer)

Return the canonical sequence of `x`.

A canonical sequence is the numerical lesser of a k-mer and its reverse complement.
This is useful in hashing/counting sequences in data that is not strand specific,
and thus observing the short sequence is equivalent to observing its reverse complement.
"""
@inline canonical(x::Kmer) = min(x, reverse_complement(x))

###
### Old Mer specific specializations of src/biosequence/transformations.jl
### - not currently transferred to new type.

# TODO: Sort this and decide on transferring to new NTuple based kmers or no.

#=
function swap(x::T, i, j) where {T<:AbstractMer}
    i = 2 * length(x) - 2i
    j = 2 * length(x) - 2j
    b = encoded_data(x)
    x = ((b >> i) ⊻ (b >> j)) & encoded_data_type(x)(0x03)
    return T(b ⊻ ((x << i) | (x << j)))
end

function Random.shuffle(x::T) where {T<:AbstractMer}
    # Fisher-Yates shuffle for mers.
    j = lastindex(x)
    for i in firstindex(x):(j - 1)
        j′ = rand(i:j)
        x = swap(x, i, j′)
    end
    return x
end
=#


text:
# AS trees

# plots
using GraphRecipes
using Plots
default(size=(700, 700))

code = :(
function mysum(list)
    out = 0
    for value in list
        out += value
    end
    out
end
)

plot(code, fontsize=12, shorten=0.01, axis_buffer=0.15, nodeshape=:rect, nodecolor=:lightgreen)

Pkg.add("AbstractTrees")

using AbstractTrees

AbstractTrees.children(d::Dict) = [p for p in d]
AbstractTrees.children(p::Pair) = AbstractTrees.children(p[2])
function AbstractTrees.printnode(io::IO, p::Pair)
    str = isempty(AbstractTrees.children(p[2])) ? string(p[1], ": ", p[2]) : string(p[1], ": ")
    print(io, str)
end

d = Dict(:a => 2,:d => Dict(:b => 4,:c => "Hello"),:e => 5.0)

using GraphRecipes
using Plots
default(size=(700, 700))

plot(TreePlot(d), method=:tree, fontsize=10, nodeshape=:ellipse)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
module Classical
using LinearAlgebra
#using StaticArrays

"""
visco()

Implements a power-law dependence for the viscosity ν on the surface density Σ,

    ν = ν₀(Σ/Σ₀)ᵝ

following Grätz, Seiß & Spahn (2018).
"""
function visco(Σrel,β)
    THRESH = 0.01 # Thresholding allows us to use β<0, which would otherwise lead to infinities.
    threshold(Σ) = Σ > THRESH ? Σ : THRESH
    Σr = threshold.(Σrel)
    return Diagonal(Σr.^β)
end

"""
Potentially confusing terminology: a planetary particulate disk is called a "ring",
not a disk.
"""
mutable struct ring
#    α::Float64
#    Δx::Float64
#    Δt::Float64
#    xi, xf ::Float64
    N::Integer
    β::Float64
    ΔtMAX::Float64
    # Defined on nodes:
    X_::Array{Float64,1}
    ΔX_::Array{Float64,1}
    U_::Array{Float64,1}
    # Defined on cells:
    x_::Array{Float64,1}
    σ_::Array{Float64,1}
    Gᵐ::Any # Tridiagonal
#    μ_::Array{Float64,1}
#    νᵐ::Any
end
function ring(N=5;β=0,xf=3.0,xi=0.01)#10*xf/N)
    xtemp_ = collect(range(xi,xf,length=N+1))
    X_ = xtemp_[2:end]
    U_ = 1.0 ./ X_.^4
    ΔX_ = diff(xtemp_)
    x_ = 0.5*(xtemp_[1:end-1]+xtemp_[2:end]) # midpoints
    σ_ = ones(N)
    νᵐ = visco(σ_,β)
    Mᵐ = Diagonal(ΔX_)
    Mⁱᵐ= Diagonal(1.0./ΔX_)
    Kdiag = [1.0 ./ diff(x_)..., 0.0]
    Ksuper = -Kdiag[1:end-1]
    Ksub = zeros(N-1)
    Kᵐ = Tridiagonal(Ksub,Kdiag,Ksuper)
    Dᵐ = Tridiagonal(ones(N-1),-ones(N),zeros(N-1))
    Gᵐ = 3 * Dᵐ * Kᵐ * νᵐ
    μ_ = Mᵐ * σ_ # mass (per azimuthal length) in annulus
    ΔtVN = 0.5 * minimum(diff(X_).^2 ./ diag(νᵐ)[1:end-1] )
    ΔtADV = minimum(ΔX_ ./ U_)
    ΔtMAX = min(ΔtVN,ΔtADV)
    return ring(N,β,ΔtMAX,X_,ΔX_,U_,x_,σ_,Gᵐ)#,μ_)
end

function adjust!(r::ring)
    νᵐ = visco(r.σ_,r.β)
    Kdiag = [1.0 ./ diff(r.x_)..., 0.0]
    Ksuper = -Kdiag[1:end-1]
    Ksub = zeros(r.N-1)
    Kᵐ = Tridiagonal(Ksub,Kdiag,Ksuper)
    Dᵐ = Tridiagonal(ones(r.N-1),-ones(r.N),zeros(r.N-1))
    Gᵐ = 3 * Dᵐ * Kᵐ * νᵐ
    r.Gᵐ = Gᵐ
end

function advect!(r::ring;Δt=r.ΔtMAX)
    σ = r.σ_; U = r.U_;; ΔX = r.ΔX_
    σ[2:end]   += Δt .* (σ[1:end-1] .* U[1:end-1] ./ ΔX[2:end] )
    σ[1:end] -= Δt .* (σ[1:end] .* U[1:end] ./ ΔX[1:end] )
end

function diffuse!(r::ring;Δt=r.ΔtMAX,α=0.5)
    Gᵐ = r.Gᵐ
    N = r.N
    Aᵐ = I(N) -    α  * Δt * Gᵐ
    Bᵐ = I(N) + (1-α) * Δt * Gᵐ
    # Solve:
    r.σ_ = Aᵐ \ (Bᵐ * r.σ_)
#    r.σ_[1]   = 0.0
    r.σ_[end] = 1.0
end

function iterate!(r::ring;Δt=r.ΔtMAX,α=0.5)
    adjust!(r)
    diffuse!(r;Δt=Δt,α=α)
    advect!(r;Δt=Δt)
end

end#module

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
using SPADL

@time test = get_events_data(:wyscout, download = false, path = "wyscout_data")
GC.gc()
# we have to update the loader first
test = get_events(test, "events_Italy.json");

@time events_df = events(test, 2576336);

@time check = convert_to_actions(events_df)

configs = ActionsConfig()
findall(x -> x == "head", configs.bodyparts)

using DataFrames

#### We need to add miliseconds in WyscoutEventFixed
view = []
for i in eachindex(check)
    tmp = (check[i].event_fixed)
    push!(view, tmp)
end
view = DataFrame(view)

# for i in 1:size(view,1)
#     if i > 1 && view[i, :period_id] == view[i-1, :period_id]
#         if view[i, :milliseconds] < view[i-1, :milliseconds]
#             println("we have a problem", i)
#         end
#     end
# end


text:
include("run_scf_and_compare.jl")
include("testcases.jl")

function run_iron_pbe(T; kwargs...)
    # These values were computed using ABINIT with the same kpoints as testcases.jl
    # and Ecut = 20. They are not yet converged and thus require the same discretisation
    # parameters to be obtained.
    ref_evals = [
        [0.0603597727989307, 0.1964963273638626, 0.196496327424440, 0.279192222553112,
         0.2791922225741613, 0.3415221335998876, 0.837882559419754, 0.883850560591423,
         0.8838505606211768, 1.3135367355436536],
        [0.1384929268069029, 0.1847168453364975, 0.223179759800174, 0.320070899985990,
         0.3500724891746176, 0.4685757607370267, 0.541752194212558, 0.751365680734661,
         0.8039132927796911, 1.3939297677405071],
        [-0.017996603976028, 0.2383855826934185, 0.238385582734711, 0.248204676138927,
         0.2509395500598295, 0.2776437400588896, 1.069915401940919, 1.088217176897224,
         1.094997859335961,  1.0949978593466851],
        [0.1102557166995405, 0.2077201723056727, 0.220685303120809, 0.289884460857327,
         0.3490062808992303, 0.3571047250832524, 0.664551132243957, 0.890354172420178,
         0.939822681382406,  1.2259972985258636],
        [0.1723514110126840, 0.1723514110181127, 0.189598224957126, 0.315084007273243,
         0.3150840073174671, 0.5487559496577702, 0.548755949657792, 0.571153866844390,
         1.0611134432316718, 1.1887518709297569],
        [0.1360541296075938, 0.1413608406233668, 0.337616953214017, 0.337616953257584,
         0.3463728840905585, 0.4304010493995122, 0.688627292839765, 0.688627292852315,
         0.885008380770321,  0.9722786718518246],
        [0.0802990962833626, 0.3488798033726516, 0.348879803416372, 0.533263624117060,
         0.560354114948579,  0.5603541149670136, 0.923281827089562, 0.967838872125574,
         0.9678388721641925, 1.300215418446228],
        [0.2341496631160049, 0.2737567834221212, 0.320646675118266, 0.590600827614029,
         0.6440928824646408, 0.6458637753212415, 0.678343515679297, 0.838647690182280,
         0.8763210347583158, 1.4092936521531203],
        [-0.002234753604747, 0.4096246186291687, 0.409624618662776, 0.434260327970128,
         0.5068101375084778, 0.5757957165012942, 1.137207834311533, 1.137826252874365,
         1.170363096833071,  1.170363096849632],
        [0.1518900787487526, 0.3293780680641614, 0.376401550325491, 0.512562269331525,
         0.5557310122303195, 0.6261449425921871, 0.794097184155989, 0.967295197092196,
         1.0000550921659532, 1.2999173820510477],
        [0.2873355363445261, 0.2873355363447599, 0.319313192152575, 0.537629072823137,
         0.5376290728591641, 0.6802062250711767, 0.704199805731151, 0.704199805731498,
         1.1322730987840155, 1.255912074880981],
        [0.2512356397409882, 0.315293666807424,  0.491297439253523, 0.4912974392811193,
         0.5558649368408816, 0.556692128645629,  0.777563890322163, 0.7775638903489546,
         0.9998569230219644, 1.1313796020728688],
    ]
    ref_etot = -18.21465922614397
    ref_magn = 2.98199463

    # Produce reference data and guess for this configuration
    Ecut = 20
    Fe = ElementPsp(iron_bcc.atnum, psp=load_psp("hgh/lda/Fe-q8.hgh"))
    magnetic_moments = [Fe => [4.0]]
    model = model_PBE(Array{T}(iron_bcc.lattice), [Fe => iron_bcc.positions],
                      temperature=0.01, magnetic_moments=magnetic_moments)
    basis = PlaneWaveBasis(model, Ecut; fft_size=[20, 20, 20],
                           kgrid=[4, 4, 4], kshift=[1/2, 1/2, 1/2])

    ρspin = guess_spin_density(basis, magnetic_moments)
    scfres = run_scf_and_compare(T, basis, ref_evals, ref_etot; ρspin=ρspin, kwargs...)

    dVol = scfres.basis.model.unit_cell_volume / prod(scfres.basis.fft_size)
    magnetisation = sum(scfres.ρspin.real) * dVol
    @test magnetisation ≈ ref_magn atol=5e-5
end

if !isdefined(Main, :FAST_TESTS) || !FAST_TESTS
    @testset "Iron PBE (Float64)" begin
        run_iron_pbe(Float64, test_tol=5e-6, scf_tol=1e-8)
    end
end


text:
function train_svdd_model(model::DataType, init_strategy, solver,
                          data::Array{Float64, 2}, labels::Vector{Symbol})
    model = model(data)
    SVDD.initialize!(model, init_strategy)
    set_adjust_K!(model, true)
    SVDD.fit!(model, solver)
    return model
end

function predict_svdd_model(model, test_data)
    return SVDD.classify.(SVDD.predict(model, test_data))
end

function evaluate_with_svdd(model::DataType, init_strategy, solver,
                            data::Array{Float64, 2}, labels::Vector{Symbol},
                            test_data::Array{Float64, 2}, test_labels::Vector{Symbol},
                            quality_metrics)
    try
        time_train = @elapsed model = train_svdd_model(model, init_strategy, solver, data, labels)
        time_pred = @elapsed pred = predict_svdd_model(model, test_data)
        scores = evaluate_prediction(quality_metrics, test_labels, pred)
        gamma = MLKernels.getvalue(model.kernel_fct.alpha)
        C = model.C
        num_support_vectors = length(SVDD.get_support_vectors(model))
        add_evaluation_stats!(scores, time_train, time_pred, gamma, C, num_support_vectors)
        return scores
    catch
        return Dict{Symbol, Any}(:eval_status => string(EvaluationException))
    end
end

function train_pwc_model(init_strategy::Dict{Symbol, Any}, data::Array{Float64, 2}, labels::Vector{Symbol})
    c = KDECache(data, init_strategy[:gamma])
    threshold = minimum(kde(c))
    pwc = PWC(c, FixedThresholdStrategy(threshold))
    calculate_threshold!(pwc, labels)
    return pwc
end

function evaluate_with_pwc(init_strategy::Dict{Symbol, Any},
                           data::Array{Float64, 2}, labels::Vector{Symbol},
                           test_data::Array{Float64, 2}, test_labels::Vector{Symbol},
                           quality_metrics)
    try
        time_train = @elapsed model = train_pwc_model(init_strategy, data, labels)
        time_pred = @elapsed pred = OneClassSampling.predict(model, test_data)
        scores = evaluate_prediction(quality_metrics, test_labels, pred)
        add_evaluation_stats!(scores, time_train, time_pred, init_strategy[:gamma])
        return scores
    catch
        return Dict{Symbol, Any}(:eval_status => string(EvaluationException))
    end
end

function add_evaluation_stats!(scores, time_train, time_pred, gamma, C=nothing, num_support_vectors=nothing)
    scores[:time_train] = time_train
    scores[:time_pred] = time_pred
    scores[:eval_status] = :success
    scores[:gamma] = gamma
    if C !== nothing
        scores[:C] = C
    end
    if num_support_vectors !== nothing
        scores[:num_support_vectors] = num_support_vectors
    end
    return nothing
end

function evaluate_prediction(quality_metrics, labels, pred)
    scores = Dict{Symbol, Any}()
    cm = ConfusionMatrix(pred, labels)
    for (metric_name, metric) in quality_metrics
        m = metric(cm)
        scores[metric_name] = m
    end
    scores[:tp] = cm.tp
    scores[:fp] = cm.fp
    scores[:tn] = cm.tn
    scores[:fn] = cm.fn
    return scores
end

function evaluate_sample_qe(params::Dict{Symbol, Any}, sample_mask::BitArray{1},
                          data::Array{Float64, 2}, labels::Vector{Symbol})
    scores = Dict{Symbol, Any}()
    c = KDECache(data, params[:gamma])
    for t in params[:threshold_strategies]
        _, inlier_mask, outlier_mask = split_masks(t, c, labels)
        dev = OneClassSampling.sample_deviation(c, sample_mask, inlier_mask, outlier_mask)
        scores[Symbol(t)] = dev
    end
    return scores
end

struct EmptySampleException <: Exception end
struct EvaluationException <: Exception end
struct TooLargeSampleException <: Exception
    n::Int
end

function evaluate(sampler::Sampler, models::Dict{Symbol, Any},
                  sample_mask::BitArray{1},
                  test_data::Array{Float64, 2}, test_labels::Vector{Symbol},
                  quality_metrics)
    data, labels = test_data[:, sample_mask], test_labels[sample_mask]
    if length(labels) <= 0
        throw(EmptySampleException())
    end
    scores = Dict{Symbol, Any}()
    for (name, init) in models[:svdd][:init_strategies]
        s = evaluate_with_svdd(eval(models[:svdd][:type]), init, models[:svdd][:solver], data, labels,
                               test_data, test_labels, quality_metrics)
        scores = push!(scores, name => s)
    end
    if :pwc in keys(models)
        s = evaluate_with_pwc(models[:pwc], data, labels,
                              test_data, test_labels, quality_metrics)
        push!(scores, :pwc => s)
    end
    scores_sample_qe = evaluate_sample_qe(models[:sample_qe], sample_mask, test_data, test_labels)
    push!(scores, :sample_qe => scores_sample_qe)
    scores_specials = evaluate_specials(sampler, models, data, labels,
                                        test_data, test_labels, quality_metrics)
    if scores_specials !== nothing
        scores = merge(scores, scores_specials)
    end
    return scores
end

function evaluate_specials(sampler::Sampler, models::Dict{Symbol, Any},
                           data::Array{Float64, 2}, labels::Vector{Symbol},
                           test_data::Array{Float64, 2}, test_labels::Vector{Symbol},
                           quality_metrics)
    return nothing
end

function evaluate_specials(sampler::Union{HDS, RAPID}, models::Dict{Symbol, Any},
                           data::Array{Float64, 2}, labels::Vector{Symbol},
                           test_data::Array{Float64, 2}, test_labels::Vector{Symbol},
                           quality_metrics)
    pwc = deepcopy(models[:pwc])
    pwc[:threshold] = sampler.threshold
    s = evaluate_with_pwc(pwc, data, labels,
                          test_data, test_labels, quality_metrics)
    return Dict{Symbol, Any}(:pwc_retuned => s)
end

max_stars_repo_name: englhardt/ocs-evaluation

max_issues_repo_name: englhardt/ocs-evaluation

max_forks_repo_name: englhardt/ocs-evaluation


text:
include("plotter_setup.jl")

using LambertW
using Statistics

function theoretical_p0(p0)
    x = (1-p0)/exp(1)/p0
    lambertw(x, 0)/lambertw(x, -1)
end

## p0 theoretical only

data_x = (1.:.1:10000)
data_p = theoretical_p0.(data_x)

cla()
figure(figsize=[4, 2.])
plot(data_x, data_p, "-k")
ylim(0, 1)
xlim(.9, 11000)
xscale("log")
xlabel(L"p_0")
ylabel("succ. prob. lower bound")
yticks([0, 0.5, 1.])
savefig("plots/hiding_search/p0_theoretical.pdf", bbox_inches="tight", metadata = Dict("CreationDate" => nothing))

## p0 real scenarios
p0s = .5:.5:2.
p(p0, n) = p0*log(n)/n
repeating = 40

data_n = 100:100:2000

fig, axes = subplots(1, length(p0s), figsize=[6, 2.])
println(p0s)
cla()
for (p0, ax) = zip(p0s, axes)
    data = npzread("hiding_search/sparse_ER_probs_and_bounds_$p0.npz")
    data_real_p = zeros(length(data_n))
    data_real_p_up = zeros(length(data_n))
    data_real_p_down = zeros(length(data_n))
    data_real_bound = zeros(length(data_n))
    data_real_bound_up = zeros(length(data_n))
    data_real_bound_down = zeros(length(data_n))
    for m =1:length(data_n)
        data_b = data[m,:,1]
        data_p = data[m,:,2]

        data_real_bound[m] = mean(data_b)
        data_real_bound_up[m] = maximum(data_b)
        data_real_bound_down[m] = minimum(data_b)

        data_real_p[m] = mean(data_p)
        data_real_p_up[m] = maximum(data_p)
        data_real_p_down[m] = minimum(data_p)
    end

    if p0 > 1
        ax.plot(data_n, data_real_bound, "-b", markersize=.8, linewidth=.8)
        ax.fill_between(data_n, data_real_bound_up, data_real_bound_down, color="b", alpha=.2)
    end
    ax.plot(data_n, data_real_p, ":r", markersize=.8, linewidth=.8)
    ax.fill_between(data_n, data_real_p_up, data_real_p_down, color="r", alpha=.2)

    if p0 > 1
        ax.plot([minimum(data_n), maximum(data_n)], fill(theoretical_p0(p0), 2), "--k", linewidth=.8)
    end
    ax.set_title(latexstring("p_0=$p0"))
    ax.set(xlabel=L"n", ylim=[0,1], yticks=0:0.5:1.)
    ax.label_outer()
end
axes[1].set(ylabel="succ. prob.")
savefig("plots/hiding_search/real_statistics_p0.pdf", bbox_inches="tight", metadata = Dict("CreationDate" => nothing))


text:
using Revise
using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using DelimitedFiles
using BenchmarkTools
using Random
using LinearAlgebra
using GLM
using Plots

#simulat data with k true predictors, from distribution d and with link l.
n = 1000
p = 10000
k = 10
d = NegativeBinomial
# l = canonicallink(d())
l = LogLink()

#set random seed
Random.seed!(2019)

#construct SnpArraym, snpmatrix, and non genetic covariate (intercept)
x = simulate_random_snparray(n, p, "test1.bed")
xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
z = ones(n, 1)

# simulate response, true model b, and the correct non-0 positions of b
true_b = zeros(p)
# true_b[1:5] = [0.01; 0.05; 0.1; 0.25; 0.5]
true_b[1:k] .= collect(0.1:0.1:1.0)
true_c = [4.0]
# true_b[1:k] = rand(Normal(0, 0.3), k)
shuffle!(true_b)
correct_position = findall(!iszero, true_b)

#simulate phenotypes (e.g. vector y)
if d == Normal || d == Bernoulli || d == Poisson
    prob = linkinv.(l, xbm * true_b)
    clamp!(prob, -20, 20)
    y = [rand(d(i)) for i in prob]
    # prob = linkinv.(l, xbm * true_b + z * true_c)
    # clamp!(prob, -20, 20)
    # y = [rand(d(i)) for i in prob]
    # k = k + 1
elseif d == NegativeBinomial
    nn = 1
    μ = linkinv.(l, xbm * true_b)
    # μ = linkinv.(l, xbm * true_b + z * true_c)
    # k = k + 1
    clamp!(μ, -20, 20)
    prob = 1 ./ (1 .+ μ ./ nn)
    y = [rand(d(nn, i)) for i in prob] #number of failtures before nn success occurs
elseif d == Gamma
    μ = linkinv.(l, xbm * true_b)
    β = 1 ./ μ # here β is the rate parameter for gamma distribution
    y = [rand(d(α, i)) for i in β] # α is the shape parameter for gamma
end
y = Float64.(y)
# histogram(y, bins=30)
# mean(y)
# var(y)

#run IHT
result = L0_reg(x, xbm, z, y, 1, k, d(), l, debias=false)

#check result
compare_model = DataFrame(
    position    = correct_position,
    true_β      = true_b[correct_position],
    estimated_β = result.beta[correct_position])
println("Total iteration number was " * string(result.iter))
println("Total time was " * string(result.time))
println("Total found predictors = " * string(length(findall(!iszero, result.beta[correct_position]))))

#clean up
rm("test1.bed", force=true)

# SIMULATION FOR FLOATING POINT MATRICES
using Revise
using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using DelimitedFiles
using BenchmarkTools
using Random
using LinearAlgebra
using GLM
using Plots

#simulat data with k true predictors, from distribution d and with link l.
n = 5000
p = 30000
k = 10
d = Normal
l = canonicallink(d())
# l = LogLink()

#set random seed
Random.seed!(1111)

#construct x matrix and non genetic covariate (intercept)
T = Float64
x = randn(T, n, p)
z = ones(T, n, 1)

# simulate response, true model b, and the correct non-0 positions of b
true_b = zeros(T, p)
true_b[1:k] .= collect(0.1:0.1:1.0)
true_c = [4.0]
shuffle!(true_b)
correct_position = findall(!iszero, true_b)

#simulate phenotypes (e.g. vector y)
if d == Normal || d == Bernoulli || d == Poisson
    prob = linkinv.(l, x * true_b)
    clamp!(prob, -20, 20)
    y = [rand(d(i)) for i in prob]
    # prob = linkinv.(l, x * true_b + z * true_c)
    # clamp!(prob, -20, 20)
    # y = [rand(d(i)) for i in prob]
    # k = k + 1
elseif d == NegativeBinomial
    nn = 10
    μ = linkinv.(l, x * true_b)
    # μ = linkinv.(l, x * true_b + z * true_c)
    # k = k + 1
    clamp!(μ, -20, 20)
    prob = 1 ./ (1 .+ μ ./ nn)
    y = [rand(d(nn, Float64(i))) for i in prob] #number of failtures before nn success occurs
elseif d == Gamma
    μ = linkinv.(l, x * true_b)
    β = 1 ./ μ # here β is the rate parameter for gamma distribution
    y = [rand(d(α, i)) for i in β] # α is the shape parameter for gamma
end
y = T.(y)
# histogram(y, bins=30)
mean(y)
var(y)

#run IHT
result = L0_reg(x, z, y, 1, k, d(), l, debias=false)
@benchmark result = L0_reg(x, z, y, 1, k, d(), l, debias=false)

#check result
compare_model = DataFrame(
    position    = correct_position,
    true_β      = true_b[correct_position],
    estimated_β = result.beta[correct_position])
println("Total iteration number was " * string(result.iter))
println("Total time was " * string(result.time))
println("Total found predictors = " * string(length(findall(!iszero, result.beta[correct_position]))))

#BELOW ARE SIMULATION FOR binomial
using Revise
using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using BenchmarkTools
using Random
using LinearAlgebra
using GLM

#simulat data with k true predictors, from distribution d and with link l.
n = 1000
p = 10000
k = 10
d = Binomial
l = CloglogLink()
nn = 10 #number of tries for binomial/negative-binomial

#set random seed
Random.seed!(2019)

#construct snpmatrix, covariate files, and true model b
x = simulate_random_snparray(n, p, "tmp.bed")
xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
z = ones(n, 1) # the intercept
true_b = zeros(p)
true_b[1:k] = randn(k)
shuffle!(true_b)
correct_position = findall(x -> x != 0, true_b)

#simulate phenotypes (e.g. vector y)
μ = linkinv.(l, xbm * true_b)
prob = μ ./ nn
y = [rand(d(nn, i)) for i in prob]
y = Float64.(y)

#run IHT
result = L0_reg(x, xbm, z, y, 1, k, d(), l, debias=false, init=false, show_info=false, convg=true)
# @benchmark L0_reg(x, xbm, z, y, 1, k, d(), l, debias=false, init=false, show_info=false, convg=true) seconds = 60

#check result
compare_model = DataFrame(
    true_β      = true_b[correct_position],
    estimated_β = result.beta[correct_position])
println("Total iteration number was " * string(result.iter))
println("Total time was " * string(result.time))
println("Total found predictors = " * string(length(findall(!iszero, result.beta[correct_position]))))

#clean up
rm("tmp.bed", force=true)

#BELOW ARE SIMULATION FOR gamma
using Revise
using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using BenchmarkTools
using Random
using LinearAlgebra
using GLM

#simulat data with k true predictors, from distribution d and with link l.
n = 1000
p = 10000
k = 10
d = Gamma
# l = canonicallink(d())
l = LogLink()
θ = 3 #scale parameter for gamma

#set random seed
Random.seed!(2019)

#construct snpmatrix, covariate files, and true model b
x = simulate_random_snparray(n, p, "tmp.bed")
xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
z = ones(n, 1) # the intercept
true_b = zeros(p)
# true_b[1:k] = randn(k)
true_b[1:k] = rand(Normal(0, 1.0), k)
shuffle!(true_b)
correct_position = findall(x -> x != 0, true_b)

#simulate phenotypes (e.g. vector y)
μ = linkinv.(l, xbm * true_b)
clamp!(μ, -20, 20)
y = [rand(d(i, θ)) for i in μ]
histogram(y)

#run IHT
result = L0_reg(x, xbm, z, y, 1, k, d(), l, debias=false, init=false, show_info=false)
# @benchmark L0_reg(x, xbm, z, y, 1, k, d(), l, debias=false, init=false, show_info=false, convg=true) seconds = 60

#check result
compare_model = DataFrame(
    true_β      = true_b[correct_position],
    estimated_β = result.beta[correct_position])
println("Total iteration number was " * string(result.iter))
println("Total time was " * string(result.time))
println("Total found predictors = " * string(length(findall(!iszero, result.beta[correct_position]))))

#clean up
rm("tmp.bed", force=true)

#BELOW ARE SIMULATION FOR inverse Gaussian
using Revise
using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using BenchmarkTools
using Random
using LinearAlgebra
using GLM

#simulat data with k true predictors, from distribution d and with link l.
n = 1000
p = 10000
k = 10
d = InverseGaussian
l = LogLink()
λ = 1 # shape parameter for inverse gaussian

#set random seed
Random.seed!(2019)

#construct snpmatrix, covariate files, and true model b
x = simulate_random_snparray(n, p, "tmp.bed")
xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
z = ones(n, 1) # the intercept
true_b = zeros(p)
true_b[1:k] = randn(k)
# true_b[1:k] = rand(Normal(0, 0.3), k)
shuffle!(true_b)
correct_position = findall(x -> x != 0, true_b)

#simulate phenotypes (e.g. vector y)
μ = linkinv.(l, xbm * true_b)
clamp!(μ, -20, 20)
mean_parameter = 1 ./ μ #mean parameter for inverse gaussian distribution
y = [rand(d(i, λ)) for i in mean_parameter]
clamp!(y, 0, 20)

#run IHT
result = L0_reg(x, xbm, z, y, 1, k, d(), l, debias=true, init=false, show_info=false)
# @benchmark L0_reg(x, xbm, z, y, 1, k, d(), l, debias=false, init=false, show_info=false, convg=true) seconds = 60

#check result
compare_model = DataFrame(
    true_β      = true_b[correct_position],
    estimated_β = result.beta[correct_position])
println("Total iteration number was " * string(result.iter))
println("Total time was " * string(result.time))
println("Total found predictors = " * string(length(findall(!iszero, result.beta[correct_position]))))

#clean up
rm("tmp.bed", force=true)

############## CROSS VALIDATION SIMULATION for normal, bernoulli, and Poisson

#first add workers
using Distributed
addprocs(4)
nprocs()

#load packages into all worker
using Revise
using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using BenchmarkTools
using Random
using LinearAlgebra
using BenchmarkTools
using GLM

#simulat data
n = 1000
p = 10000
k = 10 # number of true predictors
d = Normal
l = canonicallink(d())

#set random seed
Random.seed!(2018)

#construct snpmatrix, covariate files, and true model b
x = simulate_random_snparray(n, p, "tmp.bed")
xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
z = ones(n, 1) # the intercept
true_b = zeros(p)
d == Poisson ? true_b[1:k] = rand(Normal(0, 0.3), k) : true_b[1:k] = randn(k)
shuffle!(true_b)
correct_position = findall(x -> x != 0, true_b)

#simulate phenotypes (e.g. vector y)
y_temp = xbm * true_b
prob = linkinv.(l, y_temp)
y = [rand(d(i)) for i in prob]
y = Float64.(y)

#specify path and folds
path = collect(1:20)
num_folds = 4
folds = rand(1:num_folds, size(x, 1))

#compute cross validation
mses = cv_iht(d(), l, x, z, y, 1, path, num_folds, folds=folds, use_maf=false, debias=true, parallel=true);

#compute l0 result using best estimate for k
k_est = argmin(mses)
result = L0_reg(x, xbm, z, y, 1, k_est, d(), l, debias=false, init=false, show_info=false)

#check result
compare_model = DataFrame(
    true_β      = true_b[correct_position],
    estimated_β = result.beta[correct_position])
println("Total iteration number was " * string(result.iter))
println("Total time was " * string(result.time))
println("Total found predictors = " * string(length(findall(!iszero, result.beta[correct_position]))))

#clean up
rm("tmp.bed", force=true)

############## RUNNING A BUNCH OF MODELS for normal, bernoulli, and Poisson

#first add workers
using Distributed
addprocs(4)
nprocs()

#load packages into all worker
using Revise
using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using BenchmarkTools
using Random
using LinearAlgebra
using BenchmarkTools
using GLM

#simulat data
n = 2000
p = 10000
k = 10 # number of true predictors
d = Normal
l = canonicallink(d())

#set random seed
Random.seed!(2019)

#construct snpmatrix, covariate files, and true model b
x = simulate_random_snparray(n, p, "tmp.bed")
xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
z = ones(n, 1) # the intercept
true_b = zeros(p)
d == Poisson ? true_b[1:k] = rand(Normal(0, 1.0), k) : true_b[1:k] = randn(k)
shuffle!(true_b)
correct_position = findall(x -> x != 0, true_b)

#simulate phenotypes (e.g. vector y)
y_temp = xbm * true_b
prob = linkinv.(l, y_temp)
y = [rand(d(i)) for i in prob]
y = Float64.(y)

#specify path (i.e. all models `k` you want to test)
path = collect(1:20)

#run results
result = iht_run_many_models(d(), l, x, z, y, 1, path, parallel=true, debias=true);
@benchmark iht_run_many_models(d(), l, x, z, y, 1, path, parallel=true, debias=true) seconds=60

#clean up
rm("tmp.bed", force=true)

#RUN IHT with all CPU available
using Distributed
addprocs(8)
nprocs()

using Revise
using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using BenchmarkTools
using Random
using LinearAlgebra
using GLM

#simulat data with k true predictors, from distribution d and with link l.
n = 1000
p = 10000
k = 10
d = NegativeBinomial
# l = canonicallink(d())
l = LogLink()

#set random seed
Random.seed!(33)

#construct snpmatrix, covariate files, and true model b
T = Float32
x = simulate_random_snparray(n, p, undef)
xbm = SnpBitMatrix{T}(x, model=ADDITIVE_MODEL, center=true, scale=true);
z = ones(T, n, 1) # the intercept

# simulate response, true model b, and the correct non-0 positions of b
true_b = zeros(T, p)
# true_b[1:4] .= [0.1; 0.25; 0.5; 0.8]
true_b[1:10] .= collect(0.1:0.1:1.0)
# true_b[1:k] = rand(Normal(0, 0.3), k)
shuffle!(true_b)
correct_position = findall(!iszero, true_b)

#simulate phenotypes (e.g. vector y)
if d == Normal || d == Poisson || d == Bernoulli
    prob = linkinv.(l, xbm * true_b)
    clamp!(prob, -20, 20)
    y = [rand(d(i)) for i in prob]
elseif d == NegativeBinomial
    nn = 10
    μ = linkinv.(l, xbm * true_b)
    clamp!(μ, -20, 20)
    prob = 1 ./ (1 .+ μ ./ nn)
    y = [rand(d(nn, Float64(i))) for i in prob] #number of failtures before nn success occurs
elseif d == Gamma
    μ = linkinv.(l, xbm * true_b)
    β = 1 ./ μ # here β is the rate parameter for gamma distribution
    y = [rand(d(α, i)) for i in β] # α is the shape parameter for gamma
end
y = T.(y)
# histogram(y)
# var(y) / mean(y)

#specify path and folds
path = collect(1:20)
num_folds = 8
folds = rand(1:num_folds, size(x, 1))

# run threaded IHT
# result = iht_run_many_models(d(), l, x, z, y, 1, path);
destin = "/Users/biona001/Desktop/test/"
mses = cv_iht(d(), l, x, z, y, 1, path, num_folds, destin=destin, folds=folds, init=false, use_maf=false, debias=false, parallel=true)
mses = cv_iht_distribute_fold(d(), l, x, z, y, 1, path, num_folds, destin=destin, folds=folds, init=false, use_maf=false, debias=false, parallel=true)

#benchmarking
@benchmark cv_iht(d(), l, x, z, y, 1, path, num_folds, folds=folds, init=false, use_maf=false, debias=false, parallel=false) seconds=30 #33.800s, 135.13MiB
@benchmark cv_iht(d(), l, x, z, y, 1, path, num_folds, folds=folds, init=false, use_maf=false, debias=false, parallel=true) seconds=30  #6.055 s, 34.23 MiB

#run IHT
result = L0_reg(x, xbm, z, y, 1, argmin(mses), d(), l, debias=true, init=false, use_maf=false)
# @benchmark L0_reg(x, xbm, z, y, 1, k, d(), l, debias=false, init=false, show_info=false) seconds=60
# @code_warntype L0_reg(x, xbm, z, y, 1, k, d(), l, debias=false, init=false, show_info=false)

#check result
compare_model = DataFrame(
    position    = correct_position,
    true_β      = true_b[correct_position],
    estimated_β = result.beta[correct_position])
println("Total iteration number was " * string(result.iter))
println("Total time was " * string(result.time))
println("Total found predictors = " * string(length(findall(!iszero, result.beta[correct_position]))))

#clean up
rm("tmp.bed", force=true)

# CROSS VALIDATION ON GENERAL MATRIX
using Distributed
addprocs(8)
nprocs()

# @everywhere using Revise
using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using BenchmarkTools
using Random
using LinearAlgebra
using GLM

#simulat data with k true predictors, from distribution d and with link l.
n = 2000
p = 10000
k = 10
d = Normal
l = canonicallink(d())
# l = LogLink()

#set random seed
Random.seed!(2019)

#construct x matrix and non genetic covariate (intercept)
T = Float64
x = randn(T, n, p)
z = ones(T, n, 1)

# simulate response, true model b, and the correct non-0 positions of b
true_b = zeros(T, p)
true_b[1:k] .= collect(0.1:0.1:1.0)
true_c = [T.(4.0)]
shuffle!(true_b)
correct_position = findall(!iszero, true_b)

#simulate phenotypes (e.g. vector y)
if d == Normal || d == Bernoulli || d == Poisson
    prob = linkinv.(l, x * true_b)
    clamp!(prob, -20, 20)
    y = [rand(d(i)) for i in prob]
    # prob = linkinv.(l, x * true_b + z * true_c)
    # clamp!(prob, -20, 20)
    # y = [rand(d(i)) for i in prob]
    # k = k + 1
elseif d == NegativeBinomial
    nn = 10
    μ = linkinv.(l, x * true_b)
    # μ = linkinv.(l, x * true_b + z * true_c)
    # k = k + 1
    clamp!(μ, -20, 20)
    prob = 1 ./ (1 .+ μ ./ nn)
    y = [rand(d(nn, Float64(i))) for i in prob] #number of failtures before nn success occurs
elseif d == Gamma
    μ = linkinv.(l, x * true_b)
    β = 1 ./ μ # here β is the rate parameter for gamma distribution
    y = [rand(d(α, i)) for i in β] # α is the shape parameter for gamma
end
y = T.(y)
mean(y)
var(y)

#specify path and folds
path = collect(1:50)
num_folds = 8
folds = rand(1:num_folds, size(x, 1))

# run threaded IHT
# result = iht_run_many_models(d(), l, x, z, y, 1, path);
mses = cv_iht(d(), l, x, z, y, 1, path, num_folds, folds=folds, init=false, use_maf=false, debias=false, parallel=true)
mses = cv_iht_distribute_fold(d(), l, x, z, y, 1, path, num_folds, folds=folds, init=false, use_maf=false, debias=false, parallel=true)

@benchmark cv_iht(d(), l, x, z, y, 1, path, num_folds, folds=folds, init=false, use_maf=false, debias=false, parallel=false) seconds=30#14.409s, 889.53 MiB
@benchmark cv_iht(d(), l, x, z, y, 1, path, num_folds, folds=folds, init=false, use_maf=false, debias=false, parallel=true) seconds=30 #19.509s, 465.59 MiB

#run IHT
result = L0_reg(x, z, y, 1, argmin(mses), d(), l, debias=true, init=false, use_maf=false)
# @benchmark L0_reg(x, xbm, z, y, 1, k, d(), l, debias=false, init=false, show_info=false) seconds=60
# @code_warntype L0_reg(x, xbm, z, y, 1, k, d(), l, debias=false, init=false, show_info=false)

#check result
compare_model = DataFrame(
    position    = correct_position,
    true_β      = true_b[correct_position],
    estimated_β = result.beta[correct_position])
println("Total iteration number was " * string(result.iter))
println("Total time was " * string(result.time))
println("Total found predictors = " * string(length(findall(!iszero, result.beta[correct_position]))))

using Revise
using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using Random
using LinearAlgebra
using GLM
using CSV
using BenchmarkTools

Random.seed!(1201)
x = SnpArray(undef, 1000, 1000)
@time naive_impute(x, "hi.bed")


text:
using Unzip
using Documenter: deploydocs, makedocs

makedocs(sitename = "Unzip.jl", modules = [Unzip], doctest = false)
deploydocs(repo = "github.com/bramtayl/Unzip.jl.git")


text:
using Base.Test

@testset "Simple IntListSet Test" begin
  s = IntListSet(5)
  @show s
  push!(s, 5)
  @show s
  @test 5 in s
  @test 6 ∉ s
  @test_throws ArgumentError push!(s, 6)
  @test_throws ArgumentError push!(s, 0)
  @test_throws ArgumentError push!(s, -1)
  push!(s, 1)
  push!(s, 2)
  push!(s, 5)
  @test length(s) == 3
  @test setdiff(collect(s),[1,2,5]) == []
  delete!(s, 1)
  @test setdiff(collect(s),[2,5]) == []
  push!(s, 1)
  push!(s, 3)
  intersect!(s, [1,2])
  @test length(s) == 2
  intersect!(s, [1,3])
  @test length(s) == 1
  for i=1:5
    push!(s, i)
  end
  @test length(s) == 5
end


text:
@testset "ArrayTabular" begin
    @testset "ArraySeries" begin
        @test @inferred(Series([2,3,4,5])) isa Series

        s = Series([2,3,4,5])

        @test @inferred(indices(s)) === (Base.OneTo(4),)

        @test @inferred(s[3]) == 4
        @test_throws Exception s[1,2]
        @test @inferred(s[[3,4]]) == Series(Dict(3=>4, 4=>5))
        @test @inferred(s[:])::Series == s

        # @test @inferred(Series((2,3,4,5))) isa ArraySeries
    end

    @testset "ArrayTable" begin
        # Flat ArrayTable
        @test @inferred(Table([1 4; 2 5; 3 6])) isa Table{Matrix{Int}}
        t2 = Table([1 4; 2 5; 3 6])
        @test @inferred(indices(t2)) === (Base.OneTo(3), Base.OneTo(2))
        @test @inferred(t2[2,1]) === 2
        @test @inferred(t2[2,[1,2]]) == Series([2,5])
        @test @inferred(t2[2,:]) == Series([2,5])
        @test @inferred(t2[[1,2,3],1]) == Series([1,2,3])
        @test @inferred(t2[[1,2,3],[1,2]]) == t2
        @test @inferred(t2[[1,2,3],:]) == t2
        @test @inferred(t2[:,1]) == Series([1,2,3])
        @test @inferred(t2[:,[1,2]]) == t2
        @test @inferred(t2[:,:]) == t2

        # Nested ArrayTable
        @test @inferred(Table([[1,2,3],[4,5,6]])) isa Table{<:Vector}
        t3 = Table([[1,2,3],[4,5,6]])
        @test @inferred(indices(t3)) === (Base.OneTo(3), Base.OneTo(2))
        @test @inferred(t3[2,1]) === 2
        @test @inferred(t3[2,[1,2]]) == Series([2,5])
        @test @inferred(t3[:,1]) == Series([1,2,3])
        @test @inferred(t3[2,:]) == Series([2,5])
        @test @inferred(t3[:,[1,2]]) == t3
        @test @inferred(t3[:,:]) == t3

        # With other inner types
        # Dict inner storage
        @test @inferred(Table([Dict(:a=>1, :b=>2, :c=>3), Dict(:a=>4, :b=>5, :c=>6)])) isa Table{<:Vector}
        t4 = Table([Dict(:a=>1, :b=>2, :c=>3), Dict(:a=>4, :b=>5, :c=>6)])
        @test @inferred(indices(t4))[1] isa Base.KeyIterator # do a more specific test?
        @test @inferred(indices(t4))[2] === Base.OneTo(2)
        @test @inferred(t4[:b,1]) === 2
        @test @inferred(t4[:,1]) == Series(Dict(:a=>1, :b=>2, :c=>3))
        @test @inferred(t4[:b,:]) == Series([2,5])
        @test @inferred(t4[:,:]) == t4

        # TupleTabular inner type
        @test @inferred(Table([(l"a"=>1, l"b"=>2.0, l"c"=>3f0), (l"a"=>4, l"b"=>5.0, l"c"=>6f0)])) isa Table{<:Vector}
        t5 = Table([(l"a"=>1, l"b"=>2.0, l"c"=>3f0), (l"a"=>4, l"b"=>5.0, l"c"=>6f0)])
        @test @inferred(indices(t5)) === ((l"a", l"b", l"c"), Base.OneTo(2))
        @test @inferred(t5[l"b",1]) === 2.0
        @test @inferred(t5[:,1]) == Series(l"a"=>1, l"b"=>2.0, l"c"=>3f0)
        @test @inferred(t5[l"b",:]) == Series([2.0,5.0])
        @test @inferred(t5[:,:]) == t5

        # StructTabular inner type
        @test @inferred(Table([1+2im, 3+4im])) isa Table{Vector{Complex{Int}}}
        t6 = Table([1+2im, 3+4im])
        @test @inferred(indices(t6)) === ((l"re", l"im"), Base.OneTo(2))
        @test @inferred(t6[l"im",1]) === 2
        @test @inferred(t6[:,1]) == Series(1+2im)
        if VERSION <= v"0.6.0-rc3.0"
            # Inference works at REPL but not in test
            @test t6[l"im",:] == Series([2, 4])
            #@test_broken @inferred(t6[l"im",:]) == Series([2, 4])
        else
            @test @inferred(t6[l"im",:]) == Series([2, 4])
        end

        if VERSION <= v"0.6.0-rc3.0"
            # Inference works at REPL but not in test
            @test t6[:,:] == t6
            #@test_broken @inferred(t6[:,:]) == t6
        else
            @test @inferred(t6[:,:]) == t6
        end

        # # Tuples become arrays...
        # @test @inferred(Table(([1,2,3],[4,5,6]))) isa ArrayTable{<:Vector}
        # @test @inferred(Table(((1,2,3),(4,5,6)))) isa ArrayTable{<:Vector{<:ArraySeries}}
    end
end


text:
# mathematical operators between two Timestamp objects
for op in [:(Base.:+), :(Base.:-), :(Base.:*), :(Base.:/)]
    @eval begin
        function ($op)(ts1::Timestamp, ts2::Timestamp)
            vals   = ($op).(ts1.val, ts2.val)
            Timestamp(ts1.ts, vals)
        end # function
    end # eval
end # loop

# mathematical operators between a Timestamp object and Number
for op in [:(Base.:+), :(Base.:-), :(Base.:*), :(Base.:/)]
    @eval begin
        function ($op)(ts1::Timestamp, num::Number)
            vals = ($op).(ts1.val, num)
            Timestamp(ts1.ts, vals)
        end # function
    end # eval
end # loop

# comparison operators between two Timestamp objects
#for op in [:>, :<, :==, :>=, :<=]
for op in [:(Base.:>), :(Base.:<), :(Base.:>=), :(Base.:<=), :(Base.isless)]
    @eval begin
        function ($op)(ts1::Timestamp{D, T1}, ts2::Timestamp{D, T2}) where {D, T1, T2}
            ($op)(ts1.ts, ts2.ts)
        end # function
    end # eval
end # loop


text:
export
ExpectationPropagationRule,
@expectationPropagationRule

"""
A non-specific expectation propagation update
"""
abstract type ExpectationPropagationRule{factor_type} <: MessageUpdateRule end

messagePassingSchedule(variable::Variable) = messagePassingSchedule([variable])

function inferUpdateRule!(entry::ScheduleEntry,
                          rule_type::Type{T},
                          inferred_outbound_types::Dict{Interface, <:Type}
                         ) where T<:ExpectationPropagationRule
    # Collect inbound types
    inbound_types = collectInboundTypes(entry, rule_type, inferred_outbound_types)

    # Find outbound id
    outbound_id = something(findfirst(isequal(entry.interface), entry.interface.node.interfaces), 0)

    # Find applicable rule(s)
    applicable_rules = Type[]
    for rule in leaftypes(entry.message_update_rule)
        if isApplicable(rule, inbound_types, outbound_id)
            push!(applicable_rules, rule)
        end
    end

    # Select and set applicable rule
    if isempty(applicable_rules)
        error("No applicable $(rule_type) update for $(typeof(entry.interface.node)) node with inbound types: $(join(inbound_types, ", "))")
    elseif length(applicable_rules) > 1
        error("Multiple applicable $(rule_type) updates for $(typeof(entry.interface.node)) node with inbound types: $(join(inbound_types, ", ")): $(join(applicable_rules, ", "))")
    else
        entry.message_update_rule = first(applicable_rules)
    end

    return entry
end

function collectInboundTypes(entry::ScheduleEntry,
                             ::Type{T},
                             inferred_outbound_types::Dict{Interface, <:Type}
                            ) where T<:ExpectationPropagationRule
    inbound_message_types = Type[]
    for node_interface in entry.interface.node.interfaces
        if isDeltaConstrained(node_interface.partner)
            push!(inbound_message_types, Message{PointMass})
        else
            push!(inbound_message_types, inferred_outbound_types[node_interface.partner])
        end
    end

    return inbound_message_types
end

"""
`@expectationPropagationRule` registers a expectation propagation update
rule by defining the rule type and the corresponding methods for the `outboundType`
and `isApplicable` functions. If no name (type) for the new rule is passed, a
unique name (type) will be generated. Returns the rule type.
"""
macro expectationPropagationRule(fields...)
    # Init required fields in macro scope
    node_type = :unknown
    outbound_type = :unknown
    inbound_types = :unknown
    outbound_id = :unknown
    name = :auto # Triggers automatic naming unless overwritten

    # Loop over fields because order is unknown
    for arg in fields
        (arg.args[1] == :(=>)) || error("Invalid call to @expectationPropagationRule")

        if arg.args[2].value == :node_type
            node_type = arg.args[3]
        elseif arg.args[2].value == :outbound_type
            outbound_type = arg.args[3]
            (outbound_type.head == :curly && outbound_type.args[1] == :Message) || error("Outbound type for ExpectationPropagationRule should be a Message")
        elseif arg.args[2].value == :inbound_types
            inbound_types = arg.args[3]
            (inbound_types.head == :tuple) || error("Inbound types should be passed as Tuple")
        elseif arg.args[2].value == :outbound_id
            outbound_id = arg.args[3]
        elseif arg.args[2].value == :name
            name = arg.args[3]
        else
            error("Unrecognized field $(arg.args[2].value) in call to @expectationPropagationRule")
        end
    end

    # Assign unique name if not set already
    if name == :auto
        # Added hash ensures that the rule name is unique
        msg_types_hash = string(hash(vcat([outbound_type], inbound_types)))[1:6]
        name = Symbol("EP$(node_type)$(msg_types_hash)")
    end

    # Build validators for isApplicable
    input_type_validators = Expr[]

    push!(input_type_validators, :(length(input_types) == $(length(inbound_types.args))))
    for (i, i_type) in enumerate(inbound_types.args)
        if i_type != :Nothing
            # Only validate inbounds required for message update
            push!(input_type_validators, :(ForneyLab.matches(input_types[$i], $i_type)))
        end
    end

    expr = quote
        struct $name <: ExpectationPropagationRule{$node_type} end
        ForneyLab.outboundType(::Type{$name}) = $outbound_type
        ForneyLab.isApplicable(::Type{$name}, input_types::Vector{<:Type}, outbound_id::Int64) = begin
            $(reduce((current, item) -> :($current && $item), input_type_validators, init = :(outbound_id === $outbound_id)))
        end
    end

    return esc(expr)
end

"""
Find the inbound types that are required to compute the message for `entry`.
Returns a vector with inbound types that correspond with required interfaces.
"""
function collectInbounds(entry::ScheduleEntry, ::Type{T}) where T<:ExpectationPropagationRule
    interface_to_schedule_entry = current_inference_algorithm.interface_to_schedule_entry

    inbounds = Any[]
    for node_interface in entry.interface.node.interfaces
        inbound_interface = ultimatePartner(node_interface)
        if isClamped(inbound_interface)
            # Hard-code outbound message of constant node in schedule
            push!(inbounds, assembleClamp!(inbound_interface.node, Message))
        else
            # Collect message from previous result
            push!(inbounds, interface_to_schedule_entry[inbound_interface])
        end
    end

    return inbounds
end


text:
## work with Vector{Uint8} via I/O primitives ##

# Stateful string
type IOBuffer <: IO
    data::Vector{Uint8}
    readable::Bool
    writable::Bool
    seekable::Bool # if not seekable, free to destroy (compact) past read data
    append::Bool # add data at end instead of at pointer
    size::Int
    maxsize::Int # pre-allocated, fixed array size
    ptr::Int # read (and maybe write) pointer
    IOBuffer(data::Vector{Uint8},readable::Bool,writable::Bool,seekable::Bool,append::Bool,maxsize::Int) =
        new(data,readable,writable,seekable,append,length(data),maxsize,1)
end

function copy(b::IOBuffer)
    ret = IOBuffer(b.writable?copy(b.data):b.data,
            b.readable,b.writable,b.seekable,b.append,b.maxsize)
    ret.size = b.size
    ret.ptr  = b.ptr
    ret
end

# PipeBuffers behave like Unix Pipes. They are readable and writable, the act appendable, and not seekable.
PipeBuffer(data::Vector{Uint8},maxsize::Int) = IOBuffer(data,true,true,false,true,maxsize)
PipeBuffer(data::Vector{Uint8}) = PipeBuffer(data,typemax(Int))
PipeBuffer() = PipeBuffer(Uint8[])
PipeBuffer(maxsize::Int) = (x = PipeBuffer(Array(Uint8,maxsize),maxsize); x.size=0; x)

# IOBuffers behave like Files. They are readable and writable. They are seekable. (They can be appendable).
IOBuffer(data::Vector{Uint8},readable::Bool,writable::Bool,maxsize::Int) =
        IOBuffer(data,readable,writable,true,false,maxsize)
IOBuffer(data::Vector{Uint8},readable::Bool,writable::Bool) = IOBuffer(data,readable,writable,typemax(Int))
IOBuffer(data::Vector{Uint8}) = IOBuffer(data, true, false)
IOBuffer(str::ByteString) = IOBuffer(str.data, true, false)
IOBuffer(readable::Bool,writable::Bool) = IOBuffer(Uint8[],readable,writable)
IOBuffer() = IOBuffer(Uint8[], true, true)
IOBuffer(maxsize::Int) = (x=IOBuffer(Array(Uint8,maxsize),true,true,maxsize); x.size=0; x)

read!(from::IOBuffer, a::Array) = read_sub(from, a, 1, length(a))

function read_sub{T}(from::IOBuffer, a::Array{T}, offs, nel)
    if offs+nel-1 > length(a) || offs < 1 || nel < 0
        throw(BoundsError())
    end
    if !isbits(T)
        error("read from IOBuffer only supports bits types or arrays of bits types; got "*string(T))
    end
    read!(from, pointer(a, offs), nel*sizeof(T))
    return a
end

read!(from::IOBuffer, p::Ptr, nb::Integer) = read!(from, p, int(nb))
function read!(from::IOBuffer, p::Ptr, nb::Int)
    if !from.readable error("read failed") end
    avail = nb_available(from)
    adv = min(avail,nb)
    ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, Uint), p, pointer(from.data,from.ptr), adv)
    from.ptr += adv
    if nb > avail
        throw(EOFError())
    end
    p
end

function read(from::IOBuffer, ::Type{Uint8})
    if !from.readable error("read failed") end
    if from.ptr > from.size
        throw(EOFError())
    end
    byte = from.data[from.ptr]
    from.ptr += 1
    return byte
end

function peek(from::IOBuffer)
    if !from.readable error("read failed") end
    if from.ptr > from.size
        throw(EOFError())
    end
    return from.data[from.ptr]
end

read{T}(from::IOBuffer, ::Type{Ptr{T}}) = convert(Ptr{T}, read(from, Uint))

isreadable(io::IOBuffer) = io.readable
iswritable(io::IOBuffer) = io.writable

# TODO: IOBuffer is not iterable, so doesn't really have a length.
# This should maybe be sizeof() instead.
#length(io::IOBuffer) = (io.seekable ? io.size : nb_available(io))
nb_available(io::IOBuffer) = io.size - io.ptr + 1
position(io::IOBuffer) = io.ptr-1

function skip(io::IOBuffer, n::Integer)
    io.ptr = min(io.ptr + n, io.size+1)
    return io
end
function seek(io::IOBuffer, n::Integer)
    if !io.seekable error("seek failed") end
    io.ptr = min(n+1, io.size+1)
    return io
end
function seekend(io::IOBuffer)
    io.ptr = io.size+1
    return io
end
function truncate(io::IOBuffer, n::Integer)
    if !io.writable error("truncate failed") end
    if !io.seekable error("truncate failed") end #because absolute offsets are meaningless
    if n > io.maxsize || n < 0 error("truncate failed") end
    if n > length(io.data)
        resize!(io.data, n)
    end
    io.data[io.size+1:n] = 0
    io.size = n
    io.ptr = min(io.ptr, n+1)
    return io
end
function compact(io::IOBuffer)
    if !io.writable error("compact failed") end
    if io.seekable error("compact failed") end
    ccall(:memmove, Ptr{Void}, (Ptr{Void},Ptr{Void},Uint),
          io.data, pointer(io.data,io.ptr), nb_available(io))
    io.size -= io.ptr - 1
    io.ptr = 1
    return io
end
function ensureroom(io::IOBuffer, nshort::Int)
    if !io.writable error("ensureroom failed") end
    if !io.seekable
        if nshort < 0 error("ensureroom failed") end
        if io.ptr > 1 && io.size <= io.ptr - 1
            io.ptr = 1
            io.size = 0
        elseif (io.size+nshort > io.maxsize) ||
                (io.ptr > io.size - io.ptr > 4096) ||
                (io.ptr > 262144)
            # apply somewhat arbitrary heuristics to decide when to destroy
            # old, read data to make more room for new data
            compact(io)
        end
    end
    n = min(nshort + (io.append ? io.size : io.ptr-1), io.maxsize)
    if n > length(io.data)
        resize!(io.data, n)
    end
    return io
end
eof(io::IOBuffer) = (io.ptr-1 == io.size)
function close(io::IOBuffer)
    if io.writable
        resize!(io.data, 0)
    else
        io.data = Uint8[]
    end
    io.readable = false
    io.writable = false
    io.seekable = false
    io.size = 0
    io.maxsize = 0
    io.ptr = 1
    nothing
end
isopen(io::IOBuffer) = io.readable || io.writable || io.seekable || nb_available(io) > 0
function bytestring(io::IOBuffer)
    if !io.readable error("bytestring read failed") end
    if !io.seekable error("bytestring read failed") end
    bytestring(pointer(io.data), io.size)
end
function takebuf_array(io::IOBuffer)
    if io.seekable
        data = io.data
        if io.writable
            maxsize = (io.maxsize == typemax(Int) ? 0 : io.maxsize)
            io.data = Array(Uint8,maxsize)
        else
            data = copy(data)
        end
        resize!(data,io.size)
    else
        nbytes = nb_available(io)
        a = Array(Uint8, nbytes)
        data = read!(io, a)
    end
    if io.writable
        io.ptr = 1
        io.size = 0
    end
    data
end
takebuf_string(io::IOBuffer) = bytestring(takebuf_array(io))

write(to::IOBuffer, p::Ptr, nb::Integer) = write(to, p, int(nb))
function write(to::IOBuffer, p::Ptr, nb::Int)
    !to.writable && error("write failed")
    ensureroom(to, nb)
    ptr = (to.append ? to.size+1 : to.ptr)
    nb = min(nb, length(to.data) - ptr + 1)
    ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, Uint), pointer(to.data,ptr), p, nb)
    to.size = max(to.size, ptr - 1 + nb)
    if !to.append; to.ptr += nb; end
    nb
end

function write_sub{T}(to::IOBuffer, a::Array{T}, offs, nel)
    if offs+nel-1 > length(a) || offs < 1 || nel < 0
        throw(BoundsError())
    end
    if !isbits(T)
        error("write to IOBuffer only supports bits types or arrays of bits types; got "*string(T))
    end
    write(to, pointer(a,offs), nel*sizeof(T))
end

write(to::IOBuffer, a::Array) = write_sub(to, a, 1, length(a))

function write(to::IOBuffer, a::Uint8)
    if !to.writable error("write failed") end
    ensureroom(to, 1)
    ptr = (to.append ? to.size+1 : to.ptr)
    if ptr > to.maxsize
        return 0
    else
        to.data[ptr] = a
    end
    to.size = max(to.size, ptr)
    if !to.append; to.ptr += 1; end
    sizeof(Uint8)
end

write(to::IOBuffer, p::Ptr) = write(to, convert(Uint, p))

function readbytes!(io::IOBuffer, b::Array{Uint8}, nb=length(b))
    nr = min(nb, nb_available(io))
    if length(b) < nr
        resize!(b, nr)
    end
    read_sub(io, b, 1, nr)
    return nr
end
readbytes(io::IOBuffer) = read!(io, Array(Uint8, nb_available(io)))
readbytes(io::IOBuffer, nb) = read!(io, Array(Uint8, min(nb, nb_available(io))))

function search(buf::IOBuffer, delim)
    p = pointer(buf.data, buf.ptr)
    q = ccall(:memchr,Ptr{Uint8},(Ptr{Uint8},Int32,Csize_t),p,delim,nb_available(buf))
    nb = (q == C_NULL ? 0 : q-p+1)
end
function readuntil(io::IOBuffer, delim::Uint8)
    nb = search(io, delim)
    if nb == 0
        nb = nb_available(io)
    end
    read!(io, Array(Uint8, nb))
end

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia


text:
export assembly_dir,
    download_assembly,
    load_assembly

include("config.jl")
include("downloader.jl")
include("reader.jl")