text:
# Input binary video
# Output different ids for different sweat glands production
# Start from first frame,
# grow spots frame by frame,
# merge spots if they are too small (sometimes sweat jumps over some pixels)

@inline function kernelvote4(q, mask, output, x, y)
    q0 = @inbounds q[x, y];
    if (@inbounds mask[x, y] == 0) return 0 end

    if (q0 != 0)
        @inbounds output[x, y] = q0;
        return 0;
    end
    kernelvote4inner(q, mask, output, x, y)
end
function sort4(i1, i2, i3, i4)
    # Probably bitonic sort
    j1 = min(i1, i2);
    j2 = max(i1, i2);
    j3 = min(i3, i4);
    j4 = max(i3, i4);
    i2     = max(j1, j3);
    j1     = min(j1, j3);
    i3     = min(j2, j4);
    j4     = max(j2, j4);
    j2 = min(i2, i3);
    j3 = max(i2, i3);
    j1, j2, j3, j4
end
function kernelvote4inner(q, mask, output, x, y)
    w, h = size(q)
    i1 = x < w && mask[x + 1, y + 0] != 0 ? q[x + 1, y + 0] : 0
    i2 = x > 1 && mask[x - 1, y + 0] != 0 ? q[x - 1, y + 0] : 0
    i3 = y < h && mask[x + 0, y + 1] != 0 ? q[x + 0, y + 1] : 0
    i4 = y > 1 && mask[x + 0, y - 1] != 0 ? q[x + 0, y - 1] : 0

    j1, j2, j3, j4 = sort4(i1, i2, i3, i4)
    o = 0;
    if (j1 != 0 && j1 == j2) o = j2;
    elseif (j3 != 0) o = j3;
    elseif (j4 != 0) o = j4;
    end
    if (o != 0)
        output[x, y] = o;
        return 1
    end
    0
end

function _iter(A, mask, tmp)
    # amask = convert(Array, mask)
    # aA = convert(Array, A)
    b1 = A
    b2 = copy!(tmp, A)#copy(A)
    i = 0
    while i < 100
        changed = 0
        for y in axes(A, 2)
            for x in axes(A, 1)
                changed += kernelvote4(b1, mask, b2, x, y)
            end
        end

        if changed == 0
            break
        end
        i += 1
        b1, b2 = b2, b1
    end
    return i, b1
end

function segment_blobs(mask)
    f = Int32.(mask)
    tmp = Matrix{Int32}(undef, size(f, 1), size(f, 2))
    tmp2 = similar(tmp)
    _, clusterid = doclustering!(view(f, :,:,1))
    for frame in 2:size(f, 3)
        #println(frame)
        _, clusterid = doclustering!(view(f, :,:,frame)) #, clusterid)
        #return f[:,:,frame], f[:,:,frame-1]
        matchclusterids!(view(f, :,:,frame), view(f, :,:,frame-1), tmp, tmp2)
    end
    println()
    f
end

function matchclusterids_part1(f2, f1)
    matchclusterids_part1!(similar(f1), f2, f1)
end
function matchclusterids_part1!(f1c, f2, f1)
    #mask = @. Int32(f2 != 0)
    maxid = max(maximum(f2), maximum(f1))
    toosmall, new = computetoosmallandnew(f2, f1, maxid)
    # Remove small clusters from f1
    # Copy new clusters

    newids = Dict{Int, Int}()
    # f1c = similar(f1)
    for i in eachindex(f1c)
        a = f1[i]
        b = f2[i]
        if b != 0
            if new[b]
                c = get!(newids, b, maxid+1)
                if c == maxid+1
                    maxid += one(maxid)
                end
                f1c[i] = c
                # f1c[i] = get!(newids, b) do
                #     maxid += 1
                #     maxid
                # end
            elseif a == 0 || !toosmall[a]
                f1c[i] = a
            else
                f1c[i] = zero(eltype(f1))
            end
        else
            @assert a == 0
        end
    end
    f1c
end

function matchclusterids!(f2, f1, tmp, tmp2)
    f1c = matchclusterids_part1!(tmp, f2, f1)

    #println(sizesnew, sizesold)
    mask = f2
    i, B = _iter(f1c, mask, tmp2)
    print(i, ",") #, ":", count(!iszero, mask)
    copy!(f2, B)
end

# Count the number of pixels in each blob of A
#  # (And the number of pixels that are also in B)
# Id of the largest (too small) blob in B or 0 otherwise.
#  Too small = all but largest
function computetoosmallandnew(f2::AbstractArray{T}, f1::AbstractArray{T}, maxid::T) where T
    @assert maxid < 1_000_000
    toosmall = falses(maxid)
    new = falses(maxid)

    # Create a Dict with new cluster ids -> Dict(old cluster ids -> count)
    dict = [Dict{T, Int}() for i in 1:maxid]
    for (a,b) in zip(f1, f2)
        if a != 0 #&& b != 0
            d = dict[b]
            count = get!(d, a, 0)
            if count <= minimum_blob_size + 2# No need to count more than used in the next for loop
                d[a] = count+1
            end
        end
    end
    # The Vector{Dict} dict could be simplified to two Vector{Int}:s maybe:
    # testi = reduce(vcat, collect.(keys.(dict)))
    # @assert allunique(testi)
    for (i, d) in enumerate(dict)
        if length(d) == 0
            new[i] = true
        else
            count, key = maximum((count, k) for (k, count) in d)
            for (k, v) in d
                if v < minimum_blob_size && (v, k) < (count, key)
                    toosmall[k] = true
                end
            end
        end
    end
    toosmall, new
end


text:
"""
    log_value(logger, name, value; step=step(logger))

Logs a Floating-point variable with name `name` at step `step`.
If `value` is a complex number, the real and imaginary part are logged
separately.
"""
function log_value(logger::TBLogger, name::AbstractString, value::Real; step=nothing)
    summ = SummaryCollection(scalar_summary(name, value))
    write_event(logger.file, make_event(logger, summ, step=step))
end

function log_value(logger::TBLogger, name::AbstractString, value::Complex; step=nothing)
    log_value(logger, name*"/re", real(value), step=step)
    log_value(logger, name*"/im", imag(value), step=step)
end

function scalar_summary(name::AbstractString, value::Real)
    Summary_Value(tag=name, simple_value=value)
end

ext: jl


text:
# Krylov phiv with internal time-stepping
"""
    exp_timestep(ts,A,b[;adaptive,tol,kwargs...]) -> U

Evaluates the matrix exponentiation-vector product using time stepping

```math
u = \\exp(tA)b
```

`ts` is an array of time snapshots for u, with `U[:,j] ≈ u(ts[j])`. `ts` can
also be just one value, in which case only the end result is returned and `U`
is a vector.

The time stepping formula of Niesen & Wright is used [^1]. If the time step
`tau` is not specified, it is chosen according to (17) of Neisen & Wright. If
`adaptive==true`, the time step and Krylov subsapce size adaptation scheme of
Niesen & Wright is used, the relative tolerance of which can be set using the
keyword parameter `tol`. The delta and gamma parameter of the adaptation
scheme can also be adjusted.

Set `verbose=true` to print out the internal steps (for debugging). For the
other keyword arguments, consult `arnoldi` and `phiv`, which are used
internally.

Note that this function is just a special case of `phiv_timestep` with a more
intuitive interface (vector `b` instead of a n-by-1 matrix `B`).

[^1]: Niesen, J., & Wright, W. (2009). A Krylov subspace algorithm for
evaluating the φ-functions in exponential integrators. arXiv preprint
arXiv:0907.4631.
"""
function expv_timestep(ts::Vector{tType}, A, b; kwargs...) where {tType <: Real}
    U = Matrix{eltype(b)}(undef, size(A, 1), length(ts))
    expv_timestep!(U, ts, A, b; kwargs...)
end
function expv_timestep(t::tType, A, b; kwargs...) where {tType <: Real}
    u = Vector{eltype(b)}(undef, size(A, 1))
    expv_timestep!(u, t, A, b; kwargs...)
end
"""
    expv_timestep!(u,t,A,b[;kwargs]) -> u

Non-allocating version of `expv_timestep`.
"""
function expv_timestep!(u::AbstractVector{T}, t::tType, A, b::AbstractVector{T};
                        kwargs...) where {T <: Number, tType <: Real}
    expv_timestep!(reshape(u, length(u), 1), [t], A, b; kwargs...)
    return u
end
function expv_timestep!(U::AbstractMatrix{T}, ts::Vector{tType}, A, b::AbstractVector{T};
                        kwargs...) where {T <: Number, tType <: Real}
    B = reshape(b, length(b), 1)
    phiv_timestep!(U, ts, A, B; kwargs...)
end
"""
    phiv_timestep(ts,A,B[;adaptive,tol,kwargs...]) -> U

Evaluates the linear combination of phi-vector products using time stepping

```math
u = \\varphi_0(tA)b_0 + t\\varphi_1(tA)b_1 + \\cdots + t^p\\varphi_p(tA)b_p
```

`ts` is an array of time snapshots for u, with `U[:,j] ≈ u(ts[j])`. `ts` can
also be just one value, in which case only the end result is returned and `U`
is a vector.

The time stepping formula of Niesen & Wright is used [^1]. If the time step
`tau` is not specified, it is chosen according to (17) of Neisen & Wright. If
`adaptive==true`, the time step and Krylov subsapce size adaptation scheme of
Niesen & Wright is used, the relative tolerance of which can be set using the
keyword parameter `tol`. The delta and gamma parameter of the adaptation
scheme can also be adjusted.

Set `verbose=true` to print out the internal steps (for debugging). For the
other keyword arguments, consult `arnoldi` and `phiv`, which are used
internally.

[^1]: Niesen, J., & Wright, W. (2009). A Krylov subspace algorithm for
evaluating the φ-functions in exponential integrators. arXiv preprint
arXiv:0907.4631.
"""
function phiv_timestep(ts::Vector{tType}, A, B; kwargs...) where {tType <: Real}
    U = Matrix{eltype(B)}(undef, size(A, 1), length(ts))
    phiv_timestep!(U, ts, A, B; kwargs...)
end
function phiv_timestep(t::tType, A, B; kwargs...) where {tType <: Real}
    u = Vector{eltype(B)}(undef, size(A, 1))
    phiv_timestep!(u, t, A, B; kwargs...)
end
"""
    phiv_timestep!(U,ts,A,B[;kwargs]) -> U

Non-allocating version of `phiv_timestep`.
"""
function phiv_timestep!(u::AbstractVector{T}, t::tType, A, B::AbstractMatrix{T};
                        kwargs...) where {T <: Number, tType <: Real}
    phiv_timestep!(reshape(u, length(u), 1), [t], A, B; kwargs...)
    return u
end
function phiv_timestep!(U::AbstractMatrix{T}, ts::Vector{tType}, A, B::AbstractMatrix{T}; tau::Real=0.0,
                        m::Int=min(10, size(A, 1)), tol::Real=1e-7, opnorm=LinearAlgebra.opnorm, iop::Int=0,
                        correct::Bool=false, caches=nothing, adaptive=false, delta::Real=1.2,
                        gamma::Real=0.8, NA::Int=0, verbose=false) where {T <: Number, tType <: Real}
    # Choose initial timestep
    abstol = tol * opnorm(A, Inf)
    verbose && println("Absolute tolerance: $abstol")
    if iszero(tau)
        Anorm = opnorm(A, Inf)
        b0norm = norm(@view(B[:, 1]), Inf)
        tau = 10/Anorm * (abstol * ((m+1)/ℯ)^(m+1) * sqrt(2*pi*(m+1)) /
                          (4*Anorm*b0norm))^(1/m)
        verbose && println("Initial time step unspecified, chosen to be $tau")
    end
    # Initialization
    n = size(U, 1)
    sort!(ts); tend = ts[end]
    p = size(B, 2) - 1
    @assert length(ts) == size(U, 2) "Dimension mismatch"
    @assert n == size(A, 1) == size(A, 2) == size(B, 1) "Dimension mismatch"
    if caches == nothing
        u = Vector{T}(undef, n)              # stores the current state
        W = Matrix{T}(undef, n, p+1)         # stores the w vectors
        P = Matrix{T}(undef, n, p+2)         # stores output from phiv!
        Ks = KrylovSubspace{T}(n, m)  # stores output from arnoldi!
        phiv_cache = nothing         # cache used by phiv!
    else
        u, W, P, Ks, phiv_cache = caches
        @assert length(u) == n && size(W, 1) == n && size(P, 1) == n "Dimension mismatch"
        # W and P may be bigger than actually needed
        W = @view(W[:, 1:p+1])
        P = @view(P[:, 1:p+2])
    end
    copyto!(u, @view(B[:, 1])) # u(0) = b0
    coeffs = ones(tType, p);
    if adaptive # initialization step for the adaptive scheme
        if ishermitian(A)
            iop = 2 # does not have an effect on arnoldi!, just for flops estimation
        end
        if iszero(NA)
            _A = convert(AbstractMatrix, A)
            if isa(_A, SparseMatrixCSC)
                NA = nnz(_A)
            else
                NA = count(!iszero, _A) # not constant operation, should be best avoided
            end
        end
    end

    t = 0.0       # current time
    snapshot = 1  # which snapshot to compute next
    while t < tend # time stepping loop
        if t + tau > tend # last step
            tau = tend - t
        end
        # Part 1: compute w0...wp using the recurrence relation (16)
        copyto!(@view(W[:, 1]), u) # w0 = u(t)
        @inbounds for l = 1:p-1 # compute cl = t^l/l!
            coeffs[l+1] = coeffs[l] * t / l
        end
        @views @inbounds for j = 1:p
            mul!(W[:, j+1], A, W[:, j])
            for l = 0:p-j
                axpy!(coeffs[l+1], B[:, j+l+1], W[:, j+1])
            end
        end
        # Part 2: compute ϕp(tau*A)wp using Krylov, possibly with adaptation
        arnoldi!(Ks, A, @view(W[:, end]); tol=tol, m=m, opnorm=opnorm, iop=iop)
        _, epsilon = phiv!(P, tau, Ks, p + 1; cache=phiv_cache, correct=correct, errest=true)
        verbose && println("t = $t, m = $m, tau = $tau, error estimate = $epsilon")
        if adaptive
            omega = (tend / tau) * (epsilon / abstol)
            epsilon_old = epsilon; m_old = m; tau_old = tau
            q = m/4; kappa = 2.0; maxtau = tend - t
            while omega > delta # inner loop of Algorithm 3
                m_new, tau_new, q, kappa = _phiv_timestep_adapt(
                    m, tau, epsilon, m_old, tau_old, epsilon_old, q, kappa,
                    gamma, omega, maxtau, n, p, NA, iop, opnorm(getH(Ks), 1), verbose)
                m, m_old = m_new, m
                tau, tau_old = tau_new, tau
                # Compute ϕp(tau*A)wp using the new parameters
                arnoldi!(Ks, A, @view(W[:, end]); tol=tol, m=m, opnorm=opnorm, iop=iop)
                _, epsilon_new = phiv!(P, tau, Ks, p + 1; cache=phiv_cache, correct=correct, errest=true)
                epsilon, epsilon_old = epsilon_new, epsilon
                omega = (tend / tau) * (epsilon / abstol)
                verbose && println("  * m = $m, tau = $tau, error estimate = $epsilon")
            end
        end
        # Part 3: update u using (15)
        lmul!(tau^p, copyto!(u, @view(P[:, end - 1])))
        @inbounds for l = 1:p-1 # compute cl = tau^l/l!
            coeffs[l+1] = coeffs[l] * tau / l
        end
        @views @inbounds for j = 0:p-1
            axpy!(coeffs[j+1], W[:, j+1], u)
        end
        # Fill out all snapshots in between the current step
        while snapshot <= length(ts) && t + tau >= ts[snapshot]
            tau_snapshot = ts[snapshot] - t
            u_snapshot = @view(U[:, snapshot])
            phiv!(P, tau_snapshot, Ks, p + 1; cache=phiv_cache, correct=correct)
            lmul!(tau_snapshot^p, copyto!(u_snapshot, @view(P[:, end - 1])))
            @inbounds for l = 1:p-1 # compute cl = tau^l/l!
                coeffs[l+1] = coeffs[l] * tau_snapshot / l
            end
            @views @inbounds for j = 0:p-1
                axpy!(coeffs[j+1], W[:, j+1], u_snapshot)
            end
            snapshot += 1
        end

        t += tau
    end

return U
end
# Helper functions for phiv_timestep!
function _phiv_timestep_adapt(m, tau, epsilon, m_old, tau_old, epsilon_old, q, kappa,
                              gamma, omega, maxtau, n, p, NA, iop, Hnorm, verbose)
    # Compute new m and tau (Algorithm 4)
    if tau_old > tau
        q = log(tau/tau_old) / log(epsilon/epsilon_old) - 1
    end # else keep q the same
    tau_new = tau * (gamma / omega)^(1/(q + 1))
    tau_new = min(max(tau_new, tau/5), 2*tau, maxtau)
    if m_old < m
        kappa = (epsilon/epsilon_old)^(1/(m_old - m))
    end # else keep kappa the same
    m_new = m + ceil(Int, log(omega / gamma) / log(kappa))
    m_new = min(max(m_new, div(3*m, 4), 1), Int(ceil(4*m / 3)))
    verbose && println("  - Proposed new m: $m_new, new tau: $tau_new")
    # Compare costs of using new m vs new tau (23)
    cost_tau = _phiv_timestep_estimate_flops(m, tau_new, n, p, NA, iop, Hnorm, maxtau)
    cost_m = _phiv_timestep_estimate_flops(m_new, tau, n, p, NA, iop, Hnorm, maxtau)
    verbose && println("  - Cost to use new m: $cost_m flops, new tau: $cost_tau flops")
    if cost_tau < cost_m
        m_new = m
    else
        tau_new = tau
    end
    return m_new, tau_new, q, kappa
end
function _phiv_timestep_estimate_flops(m, tau, n, p, NA, iop, Hnorm, maxtau)
    # Estimate flops for the update of W and u
    flops_W = 2 * (p - 1) * (NA + n)
    flops_u = (2 * p + 1) * n
    # Estimate flops for arnoldi!
    if iop == 0
        iop = m
    end
    flops_matvec = 2 * m * NA
    flops_vecvec = 0
    for i = 1:m
        flops_vecvec += 3 * min(i, iop)
    end
    # Estimate flops for phiv! (7)
    MH = 44/3 + 2 * ceil(max(0.0, log2(Hnorm / 5.37)))
    flops_phiv = round(Int, MH * (m + p)^3)

    flops_onestep = flops_W + flops_u + flops_matvec + flops_vecvec + flops_phiv
    return flops_onestep * Int(ceil(maxtau / tau))
end
function _phiv_timestep_caches(u_prototype, maxiter::Int, p::Int)
    n = length(u_prototype); T = eltype(u_prototype)
    u = similar(u_prototype)                      # stores the current state
    W = Matrix{T}(undef, n, p+1)                  # stores the w vectors
    P = Matrix{T}(undef, n, p+2)                  # stores output from phiv!
    Ks = KrylovSubspace{T}(n, maxiter)            # stores output from arnoldi!
    phiv_cache = PhivCache{T}(maxiter, p+1)       # cache used by phiv! (need +1 for error estimation)
    return u, W, P, Ks, phiv_cache
end


text:
using PyCall

function gymspace2jlspace(s::PyObject)
    classname = s[:__class__][:__name__]
    if classname == "Box"
        BoxSpace(s[:low], s[:high])
    elseif classname == "Discrete"
        DiscreteSpace(s[:n])
    elseif classname == "MultiBinary"
        MultiBinarySpace(s[:n])
    elseif classname == "MultiDiscrete"
        MultiDiscreteSpace(s[:nvec])
    elseif classname == "Tuple"
        map(gymspace2jlspace, s[:spaces])
    elseif classname == "Dict"
        Dict(map((k, v) -> (k, gymspace2jlspace(v)), s[:spaces]))
    end
end


text:
import Pkg; Pkg.add("Documenter")
push!(LOAD_PATH,"../src/")
using Documenter, VTUFileHandler
makedocs(
	sitename = "VTUFileHandler.jl",
	modules = [VTUFileHandler],
	pages = [
		"Home" => "index.md"
		"Library" => "lib/lib.md"
	]
	)
deploydocs(
    repo = "github.com/baxmittens/VTUFileHandler.jl.git"
)


text:
function backwards_tfunc(@nospecialize(f), primal::IRCode, inst::Expr, @nospecialize(Δ))
    if f === Core.tuple
        return tuple_tfunc(Any[NoTangent, tuple_type_fields(Δ)...])
    elseif f == Core.getfield
        tt = argextype(inst.args[2], primal, primal.sptypes)
        ot = widenconst(tt)
        if ot <: Tuple
            if length(ot.parameters) == 1
                # Only one choice
                rt = tuple_tfunc(Any[Δ])
            else
                idxt = argextype(inst.args[3], primal, primal.sptypes)
                isa(idxt, Const) || error()
                args = Any[ZeroTangent for i = 1:length(ot.parameters)]
                args[idxt.val] = Δ
                rt = tuple_tfunc(args)
            end
        else
            idxt = argextype(inst.args[3], primal, primal.sptypes)
            if !isa(idxt, Const)
                @show idxt
                error()
            end
            if isa(idxt.val, Symbol)
                symn = idxt.val
            elseif isa(idxt.val, Int)
                @show tt
                symn = fieldname(tt, idxt.val)
            else
                error()
            end
            rt = NamedTuple{(symn,), Tuple{widenconst(Δ)}}
        end
        return tuple_tfunc(Any[NoTangent, rt, NoTangent, NoTangent])
    elseif f == Core.apply_type
        return tuple_tfunc(Any[NoTangent for i = 1:length(inst.args)])
    elseif f == Core.typeof
        return tuple_tfunc(Any[ZeroTangent for i = 1:length(inst.args)])
    elseif f == (===)
        return tuple_tfunc(Any[ZeroTangent for i = 1:length(inst.args)])
    elseif f == nfields
        return tuple_tfunc(Any[ZeroTangent for i = 1:length(inst.args)])
    elseif f == fieldtype
        return tuple_tfunc(Any[ZeroTangent for i = 1:length(inst.args)])
    elseif f === Core._apply_iterate
        ft = argextype(inst.args[3], primal, primal.sptypes)
        f = singleton_type(ft)
        @assert f === Core.tuple
        return backwards_tfunc(f, primal, inst, Δ)
    end
    @show f
    error()
end

function forward_tfunc(@nospecialize(f), primal::IRCode, inst::Expr, @nospecialize(Δ))
    if f === Core.tuple
        return tuple_tfunc(Any[tuple_type_fields(Δ)[2:end]...])
    elseif f == Core.getfield
        idxt = argextype(inst.args[3], primal, primal.sptypes)
        isa(idxt, Const) || error()
        val = getfield_tfunc(Δ, Const(2))
        fieldval = (isa(val, Const) && (isa(val.val, ZeroTangent) || isa(val.val, NoTangent))) ?
            Const(ZeroTangent()) : getfield_tfunc(val, idxt)
        return fieldval
    elseif f === Core.apply_type || f === Core.typeof
        return NoTangent
    end
    @show f
    error()
end

function getfield_prop_zero_tfunc(@nospecialize(s), @nospecialize(x))
    if isa(s, Const)
        if isa(s.val, NoTangent) || isa(s.val, ZeroTangent)
            return s
        end
    end
    if s === ZeroTangent || s === NoTangent
        return s
    end
    return getfield_tfunc(s, x)
end


text:
# TODO: test display error handling stuff, here are only basic execution tests

evaltestpath = joinpath′(@__DIR__, "fixtures", "EvalTest.jl")

@testset "eval" begin

@testset "evalall" begin
    code = read(evaltestpath, String)
    Atom.evalall(code, "Main", evaltestpath)
    @test isdefined(Main, :EvalTest)
    @test isdefined(EvalTest, :foo)
end

@testset "eval" begin
    let ret = Atom.eval("bar = :bar", 4, evaltestpath, "Main.EvalTest")
        @test ret[:text] == ":bar"
    end
    @test isdefined(EvalTest, :bar)
    let ret = Atom.eval("bar = :bar2", 4, evaltestpath, "Main.EvalTest")
        @test ret[:text] == ":bar2"
    end
    @test EvalTest.bar == :bar2
end

@testset "evalshow" begin
    # NOTE: semicolons are necessary so that they work without repl display during a test
    Atom.evalshow("baz = :baz;", 5, evaltestpath, "Main.EvalTest")
    @test isdefined(EvalTest, :baz)
    Atom.evalshow("baz = :baz2;", 5, evaltestpath, "Main.EvalTest")
    @test EvalTest.baz == :baz2
end

end

@testset "toggle docs" begin
    # basic
    let doc = Atom.docs("push!")
        @test !doc[:error]
    end

    # context module
    let doc = Atom.docs("handle", "Atom")
        @test !doc[:error]
    end

    # keyword
    let doc = Atom.docs("begin", "Atom")
        @test !doc[:error]
    end
end


text:
"""
	struct Monomer

A chemical monomer type, identified by a unique monomer ID
that accesses a field and density grid from the system.
"""
struct Monomer
	id     :: Int
	vol    :: Float64
	charge :: Float64

	function Monomer(; id::Integer, vol::Real = 1.0, charge::Real = 0.0)
		return new(id, vol, charge)
	end
end

#==============================================================================#

Base.show(io::IO, mon::Monomer) = @printf(io, "Monomer(mid = %d, vol = %.2f, charge = %.2f)", mon.id, mon.vol, mon.charge)

Base.hash(mon::Monomer, h::UInt) = hash(mon.charge, hash(mon.vol, hash(mon.id, h)))

Base.:(==)(a::Monomer, b::Monomer) = (hash(a) == hash(b))

Base.isless(a::Monomer, b::Monomer) = (a.id < b.id)


text:
# Rooted tree
#   parent of root is self
#   the kids of node i are children[kidsPtr[i] : kidsPtr[i]-numKids[i]-1 ]
#  the weigts of edges from parents to kids are indexed similarly.
#  that is, weights holds the weight of an edge to the parent.
#
#  We require that the children appears in a dfs order
#

type RootedTree{Tval,Tind}
  root::Tind
  parent::Array{Tind,1}
  children::Array{Tind,1}
  weights::Array{Tval,1}
  numKids::Array{Tind,1}
  kidsPtr::Array{Tind,1}
end # RootedTree

# the following two routines will return the children of a node,
# and the weights on the edges to them.
# but, using these can be slow as it seems to cause memory allocation

kids(tr::RootedTree, i) = tr.children[tr.kidsPtr[i]:(tr.kidsPtr[i]+tr.numKids[i] - 1)]
weights(tr::RootedTree, i) = tr.weights[tr.kidsPtr[i]:(tr.kidsPtr[i]+tr.numKids[i] - 1)]

# matToTree takes a weighted tree in mat format, and produces
# a RootedTree from it

matToTree{Tv,Ti}(mat::SparseMatrixCSC{Tv,Ti}) = matToTree(mat, 1::Ti)

function matToTree{Tv,Ti}(mat::SparseMatrixCSC{Tv,Ti}, root::Ti)

  n = mat.n
  visited = zeros(Bool,n)

  children = zeros(Ti,n)
  weights = zeros(Tv,n)
  parent = zeros(Ti,n)
  numKids = zeros(Ti,n)
  kidsPtr = zeros(Ti,n)

  children[1] = root
  visited[root] = true
  parent[root] = root
  kidsPtr[root] = 2

  degRoot = deg(mat,root)
  numKids[root] = degRoot

  children[1 + collect(1:degRoot)] = mat.rowval[mat.colptr[root]:(mat.colptr[root+1]-1)]
  weights[1 + collect(1:degRoot)] = mat.nzval[mat.colptr[root]:(mat.colptr[root+1]-1)]

  numIn = 1+degRoot

  parent[mat.rowval[mat.colptr[root]:(mat.colptr[root+1]-1)]] = root
  visited[mat.rowval[mat.colptr[root]:(mat.colptr[root+1]-1)]] = true

  ptr = 2

  while ptr <= numIn
    v = children[ptr]
    numKids[v] = deg(mat,v)-1

    if numKids[v] > 0
        kidsPtr[v] = numIn+1
    end

    for ind in mat.colptr[v]:(mat.colptr[v+1]-1)
      nbr = mat.rowval[ind]
      if !visited[nbr]
          numIn = numIn+1
          visited[nbr] = true
          children[numIn] = nbr
          weights[numIn] = mat.nzval[ind]

          parent[nbr] = v
      end # if
    end # for

    ptr = ptr+1

  end # while

  if numIn < n
      error("graph is not connected")
  end

  return RootedTree(root, parent, children, weights, numKids, kidsPtr)

end

# create a rooted tree structure, and compute the depth as we go
# this is faster than doing it independently

matToTreeDepth{Tv,Ti}(mat::SparseMatrixCSC{Tv,Ti}) = matToTreeDepth(mat, 1::Ti)

function matToTreeDepth{Tv,Ti}(mat::SparseMatrixCSC{Tv,Ti}, root::Ti)

  n = mat.n
  visited = zeros(Bool,n)

  children = zeros(Ti,n)
  weights = zeros(Tv,n)
  parent = zeros(Ti,n)
  numKids = zeros(Ti,n)
  kidsPtr = zeros(Ti,n)

  depth = zeros(Tv,n)

  children[1] = root
  visited[root] = true
  parent[root] = root
  kidsPtr[root] = 2

  degRoot = deg(mat,root)
  numKids[root] = degRoot

  kids = mat.rowval[mat.colptr[root]:(mat.colptr[root+1]-1)]
  wts =  mat.nzval[mat.colptr[root]:(mat.colptr[root+1]-1)]

  children[1 + collect(1:degRoot)] = kids
  weights[1 + collect(1:degRoot)] = wts

  numIn = 1+degRoot

  parent[kids] = root
  visited[kids] = true
  depth[kids] = 1./wts

  ptr = 2

  while ptr <= numIn
    v = children[ptr]
    numKids[v] = deg(mat,v)-1

    if numKids[v] > 0
        kidsPtr[v] = numIn+1
    end

    for ind in mat.colptr[v]:(mat.colptr[v+1]-1)
      nbr = mat.rowval[ind]
      if !visited[nbr]
          numIn = numIn+1
          visited[nbr] = true
          children[numIn] = nbr
          weights[numIn] = mat.nzval[ind]

          parent[nbr] = v
          depth[nbr] = depth[v] + 1/mat.nzval[ind]

      end # if
    end # for

    ptr = ptr+1

  end # while

  if numIn < n
      error("graph is not connected")
  end

  return RootedTree(root, parent, children, weights, numKids, kidsPtr), depth

end

import Base.pop!

# put the nodes of a tree in a dfs order,
# so the last node we see is a leaf,
# the first is the root,
# and each subtree is completed before a parallel subtree
#
# we need this for the non-recursive tarjanLCA
function dfsOrder(tr::RootedTree)
    ord = Array{Int64}(0)
    stk = Array{Int64}(0)
    push!(stk,tr.root)
    while ~isempty(stk)
        u = pop!(stk)
        push!(ord,u)
        for  vi in tr.kidsPtr[u]:(tr.numKids[u] + tr.kidsPtr[u] - 1)
            v = tr.children[vi]
            push!(stk,v)
        end
    end
    return ord

end

# put the nodes of a tree in a dfs order,
# so the last node we see is a leaf,
# the first is the root,
# and each subtree is completed before a parallel subtree
#
# this matrix should be the adj matrix of a tree
function dfsOrder{Tv,Ti}(t::SparseMatrixCSC{Tv,Ti}; start::Ti = 1)

    n = size(t,1)
    seen = zeros(Bool,n)
    ord = Array{Int64}(0)
    stk = Array{Int64}(0)
    push!(stk,start)
    while ~isempty(stk)
        u = pop!(stk)
        push!(ord,u)
        seen[u] = true
        for vi in t.colptr[u]:(t.colptr[u+1]-1)
            v = t.rowval[vi]
            if !seen[v]
                push!(stk,v)
            end
        end
    end
    return ord

end

# put the nodes of a tree in a dfs order,
# so the first node we see is a leaf,
# the last is the root,
# and each subtree is completed before a parallel subtree
#
# we need this for the non-recursive tarjanLCA
#=
function dfsOrder(t::RootedTree)

    n = size(t.children,1)

    sz = ones(Int64,n);
    for vi in n:-1:2
        v = t.children[vi]
        par = t.parent[v]
        sz[par] += sz[v]
    end

    numLeft = zeros(Int64,n)
    for ui in 1:n
        u = t.children[ui]
        cnt = numLeft[u]
        for vi in t.kidsPtr[u]:(t.numKids[u] + t.kidsPtr[u] - 1)
            v = t.children[vi]
            numLeft[v] = cnt
            cnt += sz[v]
        end
    end

    ordInd = numLeft + sz

    ord = zeros(Int64,n)
    for i in 1:n
        ord[ordInd[i]] = i
    end

    return ord

end
=#

# this is intended to be used with a tree
function bfsOrder{Tv,Ti}(mat::SparseMatrixCSC{Tv,Ti}, start::Ti)
  n = mat.n
  visited = zeros(Bool,n)

  order = zeros(Ti,n)
  order[1] = start
  ptr::Ti = 1
  numInOrder::Ti = 1
  visited[start] = true

  while ptr <= numInOrder
    v = order[ptr]

    for ind in mat.colptr[v]:(mat.colptr[v+1]-1)
      nbr = mat.rowval[ind]
      if !visited[nbr]
          numInOrder = numInOrder+1
          order[numInOrder] = nbr
          visited[nbr] = true
      end # if
    end # for

    ptr = ptr+1

  end # while

  if numInOrder < n
      error("graph is not connected")
  end

  return order

end # bfsOrder

#-----------------------------
#
# code for computing stretch

# compute the depth of every node in a rooted tree,
# taking the weight as the reciprocal of distance
function compDepth{Tv,Ti}(t::RootedTree{Tv,Ti})
    n = length(t.parent)
    depth = zeros(Tv,n)

    depth[t.root] = 0

    for v in t.children
        for i in t.kidsPtr[v]:(t.numKids[v] + t.kidsPtr[v] - 1)
            depth[t.children[i]] = depth[v] + 1/t.weights[i]
        end
    end

    return depth
end # compDepth

# tarjan stretch uses tarjans offline lca algorithm to compute stretches

# it seems that this cannot be made much faster:
# almost all of the time is in the disjoint sets code.
# and, when I just run that without the rest, it is most of the time.
# this will save time by computing stretch directly,
# by using the depth
function tarjanStretch{Tv,Ti}(t::RootedTree{Tv,Ti}, mat::SparseMatrixCSC{Tv,Ti}, depth::Array{Tv,1})
    n = length(t.parent)
    su = IntDisjointSets(n)

    ancestor = collect(1:n)

    answer = zeros(Tv,nnz(mat))

    seen = zeros(Bool, n)

    tarjanStretchSub(t.root, t, mat, ancestor, answer, seen, su, depth)

    stretches = copy(mat)
    stretches.nzval = answer
    stretches = stretches + stretches'

    return stretches

end # tarjanStretch

function tarjanStretchSub{Tv,Ti}(u::Ti, t::RootedTree{Tv,Ti}, mat::SparseMatrixCSC{Tv,Ti}, ancestor::Array{Ti,1},
                                 answer::Array{Tv,1}, seen::Array{Bool,1}, su::IntDisjointSets, depth::Array{Tv,1})

    ord = dfsOrder(t)

    # traverse nodes from leaves back to root
    for vi in size(ord,1):-1:1
        v = ord[vi]

        par = t.parent[v]

        # just for debugging
        if seen[par]
            error("saw parent!")
        end

        seen[v] = true

        # ancestor[DataStructures.find_root(su, v)] = par

        for ind in mat.colptr[v]:(mat.colptr[v+1]-1)
            w = mat.rowval[ind]
            if seen[w]
                answer[ind] = mat.nzval[ind]*(depth[v] + depth[w] - 2*depth[ancestor[DataStructures.find_root(su,w)]])
                # println(u, " ", v, " : ", answer[ind])
            end # can fill u-v query
        end # over queries

        DataStructures.union!(su, par, v)

        ancestor[DataStructures.find_root(su, par)] = par

     end # for vi, v

end # TarjanStretchSub

function compStretches{Tv,Ti}(t::RootedTree{Tv,Ti}, mat::SparseMatrixCSC{Tv,Ti})
    n = length(t.order)

    depth = compDepth(t)

    stretches = tarjanStretch(t,mat,depth)
    return stretches

end # compStretches

"""Compute the stretched of every edge in `mat` with respect to the tree `tree`.
Returns the answer as a sparse matrix with the same nonzero structure as `mat`.
Assumes that `mat` is symmetric.
`tree` should be the adjacency matrix of a spanning tree."""
function compStretches{Tv,Ti}(tree::SparseMatrixCSC{Tv,Ti}, mat::SparseMatrixCSC{Tv,Ti})

    t, depth = matToTreeDepth(tree)

    stretches = tarjanStretch(t,mat,depth)
    return stretches

end # compStretches

"""Compute the vector of depths in a tree that is in DFS order,
*with the root at the first position, and the leaves at the end*
"""
function treeDepthDFS{Tv,Ti}(tree::SparseMatrixCSC{Tv,Ti})
    n = tree.n

    depth = zeros(Tv,n)
    root = 1
    kids = tree.rowval[tree.colptr[root]:(tree.colptr[root+1]-1)]
    wts =  tree.nzval[tree.colptr[root]:(tree.colptr[root+1]-1)]
    depth[kids] = 1./wts

    for v in 2:n

        for ind in (tree.colptr[v]+1):(tree.colptr[v+1]-1)
            kid = tree.rowval[ind]
            depth[kid] = depth[v] + 1/tree.nzval[ind]
        end
    end

    return depth
end

"""Compute the stretched of every edge in `mat` with respect to the tree `tree`.
Returns the answer as a sparse matrix with the same nonzero structure as `mat`.
Assumes that `mat` is symmetric.
`tree` should be the adjacency matrix of a spanning tree,
*ordered by DFS so that every parent comes before its children in the order*"""
function compStretchesDFS{Tv,Ti}(tree::SparseMatrixCSC{Tv,Ti}, mat::SparseMatrixCSC{Tv,Ti})

    n = size(tree,1)
    su = IntDisjointSets(n)

    ancestor = collect(1:n)
    answer = zeros(Tv,nnz(mat))
    seen = zeros(Bool, n)

    depth = treeDepthDFS(tree)

    # traverse nodes from leaves back to root
    for v in n:-1:1

        if v > 1
            par = tree.rowval[tree.colptr[v]]
        else
            par = 1
        end

        # just for debugging
        if seen[par]
            error("saw parent!")
        end

        seen[v] = true

        for ind in mat.colptr[v]:(mat.colptr[v+1]-1)
            w = mat.rowval[ind]
            if seen[w]
                answer[ind] = mat.nzval[ind]*(depth[v] + depth[w] - 2*depth[ancestor[DataStructures.find_root(su,w)]])
            end # can fill u-v query
        end # over queries

        DataStructures.union!(su, par, v)

        ancestor[DataStructures.find_root(su, par)] = par

    end # for v

    stretches = copy(mat)
    stretches.nzval = answer
    stretches = stretches + stretches'

    return stretches

end # compStretchesDFS

#-----------------------------------
#
# extra stuff we don't need

# this is tarjans offline lca algorithm
# we don't use it.
#
# it seems that this cannot be made much faster:
# almost all of the time is in the disjoint sets code.
# and, when I just run that without the rest, it is most of the time.
function tarjanOLCA{Tv,Ti}(t::RootedTree{Tv,Ti}, mat::SparseMatrixCSC{Tv,Ti})
    n = length(t.order)
#    su = SetUnion(n)
    su = IntDisjointSets(n)

    ancestor = zeros(Ti,n)

    answer = zeros(Ti,nnz(mat))

    seen = zeros(Bool, n)

    tarjanSub(t.root, t, mat, ancestor, answer, seen, su)

    ansmat = copy(mat)
    ansmat.nzval = answer
    ansmat = ansmat + ansmat'

    return ansmat

end # tarjanOLCA

function tarjanSub{Tv,Ti}(u::Ti, t::RootedTree{Tv,Ti}, mat::SparseMatrixCSC{Tv,Ti}, ancestor::Array{Ti,1}, answer::Array{Ti,1}, seen::Array{Bool,1}, su::IntDisjointSets)

        ancestor[u] = u

        for v in t.children[u]
            tarjanSub(v, t, mat, ancestor, answer, seen, su)
            DataStructures.union!(su, u, v)
            ancestor[DataStructures.find_root(su, u)] = u
#            ancestor[DataStructures.find!(su, u)] = u
        end

        seen[u] = true

        for ind in mat.colptr[u]:(mat.colptr[u+1]-1)
            v = mat.rowval[ind]
            if seen[v]
#                answer[ind] = ancestor[find!(su,v)]
                answer[ind] = ancestor[DataStructures.find_root(su,v)]
                # println(u, " ", v, " : ", answer[ind])
            end # can fill u-v query
        end # over queries

end # TarjanOLCA

# union-find code, modified from https://en.wikipedia.org/wiki/Tarjan%27s_off-line_lowest_common_ancestors_algorithm

#=
type SetUnion{Ti}
  parent::Array{Ti,1}
  rank::Array{Ti,1}
end # setUnion

SetUnion(n::Int32) = SetUnion([1:n], zeros(Int32,n))
SetUnion(n::Int64) = SetUnion([1:n], zeros(Int64,n))

# recursive path compression
function find!{Ti}(su::SetUnion{Ti}, x::Ti)
     if su.parent[x] == x
        return x
     else
        su.parent[x] = find!(su,su.parent[x])
        return su.parent[x]
     end
end

function union!{Ti}(su::SetUnion{Ti}, x::Ti, y::Ti)
     xroot = find!(su, x)
     yroot = find!(su, y)
     if su.rank[xroot] > su.rank[yroot]
         su.parent[yroot] = xroot
     elseif su.rank[xroot] < su.rank[yroot]
         su.parent[xroot] = yroot
     elseif xroot != yroot
         su.parent[yroot] = xroot
         su.rank[xroot] = su.rank[xroot] + 1
     end
end

=#


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

include("../perfutil.jl")

abstract type List{T} end

mutable struct Nil{T} <: List{T}
end

mutable struct Cons{T} <: List{T}
    head::T
    tail::List{T}
end

function listn1n2(n1::Int,n2::Int)
    l1 = Nil{Int}()
    for i=n2:-1:n1
        l1 = Cons{Int}(i,l1)
    end
    l1
end

@timeit listn1n2(1,10^6) "cons" "List concatenation"
gc()

# issue #1211
include("ziggurat.jl")
a = Array{Float64}(1000000)
@timeit randn_zig!(a) "randn_zig" "Ziggurat gaussian number generator"

# issue #950
include("gk.jl")
@timeit gk(350,[0.1]) "gk" "Grigoriadis Khachiyan matrix games"

# issue #942
s = sparse(ones(280,280))
@timeit s*s "sparsemul" "Sparse matrix - sparse matrix multiplication"
s2 = sparse(rand(1:2000,10^5), kron([1:10^4;],ones(Int,10)), ones(Int,10^5), 2000, 10^4)
@timeit s2*s2' "sparsemul2" "Sparse matrix - matrix multiplication with fill-in"

# issue #938
x = 1:600000
@timeit sparse(x,x,x) "sparserange" "Construction of a sparse array from ranges"

# issue 4707
include("getdivgrad.jl")
A = getDivGrad(64,64,64)
v = rand(64^3)
@timeit A*v "matvec" "Sparse matrix - dense vector multiplication"

# issue #939
y = [500000:-1:1;]
@timeit sortperm(y) "sortperm" "Sorting of a worst-case vector"

# issue #445
include("stockcorr.jl")
@timeit stockcorr() "stockcorr" "Correlation analysis of random matrices"

include("bench_eu.jl")
@timeit bench_eu_vec(10000) "bench_eu_vec" "Vectorized Monte Carlo financial simulation"
@timeit bench_eu_devec(10000) "bench_eu_devec" "Devectorized Monte Carlo financial simulation"

# issue #1163
include("actor_centrality.jl")
@timeit actor_centrality() "actorgraph" "Centrality of actors in IMDB database"

# issue #1168
include("laplace.jl")
@timeit laplace_vec() "laplace_vec" "Vectorized Laplacian"
@timeit laplace_devec() "laplace_devec" "Devectorized Laplacian"

# issue #1169
include("go_benchmark.jl")
@timeit benchmark(10) "go_benchmark" "Simulation of random games of Go"

# issue #3142
include("simplex.jl")
@timeit doTwoPassRatioTest() "simplex" "Dual simplex algorithm for Linear Programming"

# issue #3811
include("raytracer.jl")
@timeit Raytracer(5, 256, 4) "raytracer" "raytracer"

function cmp_with_func(x::Vector, f::Function)
    count::Int = 0
    for i = 1:length(x)-1
      if f(x[i], x[i+1]) count += 1 end
    end
    count
end

x = randn(200_000)
@timeit (for n in 1:10; count = cmp_with_func(x, isless) end) "funarg" "Function argument benchmark"

arith_vectorized(b,c,d) = b.*c + d .+ 1.0

len = 1_000_000
b = randn(len)
c = randn(len)
d = randn(len)

@timeit (for n in 1:10; a = arith_vectorized(b,c,d); end) "vectorize" "Vectorized arithmetic"

writedlm("random.csv", rand(100000, 4), ',')

function parsecsv()
    for line in eachline("random.csv")
        split(line, ',')
    end
end

@timeit parsecsv() "splitline" "CSV parsing"

rm("random.csv")

include("json.jl")

_json_data = "{\"web-app\": {
  \"servlet\": [
    {
      \"servlet-name\": \"cofaxCDS\",
      \"servlet-class\": \"org.cofax.cds.CDSServlet\",
      \"init-param\": {
        \"configGlossary:installationAt\": \"Philadelphia, PA\",
        \"configGlossary:adminEmail\": \"ksm@pobox.com\",
        \"configGlossary:poweredBy\": \"Cofax\",
        \"configGlossary:poweredByIcon\": \"/images/cofax.gif\",
        \"configGlossary:staticPath\": \"/content/static\",
        \"templateProcessorClass\": \"org.cofax.WysiwygTemplate\",
        \"templateLoaderClass\": \"org.cofax.FilesTemplateLoader\",
        \"templatePath\": \"templates\",
        \"templateOverridePath\": \"\",
        \"defaultListTemplate\": \"listTemplate.htm\",
        \"defaultFileTemplate\": \"articleTemplate.htm\",
        \"useJSP\": false,
        \"jspListTemplate\": \"listTemplate.jsp\",
        \"jspFileTemplate\": \"articleTemplate.jsp\",
        \"cachePackageTagsTrack\": 200,
        \"cachePackageTagsStore\": 200,
        \"cachePackageTagsRefresh\": 60,
        \"cacheTemplatesTrack\": 100,
        \"cacheTemplatesStore\": 50,
        \"cacheTemplatesRefresh\": 15,
        \"cachePagesTrack\": 200,
        \"cachePagesStore\": 100,
        \"cachePagesRefresh\": 10,
        \"cachePagesDirtyRead\": 10,
        \"searchEngineListTemplate\": \"forSearchEnginesList.htm\",
        \"searchEngineFileTemplate\": \"forSearchEngines.htm\",
        \"searchEngineRobotsDb\": \"WEB-INF/robots.db\",
        \"useDataStore\": true,
        \"dataStoreClass\": \"org.cofax.SqlDataStore\",
        \"redirectionClass\": \"org.cofax.SqlRedirection\",
        \"dataStoreName\": \"cofax\",
        \"dataStoreDriver\": \"com.microsoft.jdbc.sqlserver.SQLServerDriver\",
        \"dataStoreUrl\": \"jdbc:microsoft:sqlserver://LOCALHOST:1433;DatabaseName=goon\",
        \"dataStoreUser\": \"sa\",
        \"dataStorePassword\": \"dataStoreTestQuery\",
        \"dataStoreTestQuery\": \"SET NOCOUNT ON;select test='test';\",
        \"dataStoreLogFile\": \"/usr/local/tomcat/logs/datastore.log\",
        \"dataStoreInitConns\": 10,
        \"dataStoreMaxConns\": 100,
        \"dataStoreConnUsageLimit\": 100,
        \"dataStoreLogLevel\": \"debug\",
        \"maxUrlLength\": 500}},
    {
      \"servlet-name\": \"cofaxEmail\",
      \"servlet-class\": \"org.cofax.cds.EmailServlet\",
      \"init-param\": {
      \"mailHost\": \"mail1\",
      \"mailHostOverride\": \"mail2\"}},
    {
      \"servlet-name\": \"cofaxAdmin\",
      \"servlet-class\": \"org.cofax.cds.AdminServlet\"},

    {
      \"servlet-name\": \"fileServlet\",
      \"servlet-class\": \"org.cofax.cds.FileServlet\"},
    {
      \"servlet-name\": \"cofaxTools\",
      \"servlet-class\": \"org.cofax.cms.CofaxToolsServlet\",
      \"init-param\": {
        \"templatePath\": \"toolstemplates/\",
        \"log\": 1,
        \"logLocation\": \"/usr/local/tomcat/logs/CofaxTools.log\",
        \"logMaxSize\": \"\",
        \"dataLog\": 1,
        \"dataLogLocation\": \"/usr/local/tomcat/logs/dataLog.log\",
        \"dataLogMaxSize\": \"\",
        \"removePageCache\": \"/content/admin/remove?cache=pages&id=\",
        \"removeTemplateCache\": \"/content/admin/remove?cache=templates&id=\",
        \"fileTransferFolder\": \"/usr/local/tomcat/webapps/content/fileTransferFolder\",
        \"lookInContext\": 1,
        \"adminGroupID\": 4,
        \"betaServer\": true}}],
  \"servlet-mapping\": {
    \"cofaxCDS\": \"/\",
    \"cofaxEmail\": \"/cofaxutil/aemail/*\",
    \"cofaxAdmin\": \"/admin/*\",
    \"fileServlet\": \"/static/*\",
    \"cofaxTools\": \"/tools/*\"},

  \"taglib\": {
    \"taglib-uri\": \"cofax.tld\",
    \"taglib-location\": \"/WEB-INF/tlds/cofax.tld\"}}}"

@timeit (for n in 1:10; a = parse_json(_json_data); end) "json" "JSON parsing"

include("indexing.jl")
x = [1:100_000;]
y = filter(iseven, 1:length(x))
logical_y = map(iseven, 1:length(x))

@timeit (for n=1:100 add1!(x,y) end) "add1" "Increment vector x at locations y"
@timeit (for n=1:100 devec_add1!(x,y) end) "devec_add1" "Devectorized increment vector x at locations y"
@timeit (for n=1:100 add1!(x,logical_y) end) "add1_logical" "Increment x_i if y_i is true"
@timeit (for n=1:100 devec_add1_logical!(x,logical_y) end) "devec_add1_logical" "Devectorized increment x_i if y_i is true"

maxrss("kernel")

max_stars_repo_name: xorJane/julia

max_issues_repo_name: xorJane/julia

max_forks_repo_name: xorJane/julia


text:
using Soss

multiple_NB_regression = @model (traps, live_in_super, log_sq_foot) begin
    alpha ∼ Normal(log(4), 1)
    beta ∼ Normal(-0.25, 1)
    beta_super ∼ Normal(-0.5, 1)
    inv_phi ∼ Normal(0, 1)

    complaints ∼ For(eachindex(traps)) do i
        # TODO: Log parametrization
        r = exp(alpha + beta * traps[i] + beta_super * live_in_super[i] + log_sq_foot[i])
        p = inv(inv_phi)
        # NegativeBinomial(r, p)
    end
end


text:
using JuMP
using GLPK
model = Model(GLPK.Optimizer)
n = 5
JuMP.@variables(model, begin
    y[1:n] >= 0
end)
rtg1_idx = [1,2,4]
rtg2_idx = [3]
rtg3_idx = [5]

a = [@constraint(model,sum(y[rtg1_idx])==12),
     @constraint(model,sum(y[rtg2_idx])==12),
     @constraint(model,sum(y[rtg3_idx])==12)
]

set_name(a[1], "efefef")
set_name(a[2], "efef")
set_name(a[3], "efefe")

optimize!(model)

model.obj_dict
value.(y)


text:
using LinearAlgebra, ForwardDiff
using Convex, SCS, ECOS

"""
	projection onto second-order cone
"""
function Πsoc(v, s)
	if norm(v) <= -s
		# @warn "below cone"
		return zero(v), 0.0
	elseif norm(v) <= s
		# @warn "in cone"
		return v, s
	elseif norm(v) > abs(s)
		# @warn "outside cone"
		a = 0.5 * (1.0 + s / norm(v))
		return a * v, a * norm(v)
	else
		@warn "soc projection error"
		return zero(v), 0.0
	end
end

function projection_difference(v, s)
	vcat(Πsoc(v, s)...) - [v; s]
end

# projection_difference(ones(2), 1.0)

"""
 min v'b
 st norm(b) <= y
"""

# parameters
v = [1.0; 1.0e-1]
ψ = norm(v)
y = 1.0

"Convex.jl"
b = Variable(2)
prob = minimize(v' * b)
prob.constraints += norm(b) <= y
@time solve!(prob, ECOS.Optimizer)

@show prob.status
@show b.value
@show prob.constraints[1].dual
prob.optval

function lagrangian(v, y, b, ψ)
	v' * b + ψ' * projection_difference(b, y) + 0.5 * projection_difference(b, y)' * projection_difference(b, y)
end

function r(z, θ)
	b = z[1:2]
	ψ = z[3:5]
	y = θ[1]
	v = θ[2:3]

	lb(w) = lagrangian(v, y, w, ψ)
	# lψ(w) = lagrangian(v, y, b, w)

	# return ForwardDiff.gradient(lb, b)
	return [ForwardDiff.gradient(lb, b);
			projection_difference(b, y)]

end

θ = [y; v]
rz(x) = r(x, θ)

sol = levenberg_marquardt(rz, ones(5))
rθ(x) = r(sol, x)

drdz = ForwardDiff.jacobian(rz, sol)
drdθ = ForwardDiff.jacobian(rθ, θ)

# norm(drdz - drz(z, θ))
# norm(drdθ - drθ(z, θ))
rank(drdz)
# eigen(drdz).values
#
# rank(drz(z, θ))
# eigen(drz(z, θ)).values

x1 = (-drdz \ drdθ)[1:2,:]
x3 = (-(drdz' * drdz) \ (drdz' * drdθ))[1:2,:]

ρ = 1.0e-8
x2 = (-(drdz' * drdz + ρ * I) \ (drdz' * drdθ))[1:2,:]
x4 = (-drdz' * ((drdz * drdz' + ρ * I) \ drdθ))[1:2,:]

# x2 = (-drz(z,θ) \ drθ(z,θ))[1:2,:]
# x4 = (-(drz(z,θ)' * drz(z,θ)) \ (drz(z, θ)' * drθ(z, θ)))[1:2, :]
#


text:
"""
    ProxAugmentedNLPBlockModel{T, S} <: AbstractNLPModel{T, S}
A data type to store augmented subproblems.
"""
mutable struct ProxAugmentedNLPBlockModel{T,S} <: AbstractNLPModel{T,S}
    meta::NLPModelMeta{T,S}
    counters::Counters
    subproblem::AbstractNLPModel
    λ::AbstractVector # dual variables
    ρ::Number # penalty Parameter
    A::AbstractMatrix # linking matrix
    b::AbstractVector
    sol::AbstractVector # current primal solution
    block_sol::AbstractVector # current block primal solution for the proximal term
    hess_info::AugmentedHessianInfo # precompute to save computation effort
    P::AbstractMatrix # proximal term's penalty parameters
end

"""
    ProxAugmentedNLPBlockModel(
        nlp::AbstractNLPModel,
        λ::AbstractVector,
        ρ::Number,
        A::AbstractMatrix,
        b::AbstractVector,
        sol::AbstractVector,
        P::AbstractMatrix
    )
Modifies a subproblem by penalizing and dualizing the linking constraints.
# Arguments
- `nlp::AbstractNLPModel`: the subproblem
- `λ::AbstractVector`: vector of dual variables
- `ρ::Number`: penalty parameter
- `A::AbstractMatrix`: full linking matrix
- `b::AbstractVector`: RHS vector for linking constraints
- `sol::AbstractVector`: vector of primal variables
- `P::AbstractMatrix`: matrix of proximal term's penalty parameters
"""
function ProxAugmentedNLPBlockModel(
    nlp::AbstractNLPModel,
    λ::AbstractVector,
    ρ::Number,
    A::AbstractMatrix,
    b::AbstractVector,
    sol::AbstractVector,
    P::AbstractMatrix,
)
    ATA = findnz(ρ .* A[:, nlp.var_idx]' * A[:, nlp.var_idx])
    block_hess_struct = hess_structure(nlp.problem_block)
    P_triplets = findnz(P)

    I = vcat(block_hess_struct[1], ATA[1], P_triplets[1])
    J = vcat(block_hess_struct[2], ATA[2], P_triplets[2])
    V = vcat(ones(nlp.meta.nnzh), ρ .* abs.(ATA[3]), abs.(P_triplets[3]))
    prox_aug_hess_struct = findnz(sparse(I, J, V, nlp.meta.nvar, nlp.meta.nvar))
    meta = update_nnzh(nlp.meta, length(prox_aug_hess_struct[1]))

    hess_struct = AugmentedHessianInfo(
        (prox_aug_hess_struct[1], prox_aug_hess_struct[2]),
        (block_hess_struct[1], block_hess_struct[2]),
        ATA,
    )
    return ProxAugmentedNLPBlockModel(
        meta,
        Counters(),
        nlp,
        λ,
        ρ,
        A,
        b,
        deepcopy(sol),
        sol[nlp.var_idx], # for the proximal term
        hess_struct,
        P,
    )
end

"""
    update_primal!(
      nlp::AugmentedNLPBlockModel,
      sol::AbstractVector
    )
Updates the primal solution estimate for the augmented nlp block `nlp`.
# Arguments
- `nlp::AugmentedNLPBlockModel`: the subproblem
- `sol::AbstractVector`: vector of primal variables
"""
function update_primal!(nlp::ProxAugmentedNLPBlockModel, sol::AbstractVector)
    nlp.sol .= sol
    nlp.block_sol .= sol[nlp.subproblem.var_idx]
end

"""
    update_dual!(
      nlp::AugmentedNLPBlockModel,
      λ::AbstractVector,
    )
Updates the dual solution estimate in-place for the augmented nlp block `nlp`.
# Arguments
- `nlp::AugmentedNLPBlockModel`: the subproblem
- `λ::AbstractVector`: vector of dual variables
"""
function update_dual!(nlp::ProxAugmentedNLPBlockModel, λ::AbstractVector)
    nlp.λ .= λ
end

function NLPModels.obj(nlp::ProxAugmentedNLPBlockModel, x::AbstractVector)
    nlp.sol[nlp.subproblem.var_idx] = x
    return obj(nlp.subproblem.problem_block, x) +
           dot(nlp.λ, nlp.A[:, nlp.subproblem.var_idx], x) +
           (nlp.ρ / 2) * norm(nlp.A * nlp.sol - nlp.b)^2 +
           1 / 2 * dot(
               (x - nlp.block_sol),
               nlp.P,
               (x - nlp.block_sol),
           )
end

function NLPModels.grad!(
    nlp::ProxAugmentedNLPBlockModel,
    x::AbstractVector,
    g::AbstractVector,
)
    nlp.sol[nlp.subproblem.var_idx] = x

    grad!(nlp.subproblem.problem_block, x, g)
    mul!(
        g,
        nlp.A[:, nlp.subproblem.var_idx]',
        (nlp.λ + nlp.ρ .* (nlp.A * nlp.sol - nlp.b)),
        1,
        1,
    )
    mul!(g, nlp.P, (x - nlp.block_sol), 1, 1)
    return g
end

function NLPModels.hess_structure!(
    nlp::ProxAugmentedNLPBlockModel,
    rows::AbstractVector{T},
    cols::AbstractVector{T},
) where {T}
    rows .= nlp.hess_info.augmented_hessian_struct[1]
    cols .= nlp.hess_info.augmented_hessian_struct[2]
    return rows, cols
end

function NLPModels.hess_coord!(
    nlp::ProxAugmentedNLPBlockModel,
    x::AbstractVector{T},
    vals::AbstractVector{T};
    obj_weight = one(T),
) where {T}
    get_augmented_hessian_coord!(nlp, x, vals, obj_weight)

    aug_hess = nlp.hess_info.augmented_hessian_struct
    P_triplets = findnz(nlp.P)
    # add the proximal term weights
    main_idx = 1
    sub_idx = 1
    for (i, j) in zip(aug_hess[1], aug_hess[2])
        if P_triplets[1][sub_idx] == i && P_triplets[2][sub_idx] == j
            vals[main_idx] += P_triplets[3][sub_idx]
            sub_idx += 1
        end
        main_idx += 1
    end
    return vals
end

function NLPModels.hess_coord!(
    nlp::ProxAugmentedNLPBlockModel,
    x::AbstractVector{T},
    y::AbstractVector{T},
    vals::AbstractVector{T};
    obj_weight = one(T),
) where {T}
    get_augmented_hessian_coord!(nlp, x, vals, obj_weight, y = y)

    aug_hess = nlp.hess_info.augmented_hessian_struct
    P_triplets = findnz(nlp.P)
    # add the proximal term weights
    main_idx = 1
    sub_idx = 1
    for (i, j) in zip(aug_hess[1], aug_hess[2])
        if P_triplets[1][sub_idx] == i && P_triplets[2][sub_idx] == j
            vals[main_idx] += P_triplets[3][sub_idx]
            sub_idx += 1
        end
        main_idx += 1
    end
    return vals
    return vals
end

function NLPModels.cons!(
    nlp::ProxAugmentedNLPBlockModel,
    x::AbstractVector,
    cx::AbstractVector,
)
    return cons!(nlp.subproblem.problem_block, x, cx)
end

function NLPModels.jac_structure!(
    nlp::ProxAugmentedNLPBlockModel,
    rows::AbstractVector{T},
    cols::AbstractVector{T},
) where {T}
    return jac_structure!(nlp.subproblem.problem_block, rows, cols)
end

function NLPModels.jac_coord!(
    nlp::ProxAugmentedNLPBlockModel,
    x::AbstractVector,
    vals::AbstractVector,
)
    return jac_coord!(nlp.subproblem.problem_block, x, vals)
end


text:
mutable struct StorageOperations{R,G,T,P}

    # Parameters

    fixedcost::Matrix{Float64}      # $/unit/investment period, g
    variablecost::Matrix{Float64}   # $/MW/hour, g

    # Variables

    energycharge::Array{VariableRef,4} # Withdraw from grid (MW, r x g x t x p)
    energydischarge::Array{VariableRef,4} # Inject to grid (MW, r x g x t x p)
    raisereserve::Array{VariableRef,4}   # Raise reserves provisioned
    lowerreserve::Array{VariableRef,4}   # Lower reserves provisioned
    startstateofcharge::Array{VariableRef,3}   # (MWh, r x g x p)

    # Expressions

    fixedcosts::Matrix{ExpressionRef} # ($, r x g)
    variablecosts::Array{ExpressionRef,3} # ($, r x g x p)
    operatingcosts::Matrix{ExpressionRef} # Operating costs ($, r x g)

    stateofcharge::Array{ExpressionRef,4}  # MWh, r x g x t x p

    ucap::Vector{ExpressionRef} # (MW, r)
    totalenergy::Array{ExpressionRef,3} # MW, r x t x p
    totalraisereserve::Array{ExpressionRef,3}
    totallowerreserve::Array{ExpressionRef,3}

    # Constraints

    mindischarge::Array{GreaterThanConstraintRef,4} # (r x g x t x p)
    maxdischarge_power::Array{LessThanConstraintRef,4}
    maxdischarge_energy::Array{LessThanConstraintRef,4}

    mincharge::Array{GreaterThanConstraintRef,4}
    maxcharge_power::Array{LessThanConstraintRef,4}
    maxcharge_energy::Array{LessThanConstraintRef,4}

    minenergy::Array{GreaterThanConstraintRef,4}   # (r x g x t x p)
    maxenergy::Array{LessThanConstraintRef,4}   # (r x g x t x p)

    minraisereserve::Array{GreaterThanConstraintRef,4}   # (r x g x t x p)
    minlowerreserve::Array{GreaterThanConstraintRef,4}   # (r x g x t x p)

    periodiccharge::Array{EqualToConstraintRef,3} # (r x g x p)

    function StorageOperations{T,P}(
        fixedcost::Matrix{Float64}, variablecost::Matrix{Float64}
    ) where {T,P}

        R, G = size(fixedcost)
        @assert size(variablecost) == (R, G)

        new{R,G,T,P}(fixedcost, variablecost)

    end

end

function StorageOperations{T,P}(
    storages::StorageDevices{G},
    regionlookup::Dict{String,Int}, storagepath::String
) where {G,T,P}

    R = length(regionlookup)
    storagedata = DataFrame!(CSV.File(joinpath(storagepath, "parameters.csv"),
                                      types=scenarios_resource_param_types))

    storagelookup = Dict(zip(storages.name, 1:G))
    fixedcost = zeros(Float64, R, G)
    variablecost = zeros(Float64, R, G)

    for row in eachrow(storagedata)
        stor_idx = storagelookup[row.class]
        region_idx = regionlookup[row.region]
        fixedcost[region_idx, stor_idx] = row.fixedcost
        variablecost[region_idx, stor_idx] = row.variablecost
    end

    return StorageOperations{T,P}(fixedcost, variablecost)

end

function setup!(
    ops::StorageOperations{R,G,T,P},
    units::StorageDevices{G},
    m::Model,
    invs::ResourceInvestments{R,G},
    periodweights::Vector{Float64},
    s::AbstractScenario
) where {R,G,T,P}

    invprob = s.investmentproblem
    Rs = invprob.regionnames
    Gs = units.name
    Ts = string.(1:T)
    Ps = invprob.periodnames

    # Variables

    ops.energydischarge = @variable(m, [1:R, 1:G, 1:T, 1:P])
    varnames!(ops.energydischarge, "energydischarge_storage_$(s.name)", Rs, Gs, Ts, Ps)

    ops.energycharge    = @variable(m, [1:R, 1:G, 1:T, 1:P])
    varnames!(ops.energycharge, "energycharge_storage_$(s.name)", Rs, Gs, Ts, Ps)

    ops.raisereserve    = @variable(m, [1:R, 1:G, 1:T, 1:P])
    varnames!(ops.raisereserve, "raisereserve_storage_$(s.name)", Rs, Gs, Ts, Ps)

    ops.lowerreserve    = @variable(m, [1:R, 1:G, 1:T, 1:P])
    varnames!(ops.lowerreserve, "lowerreserve_storage_$(s.name)", Rs, Gs, Ts, Ps)

    ops.startstateofcharge = @variable(m, [1:R, 1:G, 1:P])
    varnames!(ops.startstateofcharge, "startstateofcharge_storage_$(s.name)", Rs, Gs, Ps)

    # Expressions

    ops.fixedcosts =
        @expression(m, [r in 1:R, g in 1:G],
                    ops.fixedcost[r,g] * invs.dispatching[r,g])

    ops.variablecosts =
        @expression(m, [r in 1:R, g in 1:G, p in 1:P],
                    sum(ops.variablecost[r,g] *
                        (ops.energycharge[r,g,t,p] + ops.energydischarge[r,g,t,p])
                        for t in 1:T))

    ops.operatingcosts =
        @expression(m, [r in 1:R, g in 1:G],
                    ops.fixedcosts[r,g] +
                    sum(ops.variablecosts[r,g,p] * periodweights[p]
                        for p in 1:P))

    ops.ucap =
        @expression(m, [r in 1:R], G > 0 ? sum(
            invs.dispatching[r,g] * units.maxdischarge[g] * units.capacitycredit[g]
        for g in 1:G) : 0)

    ops.totalenergy =
        @expression(m, [r in 1:R, t in 1:T, p in 1:P],
                    G > 0 ? sum(ops.energydischarge[r,g,t,p] - ops.energycharge[r,g,t,p]
                        for g in 1:G) : 0)

    ops.totalraisereserve =
        @expression(m, [r in 1:R, t in 1:T, p in 1:P],
                    G > 0 ? sum(ops.raisereserve[r,g,t,p] for g in 1:G) : 0)

    ops.totallowerreserve =
        @expression(m, [r in 1:R, t in 1:T, p in 1:P],
                    G > 0 ? sum(ops.lowerreserve[r,g,t,p] for g in 1:G) : 0)

    ops.stateofcharge =
        @expression(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.startstateofcharge[r,g,p] +
                    sum(ops.energycharge[r,g,i,p] - ops.energydischarge[r,g,i,p]
                        for i in 1:t))

    # Constraints

    ops.mindischarge =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.energydischarge[r,g,t,p] >= 0)

    ops.maxdischarge_power =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.raisereserve[r,g,t,p] +
                    ops.energydischarge[r,g,t,p] - ops.energycharge[r,g,t,p] <=
                    invs.dispatching[r,g] * units.maxdischarge[g])

    ops.maxdischarge_energy =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.raisereserve[r,g,t,p] +
                    ops.energydischarge[r,g,t,p] - ops.energycharge[r,g,t,p] <=
                    ops.stateofcharge[r,g,t,p])

    ops.mincharge =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.energycharge[r,g,t,p] >= 0)

    ops.maxcharge_power =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.lowerreserve[r,g,t,p] +
                    ops.energycharge[r,g,t,p] - ops.energydischarge[r,g,t,p] <=
                    invs.dispatching[r,g] * units.maxcharge[g])

    ops.maxcharge_energy =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.lowerreserve[r,g,t,p] +
                    ops.energycharge[r,g,t,p] - ops.energydischarge[r,g,t,p] <=
                    invs.dispatching[r,g] * units.maxenergy[g]
                     - ops.stateofcharge[r,g,t,p])

    ops.minraisereserve =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.raisereserve[r,g,t,p] >= 0)

    ops.minlowerreserve =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.lowerreserve[r,g,t,p] >= 0)

    ops.minenergy =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.stateofcharge[r,g,t,p] >= 0)

    ops.maxenergy =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.stateofcharge[r,g,t,p] <=
                    invs.dispatching[r,g] * units.maxenergy[g])

    ops.periodiccharge =
        @constraint(m, [r in 1:R, g in 1:G, p in 1:P],
                    ops.stateofcharge[r,g,T,p] == ops.startstateofcharge[r,g,p])

end

welfare(x::StorageOperations) = -sum(x.operatingcosts)


text:
@recipe function f(pc::PCoA)
    xticks := false
    yticks := false
    xlabel --> "PCo1 ($(round(variance(pc, 1) * 100, digits = 2))%)"
    ylabel --> "PCo2 ($(round(variance(pc, 2) * 100, digits = 2))%)"
    seriestype := :scatter
    principalcoord(pc, 1), principalcoord(pc,2)
end

@userplot AbundancePlot
@recipe function f(plt::AbundancePlot; topabund=10, sort::Vector{Int}=collect(1:nsites(plt.args[1])))
    abun = plt.args[1]
    typeof(abun) <: AbstractComMatrix || error("AbundancePlot not defined for $(typeof(abun))")

    topabund = min(topabund, nfeatures(abun))
    2 <= topabund < 12 || error("n must be between 2 and 12")

    top = filterabund(abun, topabund)

    rows = featurenames(top)

    yflip := true
    bar_position := :stack
    label := rows
    GroupedBar((1:nsamples(top), Matrix(occurrences(top)[:,sort]')))
end

struct AnnotationBar
    labels::Array{<:AbstractString,1}
    colors::Array{<:Color,1}
end

annotationbar(labels::Array{<:AbstractString,1}, colors::Array{<:Color,1}) = AnnotationBar(labels, colors)
annotationbar(colors::Array{<:Color,1}) = AnnotationBar(["sample$i" for i in 1:length(colors)], colors)

@recipe function f(bar::AnnotationBar)
    xs = Int[]
    for i in 1:length(bar.colors)
        append!(xs, [0,0,1,1,0] .+ (i-1))
    end
    xs = reshape(xs, 5, length(bar.colors))
    ys = hcat([[0,1,1,0,0] for _ in bar.colors]...)

    fc = reshape(bar.colors, 1, length(bar.colors))

    seriestype := :path
    fill := (0,1)
    fillcolor := fc
    legend := false
    color --> :black
    ticks := false
    framestyle := false
    xaxis --> false
    yaxis --> false
    xs, ys
end

# From Michael K. Borregaard (posted to slack 2018-05-18)
# Usage:
#
# x = randn(100)
# y = randn(100) + 2
# y[y.<0] = 0
#
# zeroyplot(x,y, yscale = :log10, color = :red, markersize = 8)
# @userplot ZeroYPlot
# @recipe function f(h::ZeroYPlot)
#            length(h.args) != 2 && error("zeroyplot only defined for x and y input arguments")
#            x, y = h.args
#            val0 = y .== 0
#
#            layout := @layout([a
#                               b{0.05h}])
#
#            markeralpha --> 0.4
#            seriestype --> :scatter
#
#            @series begin
#                primary := false
#                subplot := 2
#                yscale := :identity
#                ylim := (-1,1)
#                yticks := ([0],[0])
#                grid := false
#                x[val0], y[val0]
#            end
#
#            subplot := 1
#            x[.!val0], y[.!val0]
# end


text:
#first add workers
using Distributed
addprocs(10)
nprocs()

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

function run_cv(n :: Int64, p :: Int64, k :: Int64, debias :: Bool, d::UnionAll, l::Link)
    #construct snpmatrix, covariate files, and true model b
    x = simulate_random_snparray(n, p, "test1.bed")
    xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
    z = ones(n, 1)

    #define true_b
    true_b = zeros(p)
    true_b[1:10] .= collect(0.1:0.1:1.0)
    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, 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)

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

    #compute cross validation
    result = @timed cv_iht(d(), l, x, z, y, 1, path, num_folds, folds=folds, use_maf=false, debias=debias, parallel=true);
    dr = result[1]
    runtime = result[2] #seconds
    memory = result[3] / 1e6  #MB

    rm("test1.bed", force=true)

    return dr, runtime, memory
end

function run()
    # set simulation parameters and seed
    repeats = 30
    n = 5000
    p = 50000
    k = 10
    d = NegativeBinomial
    l = LogLink()
    debias = false
    Random.seed!(2019)

    # initialize variables then run
    drs = zeros(20, repeats)
    run_time = zeros(repeats)
    memory = zeros(repeats)
    for i in 1:repeats
        dr, t, m = run_cv(n, p, k, debias, d, l)
        drs[:, i] .= dr
        run_time[i] = t
        memory[i] = m
    end

    writedlm("negativebinomial_nodebias_cv_drs", drs)
    writedlm("negativebinomial_nodebias_cv_run_times", run_time)
    writedlm("negativebinomial_nodebias_cv_memory", memory)
end

run()


text:
module TestVar

using Base.Test

using Bio.Seq
using Bio.Var

@testset "Counting mutations" begin

    # Create a 20bp test DNA sequence pair containing every possible transition (4),
    # every possible transversion (8), and 2 gapped sites and 2 ambiguous sites.
    # This leaves 4 sites non-mutated/conserved.
    dnas = [dna"ATTG-ACCTGGNTTTCCGAA", dna"A-ACAGAGTATACRGTCGTC"]
    m1 = seqmatrix(dnas, :seq)

    rnas = [rna"AUUG-ACCUGGNUUUCCGAA", rna"A-ACAGAGUAUACRGUCGUC"]
    m2 = seqmatrix(rnas, :seq)

    @test count_mutations(AnyMutation, dnas) == count_mutations(AnyMutation, rnas) == ([12], [16])
    @test count_mutations(AnyMutation, m1) == count_mutations(AnyMutation, m2) == ([12], [16])
    @test count_mutations(TransitionMutation, dnas) == count_mutations(TransitionMutation, rnas) == ([4], [16])
    @test count_mutations(TransitionMutation, m1) == count_mutations(TransitionMutation, m2) == ([4], [16])
    @test count_mutations(TransversionMutation, dnas) == count_mutations(TransversionMutation, rnas) == ([8], [16])
    @test count_mutations(TransversionMutation, m1) == count_mutations(TransversionMutation, m2) == ([8], [16])
    @test count_mutations(TransitionMutation, TransversionMutation, dnas) == count_mutations(TransitionMutation, TransversionMutation, rnas) == ([4], [8], [16])
    @test count_mutations(TransitionMutation, TransversionMutation, m1) == count_mutations(TransitionMutation, TransversionMutation, m2) == ([4], [8], [16])
    @test count_mutations(TransversionMutation, TransitionMutation, dnas) == count_mutations(TransversionMutation, TransitionMutation, rnas) == ([4], [8], [16])
    @test count_mutations(TransversionMutation, TransitionMutation, m1) == count_mutations(TransversionMutation, TransitionMutation, m2) == ([4], [8], [16])

    ans = Bool[false, false, true, true, false, true, true, true, false, true, true, false, true, false, true, true, false, false, true, true]
    @test flagmutations(AnyMutation, m1)[1][:,1] == ans
    @test flagmutations(AnyMutation, m2)[1][:,1] == ans

end

@testset "Distance Computation" begin

    dnas1 = [dna"ATTG-ACCTGGNTTTCCGAA", dna"A-ACAGAGTATACRGTCGTC"]
    m1 = seqmatrix(dnas1, :seq)

    dnas2 = [dna"attgaacctggntttccgaa",
             dna"atacagagtatacrgtcgtc"]
    dnas3 = [dna"attgaacctgtntttccgaa",
             dna"atagaacgtatatrgccgtc"]
    m2 = seqmatrix(dnas2, :seq)

    @test distance(Count{AnyMutation}, dnas1) == ([12], [16])
    @test distance(Count{TransitionMutation}, dnas1) == ([4], [16])
    @test distance(Count{TransversionMutation}, dnas1) == ([8], [16])
    @test distance(Count{Kimura80}, dnas1) == ([4], [8], [16])
    @test distance(Count{AnyMutation}, m1) == ([12], [16])
    @test distance(Count{TransitionMutation}, m1) == ([4], [16])
    @test distance(Count{TransversionMutation}, m1) == ([8], [16])
    @test distance(Count{Kimura80}, m1) == ([4], [8], [16])

    @test distance(Count{AnyMutation}, dnas2, 5, 5)[1][:] == [2, 4, 3, 3]
    @test distance(Count{AnyMutation}, dnas2, 5, 5)[2][:] == [5, 5, 3, 5]
    @test distance(Count{TransitionMutation}, dnas2, 5, 5)[1][:] == [0, 2, 1, 1]
    @test distance(Count{TransitionMutation}, dnas2, 5, 5)[2][:] == [5, 5, 3, 5]
    @test distance(Count{TransversionMutation}, dnas2, 5, 5)[1][:] == [2, 2, 2, 2]
    @test distance(Count{TransversionMutation}, dnas2, 5, 5)[2][:] == [5, 5, 3, 5]
    @test distance(Count{Kimura80}, dnas1) == ([4], [8], [16])

    @test distance(Count{AnyMutation}, dnas2) == ([12], [18])
    @test distance(Count{TransitionMutation}, dnas2) == ([4], [18])
    @test distance(Count{TransversionMutation}, dnas2) == ([8], [18])
    @test distance(Count{Kimura80}, dnas2) == ([4], [8], [18])
    @test distance(Count{AnyMutation}, m2) == ([12], [18])
    @test distance(Count{TransitionMutation}, m2) == ([4], [18])
    @test distance(Count{TransversionMutation}, m2) == ([8], [18])
    @test distance(Count{Kimura80}, m2) == ([4], [8], [18])

    d = distance(Proportion{AnyMutation}, dnas2, 5, 5)
    a = [0.4, 0.8, 1.0, 0.6]
    for i in 1:length(d[1])
        @test_approx_eq_eps d[1][i] a[i] 1e-4
    end
    @test d[2][:] == [5, 5, 3, 5]
    d = distance(Proportion{TransitionMutation}, dnas2, 5, 5)
    a = [0.0, 0.4, 0.333333, 0.2]
    for i in 1:length(d[1])
        @test_approx_eq_eps d[1][i] a[i] 1e-4
    end
    @test d[2][:] == [5, 5, 3, 5]
    d = distance(Proportion{TransversionMutation}, dnas2, 5, 5)
    a = [0.4, 0.4, 0.666667, 0.4]
    for i in 1:length(d[1])
        @test_approx_eq_eps d[1][i] a[i] 1e-4
    end
    @test d[2][:] == [5, 5, 3, 5]

    @test distance(Proportion{AnyMutation}, dnas1) == ([(12 / 16)], [16])
    @test distance(Proportion{TransitionMutation}, dnas1) == ([(4 / 16)], [16])
    @test distance(Proportion{TransversionMutation}, dnas1) == ([(8 / 16)], [16])
    @test distance(Proportion{AnyMutation}, m1) == ([(12 / 16)], [16])
    @test distance(Proportion{TransitionMutation}, m1) == ([(4 / 16)], [16])
    @test distance(Proportion{TransversionMutation}, m1) == ([(8 / 16)], [16])

    @test distance(Proportion{AnyMutation}, dnas2) == ([(12 / 18)], [18])
    @test distance(Proportion{TransitionMutation}, dnas2) == ([(4 / 18)], [18])
    @test distance(Proportion{TransversionMutation}, dnas2) == ([(8 / 18)], [18])
    @test distance(Proportion{AnyMutation}, m2) == ([(12 / 18)], [18])
    @test distance(Proportion{TransitionMutation}, m2) == ([(4 / 18)], [18])
    @test distance(Proportion{TransversionMutation}, m2) == ([(8 / 18)], [18])

    @test distance(JukesCantor69, dnas1) == ([Inf], [Inf]) # Returns infinity as 12/16 is 0.75 - mutation saturation.
    @test distance(JukesCantor69, m1) == ([Inf], [Inf])

    @test round(distance(JukesCantor69, dnas2)[1][1], 3) == 1.648
    @test round(distance(JukesCantor69, dnas2)[2][1], 3) == 1
    @test round(distance(JukesCantor69, m2)[1][1], 3) == 1.648
    @test round(distance(JukesCantor69, m2)[2][1], 3) == 1
    @test_throws DomainError distance(JukesCantor69, dnas2, 5, 5)
    d = distance(JukesCantor69, dnas3, 5, 5)
    a = [0.232616, 0.571605, 0.44084, 0.571605]
    v = [0.0595041, 0.220408, 0.24, 0.220408]
    for i in 1:length(d[1])
        @test_approx_eq_eps d[1][i] a[i] 1e-5
        @test_approx_eq_eps d[2][i] v[i] 1e-5
    end

    @test round(distance(Kimura80, dnas2)[1][1], 3) == 1.648
    @test round(distance(Kimura80, dnas2)[2][1], 3) == 1
    @test round(distance(Kimura80, m2)[1][1], 3) == 1.648
    @test round(distance(Kimura80, m2)[2][1], 3) == 1

end

end # module TestVar


text:
struct GeoPlot end
struct GeoZPlot end

# We only look at arrays with X, Y, Z dims here.
# Otherwise they fall back to DimensionalData.jl recipes
@recipe function f(A::AbstractGeoArray)
    ddplot(A) = DimArray(A; dims=_maybe_mapped(dims(A)))

    A = GeoArray(A)
    if !(get(plotattributes, :seriestype, :none) in (:none, :heatmap))
        DD.DimensionalPlot(), ddplot(A)
    elseif all(hasdim(A, (SpatialDim, SpatialDim)))
        # Heatmap or multiple heatmaps. Use GD recipes.
        GeoPlot(), _prepare(A)
    elseif hasdim(A, ZDim) && ndims(A) == 1
        # Z dim plot, but for spatial data we want Z on the Y axis
        GeoZPlot(), _prepare(A)
    else
        DD.DimensionalPlot(), ddplot(A)
    end
end

# Plot 3d arrays as multiple tiled plots
@recipe function f(::GeoPlot, A::GeoArray{T,3,<:Tuple{<:SpatialDim,<:SpatialDim,D}}) where {T,D}
    nplots = size(A, 3)
    if nplots > 1
        :layout --> nplots
        :title --> permutedims(string.(val(dims(A, D))))
        for i in 1:nplots
            @series begin
                seriestype := :heatmap
                aspect_ratio := 1
                subplot := i
                GeoPlot(), A[:, :, i]
                slice = A[:, :, i]
                ys, xs = map(_prepare, dims(A))
                xs, ys, parent(slice)
            end
        end
    else
        GeoPlot(), A[:, :, 1]
    end
end
# Plot a sinlge 2d map
@recipe function f(::GeoPlot, A::GeoArray{T,2,<:Tuple{<:SpatialDim,<:SpatialDim}}) where T
    # If colorbar is close to symmetric (< 25% difference) use a symmetric
    # colormap and set symmetric limits so zero shows up as a neutral color.
    A_min, A_max = extrema(skipmissing(A))
    if (A_min + A_max) / abs(A_max - A_min) < 0.25
        A_limit = max(abs(A_min), abs(A_max))
        clims = (-A_limit, A_limit)
        :seriescolor --> :balance
    else
        clims = A_min, A_max
    end

    yguide, xguide = label(dims(A))

    :seriestype --> :heatmap
    :title --> "$(_maybename(A)) $(DD.refdims_title(A))"
    :xguide --> xguide
    :yguide --> yguide
    :clims --> clims
    :colorbar_title --> name(A)

    ys, xs = map(_prepare, dims(A))

    if get(plotattributes, :seriestype, :none) == :contourf
        :linewidth --> 0
        :levels --> range(clims[1], clims[2], length=20)
        xs, ys, clamp.(A, clims[1], clims[2])
    else
        xs, ys, parent(A)
    end
end

# # Plot a vertical 1d line
@recipe function f(::GeoZPlot, A::GeoArray{T,1}) where T
    z_dim = dims(A, ZDim)
    yguide = label(z_dim)
    xguide = label(A)
    :title --> "$(_maybename(A))$(DD.refdims_title(A))"
    :xguide --> xguide
    :yguide --> yguide
    :label --> ""
    z = map(_prepare, dims(A))
    parent(A), z
end

# Plots heatmaps pixels are centered.
# So we should center, and use the projected value.
_prepare(d::Dimension) = d |> _maybe_shift |> _maybe_mapped |> index
# Convert arrays to a consistent missing value and Forward array order
_prepare(A::AbstractGeoArray) =
    _maybe_replace_missing(A) |>
    A -> reorder(A, ForwardIndex) |>
    A -> reorder(A, ForwardRelation) |>
    A -> permutedims(A, DD.commondims(>:, (ZDim, YDim, XDim, TimeDim, Dimension), dims(A)))

_maybename(A) = _maybename(name(A))
_maybename(n::Symbol) = n == Symbol("") ? "" : string(n, ": ")
_maybename(n::Name{N}) where N = _maybename(N)
_maybename(n::NoName) = ""

_maybe_replace_missing(A::AbstractArray{<:AbstractFloat}) = replace_missing(A, eltype(A)(NaN))
_maybe_replace_missing(A) = A

_maybe_shift(d) = _maybe_shift(sampling(d), d)
_maybe_shift(::Intervals, d) = DD.maybeshiftlocus(Center(), d)
_maybe_shift(sampling, d) = d

_maybe_mapped(dims::Tuple) = map(_maybe_mapped, dims)
_maybe_mapped(dim::Dimension) = _maybe_mapped(mode(dim), dim)
_maybe_mapped(mode::IndexMode, dim::Dimension) = dim
_maybe_mapped(mode::Projected, dim::Dimension) = _maybe_mapped(mappedcrs(mode), dim)
_maybe_mapped(::Nothing, dim::Dimension) = dim
_maybe_mapped(::GeoFormat, dim::Dimension) = convertmode(Mapped, dim)


text:
type LBFGSSolver <: RegressionSolver end

function solve(model::RegressionModel, method::RegERM, ::LBFGSSolver, X::AbstractMatrix, y::AbstractVector)
    function tloss_grad(theta::Vector)
      n = size(X,1)

      grad_model = gradient(model.f, X[1,:], theta)
      grad_loss = derivs(loss(method), values(model, X[1,:], theta), [y[1]])

      total = broadcast(*, grad_loss',grad_model)
      for i = 2:n
        grad_model = gradient(model.f, X[i,:], theta)
        grad_loss = derivs(loss(method), values(model, X[i,:], theta), [y[i]])

        total += broadcast(*, grad_loss',grad_model)
      end
      return vec(total)
    end

    reg_grad(theta::Vector) = gradient(regularizer(method, theta))
    obj(theta::Vector) = tloss(loss(method), values(model, method.X, theta), method.y) + value(regularizer(method, theta))
    grad!(theta::Vector, storage::Vector) = storage[:] = tloss_grad(theta) + reg_grad(theta)

    Optim.optimize(obj, grad!, model.theta, method=:l_bfgs, linesearch! = Optim.interpolating_linesearch!).minimum
end


text:
@inline function opt_diag(A::AbstractArray, B::AbstractArray=A)
    return vec(sum(A .* B; dims=2))
end

@inline function opt_trace(A::AbstractMatrix{<:Real}, B::AbstractMatrix{<:Real})
    return dot(A, B)
end

XXt(X) = X * X'

PDMats.invquad(a::Symmetric, b::AbstractVector) = dot(b, a \ b)

function ParameterHandling.flatten(L::LowerTriangular)
    reduce()
end


text:
# ------------------------------------------------------------------------------------------
# # Arrays and Loops
#
# We'll cover:
# 1. Array literals
# 2. Concatenation
# 3. For loops
# 4. Comprehensions
# 5. Element types
# 6. Dequeues
# 7. The `bang!` convention
# 8. Variable names vs. copies
#
# # Arrays
#
# Julia has highly efficient multidimensional arrays, both constructed and indexed with
# square brackets.
#
# Syntax: <br>
# ```julia
# [item1, item2, ...]
# ```
# ------------------------------------------------------------------------------------------

squares = [1, 4, 9, 15, 25, 36, 49, 64]

squares[1]

squares[1:3]

squares[end]

squares[4] = 16

squares

typeof(squares)

# ------------------------------------------------------------------------------------------
# ## Concatenation
#
# If, instead of commas, you just use spaces, then the values are concatenated horizontally.
# ------------------------------------------------------------------------------------------

cubes = [1, 8, 27, 64, 125, 216, 343, 512]

powers = [1:8 squares cubes]

powers[4, 2]

powers[:, 3]

powers[7, :]

typeof(powers)

# ------------------------------------------------------------------------------------------
# Semicolon separators perform vertical concatenation:
# ------------------------------------------------------------------------------------------

[squares; cubes]

# ------------------------------------------------------------------------------------------
# Whereas commas would simply create an array of arrays:
# ------------------------------------------------------------------------------------------

nested_powers = [[1,2,3,4,5,6,7,8], squares, cubes]

nested_powers[2]

# ------------------------------------------------------------------------------------------
# Horizontal and vertical concatenation can be used together to as a simple syntax for
# matrix literals:
# ------------------------------------------------------------------------------------------

[1 3 5; 2 4 6]

# ------------------------------------------------------------------------------------------
# # Loops
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Of course, we could construct this programmatically with a for-loop. The syntax for a
# `for` loop is
#
# ```julia
# for *var* in *loop iterable*
#     *loop body*
# end
# ```
#
#
# ------------------------------------------------------------------------------------------

A = fill(0, (8, 3)) # Allocate an 8x3 matrix to store the values into
for pow in 1:3
    for value in 1:8
        A[value, pow] = value ^ pow
    end
end
A

A == powers

# ------------------------------------------------------------------------------------------
# ## Array Comprehensions
# ------------------------------------------------------------------------------------------

squares = [value^2 for value in 1:8]

cubes = [value^3 for value in 1:8]

powers = [value^pow for value in 1:8, pow in 1:3]

# ------------------------------------------------------------------------------------------
# # The element type
#
# Note that every time an array prints out, it is displaying its element type and
# dimensionality, for example `Array{Int64, 2}`. This describes what it can store — and thus
# what it can return upon indexing.
# ------------------------------------------------------------------------------------------

typeof(powers)

typeof(powers[1, 1])

# ------------------------------------------------------------------------------------------
# Further, the array will try to convert any new values assigned into it to its element
# type:
# ------------------------------------------------------------------------------------------

powers[1, 1] = 1.6

powers[1, 1] = -5.0 # This can be losslessly converted to an integer

powers

# ------------------------------------------------------------------------------------------
# Arrays that have an exact and concrete element type are generally significantly faster, so
# Julia will try to find an amenable element type for you in its literal construction
# syntax:
# ------------------------------------------------------------------------------------------

fortytwosarray = [42, 42.0, 4.20e1, 4.20f1, 84//2, 0x2a]

for x in fortytwosarray
    show(x)
    println("\tisa $(typeof(x))")
end

# ------------------------------------------------------------------------------------------
# The `Any` array can be helpful for disabling these behaviors and allowing all kinds of
# different objects:
# ------------------------------------------------------------------------------------------

anyfortytwos = Any[42, 42.0, 4.20e1, 4.20f1, 84//2, 0x2a]

anyfortytwos[1] = "FORTY TWO"
anyfortytwos

# ------------------------------------------------------------------------------------------
# # Vectors as dequeues
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# One-dimensional arrays can be appended to and have items removed from them:
# ------------------------------------------------------------------------------------------

fib = [1, 1, 2, 3, 5, 8, 13]

push!(fib, 21)

push!(fib, sum(fib[end-1:end]))

pop!(fib)

pushfirst!(fib, 0)

popfirst!(fib)

# ------------------------------------------------------------------------------------------
# ## Aside: why so shouty!?
#
# Why are there exclamations in all the above function calls? `push!(fib, ...)`,
# `pop!(fib)`, etc.?
# ------------------------------------------------------------------------------------------

push!

# ------------------------------------------------------------------------------------------
# This is entirely a convention. It's a signal to the caller that one of the argument is
# going to be _mutated_. It's perhaps easiest to demonstrate with an example:
# ------------------------------------------------------------------------------------------

A = rand(0:10, 10)

sort(A)

A

sort!(A)

A # That changed A!

# ------------------------------------------------------------------------------------------
# ## Aside: names vs. copies
#
# Remember that variables are just names we give our objects. So watch what happens if we
# give the same object two different names:
# ------------------------------------------------------------------------------------------

fibonacci = [1, 1, 2, 3, 5, 8, 13]

some_numbers = fibonacci
some_numbers[1] = 404
some_numbers

fibonacci

fibonacci[1] = 1
some_numbers = copy(fibonacci)
some_numbers[2] = 404
fibonacci

max_stars_repo_path: .nbexports/introductory-tutorials/intro-to-julia/long-version/020 Arrays and Loops.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
include("poly_factorization_project.jl")

"""
Takes some input polynomial (of type Polynomial) and converts it to a string whilst also simplifying it to make it more human readable.
"""
function pretty_printer(input::Polynomial)::String
    manipulated_string = replace(string(input), "+ -" => "- ")
    manipulated_string = replace(manipulated_string, "x^1" => "x")
    return replace(manipulated_string, "⋅x^0" => "")
end

x = x_poly()

p1 = 11x^2 + 52x^5
p2 = 103x^5 + 7x^3
prime = 5

println("To showcase some of the functions within this repo, we will use the two polynomials p1 = ", pretty_printer(p1), " and p2 = ", pretty_printer(p2), "\n")

println("The scripts in this repo allow you to:\n")
println("Add polynomials together:")
println("p1 + p2 = ", pretty_printer(p1+p2), "\n")

println("Multiply polynomials and raise them by a power:")
println("p1 * p2 = ", pretty_printer(p1*p2))
println("p1^3 = ", pretty_printer(p1^prime), "\n")

println("Divide one polynomial by another:")
println("p1/p2 = ", pretty_printer((p2÷p1)(prime)), "\n")

println("Find the remainder of this polynomial division")
println(rem(p1,p2)(prime), "\n")

println("Take the derivative of a polynomial:")
println("p1' * p2 = ", pretty_printer(derivative(p1*p2)))
println("p1' * p2 + p1 * p2' = ", derivative(p1)*p2 + p1*derivative(p2), "\n")

println("Find the GCD of two polynomials modulo prime:")
println("GCD(p1,p2)(mod 3) = ", pretty_printer(gcd(p1,p2,prime)), "\n")

println("Factorize polynomials:")
p = mod(p1,prime)
println("Will factor p1(mod $prime): ", p)
factorization = factor(p,prime)
println("Here is the factorization: ", factorization)

pr = mod(expand_factorization(factorization),prime)
println("Reconstructing: ", pretty_printer(pr), "\n")

println("Make a polynomial square free:")
println("The square free version of the polynomial, ", pretty_printer(p1), " is, ", pretty_printer(square_free(p1, prime)), "\n")

println("Construct a cyclotonic polynomial given a prime integer input:")
println("The following is a cyclotonic polynomial, ", pretty_printer(cyclotonic_polynomial(prime)), "\n")

println("Construct a linear monic polynomial given an integer input:")
println("The following is a linear monic polynomial generated with input $prime, ", pretty_printer(linear_monic_polynomial(prime)), "\n")

println("Make a polynomial primitive:")
println("The primitive part of the polynomial ", pretty_printer(p1), " is ", pretty_printer(prim_part(p1)(prime)), "\n")

println("Take the modulus of a polynomial:")
println(pretty_printer(p1), " modulo $prime is congruent to ", pretty_printer(mod(p1, prime)), "\n")

println("Find the degree of a polynomial:")
println("The degree of the polynomial ", pretty_printer(p1), " is equal to ", degree(p1), "\n")

println("Find the leading term of any given polynomial:")
println("The leading term of the polynomial ", pretty_printer(p1), " is equal to ", pretty_printer(Polynomial(leading(p1))), "\n")

println("Aside from the functions showcased, this software package has much more functionality! \nIt can also perform arithmetic on individual polynomial terms, manipulate polynomials \nprogrammatically with pop! and push! and much, much more!")