text:
module KissThreading

using Compat.Random: MersenneTwister
@static if VERSION >= v"0.7-"
    using Future: randjump
end

export trandjump, TRNG, tmap!, tmapreduce, getrange

@static if VERSION < v"0.7-"
    const _randjump = randjump

    function trandjump(rng = MersenneTwister(0), gpt=1)
        n = Threads.nthreads()
        # create gpt+1 RNGs, but leave gpt
        # generator number gpt+1 is dropped to have rngs for each thread separated in memory
        rngjmp = randjump(rng, n*(gpt+1))
        reshape(rngjmp, (gpt+1, n))[1:gpt, :]
    end
else
    _randjump(rng, n, jump=big(10)^20) = accumulate(randjump, [jump for i in 1:n], init = rng)

    function trandjump(rng = MersenneTwister(0); jump=big(10)^20)
        n = Threads.nthreads()
        rngjmp = Vector{MersenneTwister}(undef, n)
        Threads.@threads for i in 1:n
            rngjmp[i] = randjump(rng, big(10)^20*i)
        end
        rngjmp
    end
end

const TRNG = trandjump()

default_batch_size(n) = min(n, round(Int, 10*sqrt(n)))

struct Mapper
    atomic::Threads.Atomic{Int}
    len::Int
end

@inline function (mapper::Mapper)(batch_size, f, dst, src...)
    ld = mapper.len
    atomic = mapper.atomic
    Threads.@threads for j in 1:Threads.nthreads()
        while true
            k = Threads.atomic_add!(atomic, 1)
            batch_start = 1 + (k-1) * batch_size
            batch_end = min(k * batch_size, ld)
            batch_start > ld && break
            batch_map!(batch_start:batch_end, f, dst, src...)
        end
    end
    dst
end

@inline function batch_map!(range, f, dst, src...)
    @inbounds for j in range
        dst[j] = f(getindex.(src, j)...)
    end
end

function tmap!(f, dst::AbstractArray, src::AbstractArray...; batch_size=1)
    ld = length(dst)
    if (ld, ld) != extrema(length.(src))
        throw(ArgumentError("src and dst vectors must have the same length"))
    end
    atomic = Threads.Atomic{Int}(1)
    mapper = Mapper(atomic, ld)
    mapper(batch_size, f, dst, src...)
end

struct MapReducer{T}
    r::Base.RefValue{T}
    atomic::Threads.Atomic{Int}
    lock::Threads.SpinLock
    len::Int
end

@inline function (mapreducer::MapReducer{T})(batch_size, f, op, src...) where T
    atomic = mapreducer.atomic
    lock = mapreducer.lock
    len = mapreducer.len
    Threads.@threads for j in 1:Threads.nthreads()
        k = Threads.atomic_add!(atomic, batch_size)
        k > len && continue
        r = T(f(getindex.(src, k)...))
        range = (k + 1) : min(k + batch_size - 1, len)
        r = batch_mapreduce(r, range, f, op, src...)
        k = Threads.atomic_add!(atomic, batch_size)
        while k ≤ len
            range = k : min(k + batch_size - 1, len)
            r = batch_mapreduce(r, range, f, op, src...)
            k = Threads.atomic_add!(atomic, batch_size)
        end
        Threads.lock(lock)
        mapreducer.r[] = op(mapreducer.r[], r)
        Threads.unlock(lock)
    end
    mapreducer.r[]
end

@static if VERSION < v"0.7-"
    # we assume that f.(src) and init are a subset of Abelian group with op
    function tmapreduce(f, op, init, src::AbstractArray...; batch_size=default_batch_size(length(src[1])))
        T = get_reduction_type(init, f, op, src...)
        _tmapreduce(T, init, batch_size, f, op, src...)
    end

    function tmapreduce(::Type{T}, f, op, init, src::AbstractArray...; batch_size=default_batch_size(length(src[1]))) where T
        _tmapreduce(T, init, batch_size, f, op, src...)
    end
else
    # we assume that f.(src) and init are a subset of Abelian group with op
    function tmapreduce(f, op, src::AbstractArray...; init, batch_size=default_batch_size(length(src[1])))
        T = get_reduction_type(init, f, op, src...)
        _tmapreduce(T, init, batch_size, f, op, src...)
    end

    function tmapreduce(::Type{T}, f, op, src::AbstractArray...; init, batch_size=default_batch_size(length(src[1]))) where T
        _tmapreduce(T, init, batch_size, f, op, src...)
    end
end

@inline function _tmapreduce(::Type{T}, init, batch_size, f, op, src...) where T
    lss = extrema(length.(src))
    lss[1] == lss[2] || throw(ArgumentError("src vectors must have the same length"))

    atomic = Threads.Atomic{Int}(1)
    lock = Threads.SpinLock()
    len = lss[1]
    mapreducer = MapReducer{T}(Base.RefValue{T}(init), atomic, lock, len)
    return mapreducer(batch_size, f, op, src...)
end

@static if VERSION < v"0.7-"
    const return_type = Core.Inference.return_type
else
    const return_type = Core.Compiler.return_type
end

@inline function get_reduction_type(init, f, op, src...)
    Tx = return_type(f, Tuple{eltype.(src)...})
    Trinit = return_type(op, Tuple{typeof(init), Tx})
    Tr = return_type(op, Tuple{Trinit, Tx})
    Tr === Union{} ? typeof(init) : Tr
end

@inline function batch_mapreduce(r, range, f, op, src...)
    @inbounds for i in range
        r = op(r, f(getindex.(src, i)...))
    end
    r
end

function getrange(n)
    tid = Threads.threadid()
    nt = Threads.nthreads()
    d , r = divrem(n, nt)
    from = (tid - 1) * d + min(r, tid - 1) + 1
    to = from + d - 1 + (tid ≤ r ? 1 : 0)
    from:to
end

end


text:
# ------------------------------------------------------------------------------

function generator_var_change(generator, var::P, numer::P, denom::P) where P <: MPolyElem
    return IterTools.imap(
        point -> begin
            result = copy(point)
            result[var] = eval_at_dict(numer, point) // eval_at_dict(denom, point)
            return result
    end, Iterators.filter(point -> eval_at_dict(denom, point) != 0, generator)
    )
end

# ------------------------------------------------------------------------------

function diff_poly(poly::P, derivation::Dict{P,P}) where P <: MPolyElem
    return sum(derivative(poly, x) * xd for (x, xd) in derivation)
end

# ------------------------------------------------------------------------------

function generate_io_equation_problem(ode::ODE{P}) where P <: MPolyElem{<: FieldElem}
    dim_x = length(ode.x_vars)

    # Creating a ring
    old_vars = map(string, gens(ode.poly_ring))
    var_names = vcat(
        old_vars,
        [var_to_str(x) * "_dot" for x in ode.x_vars],
        [var_to_str(y) * "_$i" for i in 0:dim_x for y in ode.y_vars],
        [var_to_str(u) * "_$i" for i in 0:dim_x for u in ode.u_vars],
        ["rand_proj_var"]
    )
    ring, ring_vars = Nemo.PolynomialRing(base_ring(ode.poly_ring), var_names)

    # Definiting a (partial) derivation on it
    derivation = Dict{P,P}()
    for x in ode.x_vars
        derivation[switch_ring(x, ring)] = str_to_var(var_to_str(x) * "_dot", ring)
    end
    for i in 0:(dim_x - 1)
        for y in ode.y_vars
            derivation[str_to_var(var_to_str(y) * "_$i", ring)] = str_to_var(var_to_str(y) * "_$(i + 1)", ring)
        end
    end
    for u in ode.u_vars
        for i in 0:(dim_x - 1)
            derivation[str_to_var(var_to_str(u) * "_$i", ring)] = str_to_var(var_to_str(u) * "_$(i + 1)", ring)
        end
    end

    # Generate equations
    old_us = [parent_ring_change(u, ring) for u in ode.u_vars]
    new_us = [str_to_var(var_to_str(u) * "_0", ring) for u in ode.u_vars]
    function to_new_ring(p)
        return evaluate(parent_ring_change(ode.poly_ring(p), ring), old_us, new_us)
    end
    x_equations = Dict{P,P}()
    for x in ode.x_vars
        x_lifted = parent_ring_change(x, ring)
        num, den = map(to_new_ring, unpack_fraction(ode.x_equations[x]))
        x_equations[x_lifted] = den * derivation[x_lifted] - num
    end
    y_equations = Dict{P,P}()
    for y in ode.y_vars
        y_lifted = str_to_var(var_to_str(y) * "_0", ring)
        g_num, g_den = map(to_new_ring, unpack_fraction(ode.y_equations[y]))
        y_equations[y_lifted] = g_den * y_lifted - g_num
    end

    generic_point_generator = ODEPointGenerator{P}(ode, ring)

    return (ring, derivation, x_equations, y_equations, generic_point_generator)
end

# ------------------------------------------------------------------------------

"""
    find_ioequations(ode, [var_change_policy=:default])

Finds the input-output equations of an ODE system
Input:
- `ode` - the ODE system
- `var_change_policy` - whether to perform automatic variable change, can be one of `:default`, `:yes`, `:no`

Output:
- a dictionary from "leaders" to the corresponding input-output equations
"""
function find_ioequations(
        ode::ODE{P};
        var_change_policy=:default
    ) where P <: MPolyElem{<: FieldElem}
    # Setting the var_change policy
    if (var_change_policy == :yes) || (var_change_policy == :default && length(ode.y_vars) < 3)
        auto_var_change = true
    elseif (var_change_policy == :no) || (var_change_policy == :default && length(ode.y_vars) >= 3)
        auto_var_change = false
    else
        @error "Unknown var_change policy $var_change_policy"
        return
    end

    # Initialization
    ring, derivation, x_equations, y_equations, point_generator = generate_io_equation_problem(ode)
    y_orders = Dict(y => 0 for y in keys(y_equations))

    while true
        var_degs = [(y, [Nemo.degree(eq, x) for x in keys(x_equations) if Nemo.degree(eq, x) > 0]) for (y, eq) in y_equations]
        filter!(d -> length(d[2]) > 0, var_degs)
        if isempty(var_degs)
            break
        end
        @debug "Current degrees of io-equations $var_degs"
        @debug "Orders: $y_orders"
        @debug "Sizes: $(Dict(y => length(eq) for (y, eq) in y_equations))"

        # choosing the output to prolong
        outputs_with_scores = [
            (
                min(d[2]...) * length(y_equations[d[1]]),
                min(d[2]...),
                -count(x -> x == min(d[2]...), d[2]) + length(y_equations[d[1]]) // 30,
                length(y_equations[d[1]]),
                d[1]
            ) for d in var_degs
        ]
        @debug "Scores: $outputs_with_scores"
        y_prolong = sort(outputs_with_scores)[1][end]
        y_orders[y_prolong] += 1
        @debug "Prolonging output $y_prolong"
        flush(stdout)

        # Calculate the Lie derivative of the io_relation
        @debug "Prolonging"
        flush(stdout)
        next_y_equation = diff_poly(y_equations[y_prolong], derivation)
        for (x, eq) in x_equations
            @debug "Eliminating the derivative of $x"
            flush(stdout)
            next_y_equation = eliminate_var(eq, next_y_equation, derivation[x], point_generator)
        end
        for (y, eq) in y_equations
            if y != y_prolong
                @debug "Eliminating the leader of the equation for $y"
                flush(stdout)
                # an ugly way of gettin the leader, to replace
                next_y_equation = eliminate_var(
                    next_y_equation, eq,
                    str_to_var(var_to_str(y)[1:end - 2] * "_$(y_orders[y])", ring),
                    point_generator
                )
            end
        end

        # Choose variable to eliminate
        var_degs_next = [(Nemo.degree(y_equations[y_prolong], x), Nemo.degree(next_y_equation, x), x) for x in keys(x_equations) if Nemo.degree(y_equations[y_prolong], x) > 0]
        our_choice = sort(var_degs_next)[1]
        var_elim_deg, var_elim = our_choice[1], our_choice[3]

        @debug "Elimination of $var_elim, $(length(x_equations)) left; $(Dates.now())"
        flush(stdout)

        # Possible variable change for Axy + Bx + p(y) (x = var_elim)
        if auto_var_change && (var_elim_deg == 1)
            Ay_plus_B = coeff(y_equations[y_prolong], [var_elim], [1])
            for x in setdiff(keys(x_equations), [var_elim])
                if Nemo.degree(Ay_plus_B, x) == 1
                    A, B = divrem(Ay_plus_B, x)
                    A, B = simplify_frac(A, B)
                    if isempty(filter!(v -> (v in keys(x_equations)), vars(A))) && (B != 0) # && (length(coeffs(A))==1)
                        # variable change x_i' --> x_i' - (B/A)', x_i --> x_i - B/A
                        @debug "\t Applying variable change: $(x) --> $(x) - ( $B )/( $A ); $(Dates.now())"
                        flush(stdout)
                        dB = diff_poly(B, derivation)
                        dA = diff_poly(A, derivation)
                        numer_d, denom_d = simplify_frac(A * dB - dA * B, A * A)
                        # !!!push x_i_dot first, otherwise there will be problem with plugin!!!
                        point_generator = generator_var_change(
                            point_generator,
                            derivation[x],
                            derivation[x] * denom_d + numer_d, denom_d
                        )
                        point_generator = generator_var_change(point_generator, x, A * x + B, A)

                        # Change current system
                        @debug "Change in the system"
                        flush(stdout)
                        x_equations[x] = make_substitution(x_equations[x], derivation[x], denom_d * derivation[x] - numer_d, denom_d)
                        for xx in keys(x_equations)
                            @debug "\t Change in the equation for $xx"
                            flush(stdout)
                            x_equations[xx] = make_substitution(x_equations[xx], x, A * x - B, A)
                        end
                        @debug "Change in the outputs"
                        flush(stdout)
                        for y in keys(y_equations)
                            @debug "\t Change in the output $y"
                            flush(stdout)
                            y_equations[y] = make_substitution(y_equations[y], x, A * x - B, A)
                        end
                        @debug "\t Change in the prolonged equation"
                        flush(stdout)
                        next_y_equation = make_substitution(next_y_equation, x, A * x - B, A)
                        # recalibrate system
                        @debug "Unmixing the derivatives"
                        flush(stdout)
                        for xx in setdiff(keys(x_equations), [x])
                            @debug "\t Unmixing $xx"
                            flush(stdout)
                            x_equations[x] = eliminate_var(x_equations[x], x_equations[xx], derivation[xx], point_generator)
                        end
                        @debug "Change of variables performed"
                        flush(stdout)
                        break
                    end
                end
            end
        end

        # Eliminate var_elim from the system
        delete!(x_equations, var_elim)
        @debug "Elimination in states"
        flush(stdout)
        for (x, eq) in x_equations
            @debug "\t Elimination in the equation for $x"
            flush(stdout)
            x_equations[x] = eliminate_var(eq, y_equations[y_prolong], var_elim, point_generator)
        end
        @debug "Elimination in y_equations"
        flush(stdout)
        for y in keys(y_equations)
            if y != y_prolong
                @debug "Elimination in the output $y"
                flush(stdout)
                y_equations[y] = eliminate_var(y_equations[y], y_equations[y_prolong], var_elim, point_generator)
            end
        end
        @debug "\t Elimination in the prolonged equation"
        flush(stdout)
        y_equations[y_prolong] = eliminate_var(y_equations[y_prolong], next_y_equation, var_elim, point_generator)
    end

    io_equations = Dict(str_to_var(var_to_str(y)[1:end - 2] * "_$(y_orders[y])", ring) => p for (y, p) in y_equations)

    @debug "Check whether the original projections are enough"
    if length(io_equations) == 1 || check_primality(io_equations)
        @debug "The projections generate an ideal with a single components of highest dimension, returning"
        return io_equations
    end

    sampling_range = 5
    while true
        @debug "There are several components of the highest dimension, trying to isolate one"
        extra_var = str_to_var("rand_proj_var", ring)
        extra_var_rhs = sum(rand(1:sampling_range) * v for v in keys(io_equations))
        @debug "Extra projections: $extra_var_rhs"
        point_generator = generator_var_change(point_generator, extra_var, extra_var_rhs, one(ring))
        extra_poly = extra_var - extra_var_rhs
        for (y, io_eq) in io_equations
            extra_poly = eliminate_var(extra_poly, io_eq, y, point_generator)
        end
        extra_poly = evaluate(extra_poly, [extra_var], [extra_var_rhs])
        if check_primality(io_equations, [extra_poly])
            @debug "Single component of highest dimension isolated, returning"
            io_equations[extra_var] = extra_poly
            break
        end
        sampling_range = 2 * sampling_range
    end
    return io_equations
end

ext: jl


text:
using LinearAlgebra
using FluxDiffUtils
using ForwardDiff
using SparseArrays
using StaticArrays
using Test

@testset "Flux diff tests" begin
    u = collect(LinRange(-1,1,10))
    U = (u,2*u)
    function flux1D(UL,UR)
        uL,vL = UL
        uR,vR = UR
        avg(x,y) = @. .5*(x+y)
        Fx = SVector{2}(avg(uL,uR),avg(vL,vR))
        return Fx
    end
    A = randn(10,10)
    A[8:10,8:10] .= 0 # for testing skip_index
    skip_index(i,j) = (i>8)&&(j>8)
    A = A-A' # make skew for exact formula testing
    ATr = (A->Matrix(transpose(A)))(A)

    rhs = hadamard_sum((A,),tuple∘flux1D,U,skip_index=skip_index)
    @test (A*U[1] + U[1].*sum(A,dims=2)) ≈ 2*rhs[1]
    @test (A*U[2] + U[2].*sum(A,dims=2)) ≈ 2*rhs[2]
end

# make 3-field solution
u = collect(LinRange(-1,1,10))
U = SVector{3}(u,2*u,3*u)

# define flux function and its derivative
function flux_function(UL,UR)
    uL,vL,wL = UL
    uR,vR,wR = UR
    avg(x,y) = @. .5*(x+y)
    Fx = avg(uL,uR),avg(vL,vR),avg(wL,wR)
    Fy = avg(uL,uR),avg(vL,vR),avg(wL,wR)
    return SVector{3}(Fx),SVector{3}(Fy)
end

dfx(UL,UR) = ForwardDiff.jacobian(UR->flux_function(UL,UR)[1], SVector{3}(UR))
dfy(UL,UR) = ForwardDiff.jacobian(UR->flux_function(UL,UR)[2], SVector{3}(UR))
df(UL,UR) = dfx(UL,UR),dfy(UL,UR)

A = randn(10,10)
A[8:10,8:10] .= 0 # for testing skip_index
skip_index(i,j) = (i>8)&(j>8)
A = ntuple(x->A,2)
A = (A->A-A').(A) # make skew for exact formula testing

@testset "Flux diff tests w/StaticArrays" begin

    rhs = hadamard_sum(A,flux_function,U)
    @test (A[1]*U[1] + U[1].*sum(A[1],dims=2)) ≈ rhs[1]
    @test (A[1]*U[2] + U[2].*sum(A[1],dims=2)) ≈ rhs[2]
    @test (A[1]*U[3] + U[3].*sum(A[1],dims=2)) ≈ rhs[3]

    rhs_sparse = hadamard_sum(sparse.(A),flux_function,U)
    @test all(rhs .≈ rhs_sparse)

    rhs_skipped = hadamard_sum(sparse.(A),flux_function,U;
                               skip_index=skip_index)
    @test all(rhs .≈ rhs_skipped)

    # test tuple args too
    v = collect(LinRange(-1,1,10))
    V = (v,2*v,3*v)
    # define flux function and its derivative
    function flux_function_tuple(UL,UR)
        uL,vL,wL = UL
        uR,vR,wR = UR
        avg(x,y) = @. .5*(x+y)
        Fx = avg(uL,uR),avg(vL,vR),avg(wL,wR)
        Fy = avg(uL,uR),avg(vL,vR),avg(wL,wR)
        return SVector{3}(Fx),SVector{3}(Fy)
    end
    rhs2 = hadamard_sum(A,flux_function_tuple,V)
    @test all(rhs .≈ rhs2)

end

@testset "Jacobian tests" begin
    tol = 5e2*eps()

    Jdense = hadamard_jacobian(A,df,U)
    J = hadamard_jacobian(sparse.(A),df,U)
    # no factor of .5 since Alist = (A,A) and exact jac is ∑_i (A_i - diag(sum(A_i,2)))
    @test J[1,1] ≈ first(A) - diagm(vec(sum(first(A),dims=1)))
    @test J[2,2] ≈ first(A) - diagm(vec(sum(first(A),dims=1)))
    @test J[3,3] ≈ first(A) - diagm(vec(sum(first(A),dims=1)))
    accum_norm = 0.0
    for i = 1:3, j = 1:3
        if i!=j
            accum_norm += norm(J[i,j])
        end
    end
    @test accum_norm < tol

    Jskip = hadamard_jacobian(A,df,U,skip_index=skip_index)
    test_dense_eq_sparse = Bool[]
    test_skip = Bool[]
    for i=1:3,j=1:3
        push!(test_dense_eq_sparse, J[i,j] ≈ sparse(Jdense[i,j]))
        push!(test_skip, J[i,j] ≈ sparse(Jskip[i,j]))
    end
    @test all(test_dense_eq_sparse)
    @test all(test_skip)

    # for banded diagonal matrix test
    g(u) = SMatrix{2,2}([u[1] 0
                           0 u[2]])
    G = banded_function_evals(g,U[1:2])
    @test G[1,1] ≈ spdiagm(0=>U[1])
    @test norm(G[1,2]) + norm(G[2,1]) < tol
    @test G[2,2] ≈ spdiagm(0=>U[2])

    # test blockcat
    B = [[i] for i=1:2,j=1:2]
    @test blockcat(size(B,2),B) ≈ [1 1;2 2]

    # test flattening of matrix
    @test blockcat(size(J,2),J) ≈ kron(I(3),first(A) - diagm(vec(sum(first(A),dims=1))))
end

@testset "Extra sparse/dense tests" begin
    tol = 5e2*eps()
    A,B = ntuple(x->sprandn(25,25,.1),2)
    A_list = (A,B)
    function flux2D(UL,UR)
        uL,vL = UL
        uR,vR = UR
        return (uL+uR),2*(vL+vR)
    end
    U = randn(25),randn(25)
    rhs_exact = (A*U[1] + U[1].*sum(A,dims=2), 2*B*U[2] + 2*U[2].*sum(B,dims=2))
    rhs1 = zero.(U)
    rhs2 = zero.(U)
    hadamard_sum_ATr!(rhs1,A_list,flux2D,U)
    hadamard_sum_ATr!(rhs2,Matrix.(A_list),flux2D,U) # dense version
    @test sum(norm.(rhs1 .- rhs2)) < tol
    tol = 5e2*eps()
    N = 10

    # scalar test
    A = sparse(randn(N,N))
    ATr_list = (Matrix(transpose(A)),)
    function flux1(UL,UR)
        return UL .+ UR
    end
    U = (randn(N),)
    rhs_exact = (A*U[1] + U[1].*sum(A,dims=2),)
    rhs1, rhs2 = zero.(U), zero.(U)
    hadamard_sum_ATr!(rhs1, ATr_list,flux1,U)
    hadamard_sum_ATr!(rhs2, sparse.(ATr_list),flux1,U)
    @test sum(norm.(rhs1 .- rhs_exact)) < tol
    @test sum(norm.(rhs_exact .- rhs2)) < tol
    @test sum(norm.(rhs1 .- rhs2)) < tol

    # 2-component 2D flux test
    A,B = randn(N,N),randn(N,N)
    ATr_list = map(A->Matrix(transpose(A)),(A,B))
    function flux2(UL,UR)
        uL,vL = UL
        uR,vR = UR
        return SVector{2}(uL+uR,0),SVector{2}(0,2*(vL+vR))
    end
    u,v = (randn(N),randn(N))
    U = (u,v)

    # computes sum(A.*F1 + B.*F2,dims=2)
    rhs_exact = (A*u + u.*sum(A,dims=2),2*B*v + 2*v.*sum(B,dims=2))
    rhs1, rhs2 = zero.(U), zero.(U)
    hadamard_sum_ATr!(rhs1,ATr_list,flux2,U)
    hadamard_sum_ATr!(rhs2,sparse(ATr_list[1]),(uL,uR)->flux2(uL,uR)[1],U)
    hadamard_sum_ATr!(rhs2,sparse(ATr_list[2]),(uL,uR)->flux2(uL,uR)[2],U)

    @test sum(norm.(rhs1 .- rhs_exact)) < tol
    @test sum(norm.(rhs_exact .- rhs2)) < tol
    @test sum(norm.(rhs1 .- rhs2)) < tol
end


text:
# Array Wrapper around a persistent memory region.
"""
    PersistentArray{T,N} <: AbstractArray{T,N}

`N`-dimensional dense array with element of type `T` - backed by persistent memory NVDIMMs
instead of DRAM.

```julia
PersistentArray(A::AbstractArray{T,N}) -> PersistentArray{T,N}
```

Construct a `PersistentArray` copy of `A`.

```julia
PersistentArray{T}(undef, dims)
PersistentArray{T,N}(undef, dims)
```

Construct an unitialized `N`-dimensional `PersistentArray` containing elements of type `T`.
"""
mutable struct PersistentArray{T,N} <: DenseArray{T,N}
    ptr::Ptr{T}
    size::NTuple{N, Int}
    file::String

    function PersistentArray{T,N}(::UndefInitializer, d::NTuple{N,Int}) where {T,N}
        # Make sure `T` is an `isbitstype` so that the backing data is actually stored in
        # the persistent memory region.
        if !isbitstype(T)
            err = ArgumentError("""
                Expected elements of PersistentArrays to be `isbitstype`.
                Got an eltype of: $(T)
                """
            )
            throw(err)
        end

        # Determine the number of bytes needed for this array
        bytes = prod(d) * sizeof(T)
        (ptr_nothing, file) = alloc(bytes)
        ptr = Base.unsafe_convert(Ptr{T}, ptr_nothing)

        array = new{T,N}(ptr, d, file)
        finalizer(destroy, array)
        return array
    end
end

PersistentArray{T}(::UndefInitializer, d::NTuple{N,Int}) where {T,N} =
    PersistentArray{T,N}(undef, d)

PersistentArray{T}(::UndefInitializer, a::Integer, b...) where {T} =
    PersistentArray{T}(undef, convert.(Int, (a, b...)))

PersistentArray{T}(::UndefInitializer, d::Integer) where {T} =
    PersistentArray{T,1}(undef, (Int(d),))

function PersistentArray(A::AbstractArray{T,N}) where {T,N}
    P = PersistentArray{T,N}(undef, size(A))
    P .= A
    return P
end

# Basic accessor.
Base.pointer(P::PersistentArray) = P.ptr

# needed to support views? Why do the fallbacks not work??
Base.elsize(::PersistentArray{T}) where {T} = sizeof(T)

# Clean up the backing file as well.
destroy(P::PersistentArray) = free(Base.unsafe_convert(Ptr{Nothing}, P.ptr), true)

# Implement array interface.
Base.size(P::PersistentArray) = P.size
Base.sizeof(P::PersistentArray{T}) where {T} = length(P) * sizeof(T)
Base.@propagate_inbounds function Base.getindex(P::PersistentArray, i::Integer)
    @boundscheck checkbounds(P, i)
    return unsafe_load(pointer(P), i)
end
Base.@propagate_inbounds function Base.setindex!(P::PersistentArray, v, i::Integer)
    @boundscheck checkbounds(P, i)
    return unsafe_store!(pointer(P), v, i)
end

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

# Allow conversion to pointer to dispatch to GEMM.
Base.unsafe_convert(::Type{Ptr{T}}, P::PersistentArray{T}) where {T} = pointer(P)


text:
function runLocalCommTests(origComm::Comm{GID, PID, LID}) where{GID, PID, LID}
    localComm = LocalComm(origComm)

    @test_throws InvalidStateError barrier(localComm)
    @test_throws InvalidStateError broadcastAll(localComm, [1, 2, 3], 1)
    @test_throws InvalidStateError broadcastAll(localComm, 18, 1)
    @test_throws InvalidStateError createDistributor(localComm)

    @test_throws InvalidStateError gatherAll(localComm, [4, 5, 6, 7])
    @test_throws InvalidStateError gatherAll(localComm, 8)

    @test_throws InvalidStateError sumAll(localComm, [4, 5, 6, 7])
    @test_throws InvalidStateError sumAll(localComm, 8)

    @test_throws InvalidStateError minAll(localComm, [4, 5, 6, 7])
    @test_throws InvalidStateError minAll(localComm, 8)

    @test_throws InvalidStateError maxAll(localComm, [4, 5, 6, 7])
    @test_throws InvalidStateError maxAll(localComm, 8)

    @test_throws InvalidStateError scanSum(localComm, [4, 5, 6, 7])
    @test_throws InvalidStateError scanSum(localComm, 8)

    @test myPid(origComm) == myPid(localComm)
    @test numProc(origComm) == numProc(localComm)
end


text:
#=

# 4: Lid-driven cavity

=#

module cavity_flow

using Printf
using SmoothedParticles

#=
Declare const parameters (all dims in SI)
=#

##geometrical/physical parameters
const Re = 100.0                #Reynolds number
const llid = 0.2                #length of the lid
const mu = 8.4e-4               #viscosity of water
const rho0 = 1.0e+3             #density of water
const vlid = mu*Re/(rho0*llid)     #flow speed of the lid
const dr = llid/70 		        #average particle distance
const h = 2.2*dr		        #size of kernel support
const m = rho0*dr^2             #particle mass
const c = 10*vlid			#numerical speed of sound
const wwall = h

##temporal parameters
const dt = 0.2*h/c              #numerical time-step
const dt_frame = 20.           #how often save data
const t_end = 2000.            #end of simulation

##particle types
const FLUID = 0.
const WALL = 1.
const LID = 2.

#=
Declare variables to be stored in a Particle
=#

mutable struct Particle <: AbstractParticle
	x::RealVector	#position
	v::RealVector #velocity
	a::RealVector #acceleratation
	rho::Float64 #density
	Drho::Float64 #rate of density
	P::Float64 #pressure
	type::Float64 #particle type
	Particle(x::RealVector, type::Float64) = begin
		speed = (type == LID ? vlid : 0.)
		return new(x, speed*VECX, VEC0, rho0, 0., 0., type)
	end
end

#=
Define geometry and create particles
=#

function make_system()
	grid = Grid(dr, :hexagonal)
	box = Rectangle(0., 0., llid, llid)
	wall = BoundaryLayer(box, grid, wwall)
	sys = ParticleSystem(Particle, box + wall, h)

	lid   = Specification(wall, x -> x[2] > llid)
	wall = Specification(wall, x -> x[2] <= llid)

	generate_particles!(sys, grid, box, x -> Particle(x, FLUID))
	generate_particles!(sys, grid, lid, x -> Particle(x, LID))
	generate_particles!(sys, grid, wall, x -> Particle(x, WALL))
	return sys
end

#=
Define interactions between particles
=#

#Define interactions between particles

@inbounds function balance_of_mass!(p::Particle, q::Particle, r::Float64)
	if p.type == FLUID
		ker = m*rDwendland2(h,r)
		p.Drho += ker*(dot(p.x-q.x, p.v-q.v))
	end
end

function find_pressure!(p::Particle)
	p.rho += p.Drho*dt
	p.Drho = 0.0
	p.P = c^2*(p.rho-rho0)
end

@inbounds function internal_force!(p::Particle, q::Particle, r::Float64)
	ker = m*rDwendland2(h,r)
	p.a += -ker*(p.P/rho0^2 + q.P/rho0^2)*(p.x - q.x)
	p.a += +2*ker*(mu/rho0^2)*(p.v - q.v)
end

function move!(p::Particle)
	p.a = VEC0
	if p.type == FLUID
		p.x += dt*p.v
	end
end

function accelerate!(p::Particle)
	if p.type == FLUID
		p.v += 0.5*dt*p.a
	end
end

#=
Time iteration
=#
function main()
	sys = make_system()
	out = new_pvd_file("results/cavity_flow")
	@time for k = 0 : Int64(round(t_end/dt))
		apply!(sys, move!)
		create_cell_list!(sys)
		apply!(sys, balance_of_mass!)
		apply!(sys, find_pressure!)
		apply!(sys, internal_force!)
		apply!(sys, accelerate!)
		if (k % Int64(round(dt_frame/dt)) == 0) #save the frame
			t = k*dt
			@show t
			save_frame!(out, sys, :P, :v, :type)
		end
		apply!(sys, accelerate!)
	end
	save_pvd_file(out)
end

end


text:
using MUC16Analysis
using BioSequences
using DataFrames
using Test
using CSV

@testset "MUC16Analysis.jl" begin
    @test 10 == 10 #just to make sure test is working

    @test letter_list(5) == ["A", "B", "C", "D", "E"]

    SUBMAT_ALI = sub_mat_ali() #create the ali et al substitution matrix
    @test SUBMAT_ALI[AA_A,AA_R] === 90 #check random values in the substitution matrix
    @test SUBMAT_ALI[AA_G,AA_D] === 65
    @test SUBMAT_ALI[AA_P,AA_P] === 0

    #run upgma on example dataframe
    example_dmdf = CSV.read(joinpath(dirname(@__FILE__), "example.csv"))
    add_missing(example_dmdf)
    example_newick = upgma(example_dmdf, output = "example.tree", verbose = false)
    println(example_newick)
    @test example_newick == "((((A:4.0,D:4.0):4.25,((B:0.5,F:0.5):5.75,G:6.25):2.0):6.25,C:14.5):2.5,E:17.0);"

    #run upgma on example distance matrix by different method
    example_dmdf2 = create_dm_df(joinpath(dirname(@__FILE__), "example.csv"))
    example_newick2 = upgma(example_dmdf2, output = "example2.tree", verbose = true)
    @test example_newick2 == "((((A:4.0,D:4.0):4.25,((B:0.5,F:0.5):5.75,G:6.25):2.0):6.25,C:14.5):2.5,E:17.0);"

    #create distance matrix from sample fasta. Fasta sequences for test were randomly generated.
    example_dm_fromfasta = dist_mat_gen(joinpath(dirname(@__FILE__), "example.fasta"), "example_dm_fromfasta.csv")
    @test nrow(example_dm_fromfasta) == 10
    add_missing(example_dm_fromfasta)
    example_newick3 = upgma(example_dm_fromfasta, output = "example3.tree")
    @test example_newick3 == "((((1:2068.0,6:2068.0):138.70833333333348,(((2:1948.5,4:1948.5):155.0,10:2103.5):66.77777777777783,(3:2102.5,(5:2001.0,8:2001.0):101.5):67.77777777777783):36.43055555555566):12.229166666666515,7:2218.9375):63.17361111111131,9:2282.1111111111113);"
end


text:
function to_rna(dna::AbstractString)
    if length(split(dna, r"A|C|G|T", keep=false)) >= 1
        return throw(ErrorException("Error"))
    else
        tmp = ""
        for i in dna
            if i == 'G'
                tmp = string(tmp, "C")
            elseif i == 'C'
                tmp = string(tmp, "G")
            elseif i == 'T'
                tmp = string(tmp, "A")
            elseif i == 'A'
                tmp = string(tmp, "U")
            end
        end
        return tmp
    end
end

to_rna("A")

max_stars_repo_name: Devlious/LearningJulia

max_issues_repo_name: Devlious/LearningJulia

max_forks_repo_name: Devlious/LearningJulia


text:
module classification

#import Plots
import Statistics
using Requires
using ..utils

include("confusion_matrix.jl"); export confusion_params, confusion_matrix, class_confusion
include("metrics.jl"); export classification_report, condition_positive, condition_negative, predicted_positive,predicted_negative, correctly_classified, incorrectly_classified, sensitivity_score, recall_score, specificity_score, precision_score, positive_predictive_value, accuracy_score, balanced_accuracy_score, negative_predictive_value, false_negative_rate, false_positive_rate, false_discovery_rate, false_omission_rate, f1_score, prevalence_threshold, threat_score, matthews_correlation_coeff, fowlkes_mallows_index, informedness, markedness, cohen_kappa_score, hamming_loss, jaccard_score, confusion_params
include("roc_auc.jl"); export auc, RocCurve, PrecisionRecallCurve, ratio_at_threshold, fpr_at_threshold, tpr_at_threshold, auc_at_threshold, precision_at_threshold, recall_at_threshold

function __init__()
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" @eval begin
        include("visualization.jl")
        export visualize
    end
end

end


text:
function timer(f::Function, n::Integer)
    # Call once to force JIT compilation
    f()

    times = Array(Float64, n)
    for itr in 1:n
        times[itr] = @elapsed f()
    end

    return times
end

function benchmark(f::Function,
                   category::AbstractString,
                   name::AbstractString,
                   n::Integer)
    times = timer(f, n)

    df = DataFrame()
    df[:Category] = category
    df[:Benchmark] = name
    df[:Iterations] = n
    # TODO: Remove
    df[:TotalWall] = sum(times)
    df[:AverageWall] = mean(times)
    df[:MaxWall] = maximum(times)
    df[:MinWall] = minimum(times)
    df[:Timestamp] = Libc.strftime("%Y-%m-%d %H:%M:%S", round(Int, time()))
    df[:JuliaHash] = Base.GIT_VERSION_INFO.commit
    if isdir(".git")
        df[:CodeHash] = readchomp(`git rev-parse HEAD`)
    else
        df[:CodeHash] = NA
    end
    df[:OS] = string(OS_NAME)
    df[:CPUCores] = CPU_CORES

    return df
end

benchmark(f::Function, name::AbstractString, n::Integer) = benchmark(f, name, name, n)

function benchmarks(marks::Vector)
    df = DataFrame()

    for mark in marks
        df = vcat(df, benchmark(mark[1], mark[2], mark[3], mark[4]))
    end

    return df
end


text:
"""
    struct NystromMesh{N,T} <: AbstractMesh{N,T}

A mesh data structure for solving boundary integral equation using Nyström
methods.

A `NystromMesh` can be constructed from a `mesh::AbstractMesh` and a dictionary
`etype2qrule` mapping element types in `mesh` (the keys) to an appropriate
quadrature rule for integration over elements of that type (the value).

The degrees of freedom of a `NystromMesh` are associated to its quadrature nodes
`qnodes::Vector{SVector{N,T}}`.
"""
Base.@kwdef struct NystromMesh{N,T} <: AbstractMesh{N,T}
    elements::OrderedDict{DataType,Any} = OrderedDict{DataType,Any}()
    etype2qrule::OrderedDict{DataType,AbstractQuadratureRule}= OrderedDict{DataType,AbstractQuadratureRule}()
    # quadrature info
    qnodes::Vector{SVector{N,T}} = Vector{SVector{N,T}}()
    qweights::Vector{T} = Vector{T}()
    qnormals::Vector{SVector{N,T}} = Vector{SVector{N,T}}()
    # mappings from elements --> dofs and entities --> elements
    elt2dof::OrderedDict{DataType,Matrix{Int}} = OrderedDict{DataType,Matrix{Int}}()
    ent2elt::OrderedDict{AbstractEntity,OrderedDict{DataType,Vector{Int}}} = OrderedDict{AbstractEntity,OrderedDict{DataType,Vector{Int}}}()
end

function Base.append!(msh1::NystromMesh,msh2::NystromMesh)
    Ω = domain(msh2)
    @assert ambient_dimension(msh1) == ambient_dimension(msh2)
    @assert eltype(msh1) == eltype(msh2)
    # extract relevant quadrature rules
    etype2dof = OrderedDict{DataType,Vector{Int}}()
    for (E,q) in etypes(msh2,Ω)
        if haskey(msh1.elements,E)
            etype2dof[E]        = msh1.elt2dof[E] |> vec
            @assert msh1.etype2qrule[E] == q
        else
            etype2dof[E]        = Int[]
            msh1.elements[E]    = E[]
            msh1.etype2qrule[E] = q
        end
    end
    # loop over entities and extract the information
    for ent in entities(msh2)
        haskey(msh1.ent2elt,ent) && continue # entity already meshed
        msh1.ent2elt[ent] = OrderedDict{DataType,Vector{Int}}()
        for (E,v) in ent2elt(msh2,ent)
            # add dofs information
            msh2_dofs  = msh2.elt2dof[E][:,v]
            idx_start = length(msh1.qnodes) + 1
            append!(msh1.qnodes,msh2.qnodes[msh2_dofs])
            append!(msh1.qweights,msh2.qweights[msh2_dofs])
            if !isempty(qnormals(msh2)) # otherwise assume it did not have normal stored
                append!(msh1.qnormals,msh2.qnormals[msh2_dofs])
            end
            idx_end = length(msh1.qnodes) # end index of dofs added
            append!(etype2dof[E],collect(idx_start:idx_end))
            # append elements
            idx_start = length(msh1.elements[E]) + 1
            append!(msh1.elements[E],msh2.elements[E][v])
            idx_end   = length(msh1.elements[E])
            # and mapping ent -> elt
            msh1.ent2elt[ent][E] = collect(idx_start:idx_end)
        end
    end
    # finally do some rehshaping
    for (E,v) in etype2dof
        dof_per_el = size(msh2.elt2dof[E],1)
        msh1.elt2dof[E] = reshape(v,dof_per_el,:)
    end
    return msh1
end

# getters
qnodes(m::NystromMesh)   = m.qnodes
qweights(m::NystromMesh) = m.qweights
qnormals(m::NystromMesh) = m.qnormals
elements(m::NystromMesh) = m.elements

# various mappings
elt2dof(m::NystromMesh)  = m.elt2dof
elt2dof(m::NystromMesh,E::DataType) = m.elt2dof[E]
ent2elt(m::NystromMesh) = m.ent2elt
ent2elt(m::NystromMesh,ent::AbstractEntity) = m.ent2elt[ent]
etype2qrule(m)   = m.etype2qrule
etype2qrule(m,E) = etype2qrule(m)[E]

Geometry.domain(mesh::NystromMesh)   = Domain(entities(mesh))
Geometry.entities(mesh::NystromMesh) = collect(keys(mesh.ent2elt))

function dom2elt(mesh::NystromMesh,domain::Domain)
    dict = OrderedDict{DataType,Vector{Int}}()
    for ent in entities(domain)
        other = ent2elt(mesh,ent)
        mergewith!(append!,dict,other)
    end
    return dict
end

function dom2dof(mesh::NystromMesh,domain::Domain)
    idxs = Int[]
    for ent in entities(domain)
        dict = ent2elt(mesh,ent)
        for (E,tags) in dict
            append!(idxs,view(elt2dof(mesh,E),:,tags))
        end
    end
    return idxs
end
dom2dof(mesh,ent::AbstractEntity) = dom2dof(mesh,Domain(ent))
ent2dof(mesh,ent::AbstractEntity) = dom2dof(mesh,ent)

Base.length(m::NystromMesh) = length(qnodes(m))

function etypes(m::NystromMesh,Ω::Domain=domain(m))
    ee = Tuple{DataType,AbstractQuadratureRule}[]
    for ent in entities(Ω)
        dict = ent2elt(m,ent)
        EE    = keys(dict)
        for E in EE
            q = etype2qrule(m,E)
            push!(ee,(E,q))
        end
    end
    unique!(ee)
    return ee
end

function NystromMesh(mesh::GenericMesh,Ω::Domain,e2qrule,compute_normal::Bool=true)
    N,T        = ambient_dimension(mesh), eltype(mesh)
    # initialize empty fields
    nys_mesh = NystromMesh{N,T}(;etype2qrule=e2qrule)
    # loop over entities
    for ent in entities(Ω)
        dict    = OrderedDict{DataType,Vector{Int}}()
        submesh = view(mesh,ent)
        for E in etypes(submesh)
            haskey(e2qrule,E) || error("no quadrature rule found for element of type $E")
            qrule = e2qrule[E]
            if haskey(nys_mesh.elements,E)
                istart = length(nys_mesh.elements[E])+1
            else
                istart = 1
            end
            _build_nystrom_mesh!(nys_mesh,submesh,E,qrule,compute_normal)
            iend   = length(nys_mesh.elements[E])
            dict[E] = collect(istart:iend)
        end
        nys_mesh.ent2elt[ent] = dict # new entry
    end
    return nys_mesh
end

@noinline function _build_nystrom_mesh!(nys_mesh,submesh,E,qrule,compute_normal)
    els::Vector{E}         = get!(nys_mesh.elements,E,Vector{E}())
    el2qnodes   = get!(nys_mesh.elt2dof,E,Matrix{Int}(undef,0,0))
    el2qnodes   = vec(el2qnodes)
    x̂,ŵ         = qrule() #nodes and weights on reference element
    for el in ElementIterator(submesh,E)
        # compute quadrature on element
        x = map(el,x̂)
        w = map(zip(x̂,ŵ)) do (x̂,ŵ)
            ŵ*measure(el,x̂)
        end
        # add element
        push!(els,el)
        # append quadrature information
        istart = length(qnodes(nys_mesh)) + 1
        append!(qweights(nys_mesh),w)
        append!(qnodes(nys_mesh),x)
        if compute_normal
            ν = map(x->normal(el,x),x̂)
            append!(qnormals(nys_mesh),ν)
        end
        iend   = length(qnodes(nys_mesh))
        # add information to recover qnodes from index of el
        append!(el2qnodes,collect(istart:iend))
    end
    el2qnodes = reshape(el2qnodes,length(x̂),:)
    nys_mesh.elt2dof[E] = el2qnodes
    return nys_mesh
end

NystromMesh(submesh::SubMesh,args...) = NystromMesh(submesh.mesh,submesh.domain,args...)

function NystromMesh(mesh::SubMesh;order,compute_normal::Bool=true)
    etype2qrule = Mesh._qrule_for_mesh(mesh,order)
    NystromMesh(mesh,etype2qrule,compute_normal)
end

function NystromMesh(mesh::GenericMesh,Ω::Domain;quad_rule,compute_normal::Bool=true)
    etype2qrule = OrderedDict(E=>quad_rule for E in etypes(view(mesh,Ω)))
    NystromMesh(mesh,Ω,etype2qrule,compute_normal)
end

# construct a mesh as a restriction to a domain
function Base.getindex(msh2::NystromMesh,Ω::Domain)
    N,T      = ambient_dimension(msh2), eltype(msh2)
    # initialize empty fields
    msh1 = NystromMesh{N,T}()
    # extract relevant quadrature rules
    etype2dof = OrderedDict{DataType,Vector{Int}}()
    for (E,q) in etypes(msh2,Ω)
        msh1.etype2qrule[E] = q
        etype2dof[E] = Int[]
        msh1.elements[E] = E[]
    end
    # loop over entities and extract the information
    for ent in Ω
        @assert ent in domain(msh2)
        msh1.ent2elt[ent] = OrderedDict{DataType,Vector{Int}}()
        for (E,v) in msh2.ent2elt[ent]
            # add dofs information
            msh2_dofs  = msh2.elt2dof[E][:,v]
            idx_start = length(msh1.qnodes) + 1
            append!(msh1.qnodes,msh2.qnodes[msh2_dofs])
            append!(msh1.qweights,msh2.qweights[msh2_dofs])
            if !isempty(qnormals(msh2))
                append!(msh1.qnormals,msh2.qnormals[msh2_dofs])
            end
            idx_end = length(msh1.qnodes)
            append!(etype2dof[E],collect(idx_start:idx_end))
            # append elements
            idx_start = length(msh1.elements[E]) + 1
            append!(msh1.elements[E],msh2.elements[E][v])
            idx_end   = length(msh1.elements[E])
            # and mapping ent -> elt
            msh1.ent2elt[ent][E] = collect(idx_start:idx_end)
        end
    end
    # finally do some rehshaping
    for (E,v) in etype2dof
        dof_per_el = size(msh2.elt2dof[E],1)
        msh1.elt2dof[E] = reshape(v,dof_per_el,:)
    end
    return msh1
end
Base.getindex(m::NystromMesh,ent::AbstractEntity) = m[Domain(ent)]

function _elt2dof(m::NystromMesh,Ω::Domain)
    out = OrderedDict{DataType,Matrix{Int}}()
    for ent in Ω
        dict = ent2elt(m,ent)
        for (E,v) in dict
            dofs = get!(out,E,Matrix{Int}(undef,0,0))
            # make into vec to append
            dofs = vec(dofs)
            dof_per_el = size(elt2dof(m,E),1)
            append!(dofs,elt2dof(m,E)[:,v])
            # convert back to matrix shape
            dofs = reshape(dofs,dof_per_el,:)
            out[E] = dofs
        end
    end
    return out
end

"""
    near_interaction_list(pts,Y;atol)

Given a target points `pts` and a `Y::NystromMesh`, return a dictionary with keys
given by element types of the source mesh `Y`. To each key, which represents an
element type, we associate a vector whose `i`-th entry encodes
information about the points in `X` which are close to a given element of
the key type.
"""
function near_interaction_list(pts,Y::NystromMesh;atol)
    dict = OrderedDict{DataType,Vector{Vector{Tuple{Int,Int}}}}()
    for (E,Q) in etypes(Y)
        ilist = _etype_near_interaction_list(pts,Y,E,atol)
        push!(dict,E=>ilist)
    end
    return dict
end

function _etype_near_interaction_list(pts,Y,E,atol)
    ilist = Vector{Vector{Tuple{Int,Int}}}()
    e2n   = elt2dof(Y,E)
    npts,nel = size(e2n)
    for n in 1:nel
        ynodes = qnodes(Y)[e2n[:,n]]
        inear  = _near_interaction_list(pts,ynodes,atol)
        push!(ilist,inear)
    end
    ilist
end

function _near_interaction_list(pts,ynodes,atol)
    ilist    = Tuple{Int,Int}[]
    for (i,x) in enumerate(pts)
        d,j = findmin([norm(x-y) for y in ynodes])
        if d ≤ atol
            push!(ilist,(i,j))
        end
    end
    return ilist
end

elements(m::NystromMesh,E::DataType) = m.elements[E]

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
# Creation macros #
# ========================== #
"""
    @scenario(def)

Define a scenario type compatible with StochasticPrograms using the syntax
```julia
@scenario name = begin
    ...structdef...

    [@zero begin
        ...
        return zero(scenario)
    end]

    [@expectation begin
        ...
        return expected(scenarios)
     end]
end
```
The generated type is referenced through `name` and a default constructor is always generated. This constructor accepts the keyword `probability` to set the probability of the scenario occuring. Otherwise, any internal variables and specialized constructors are defined in the @scenario block as they would be in any Julia struct.

If possible, a `zero` method and an `expected` method will be generated for the defined type. Otherwise, or if the default implementation is not desired, these can be user provided through [`@zero`](@ref) and [`@expectation`](@ref).

The defined scenario type will be available on all Julia processes.

## Examples

The following defines a simple scenario ``ξ`` with a single value.

```jldoctest
@scenario ExampleScenario = begin
    ξ::Float64
end

ExampleScenario(1.0, probability = 0.5)

# output

ExampleScenario with probability 0.5
  ξ: 1.0

```

See also: [`@zero`](@ref), [`@expectation`](@ref), [`@sampler`](@ref)
"""
macro scenario(arg)
    @capture(arg, scenarioname_Symbol = scenariodef_) || error("Invalid syntax. Expected: scenarioname = begin scenariodef end")
    vars = Vector{Symbol}()
    vartypes = Vector{Union{Expr,Symbol}}()
    vardefs = Vector{Expr}()
    zerodefs = Vector{Expr}()
    expectdefs = Vector{Expr}()
    def = postwalk(prettify(scenariodef)) do x
        if @capture(x, var_Symbol::vartype_)
            push!(vars, var)
            push!(vartypes, vartype)
            push!(vardefs, x)
            return x
        end
        if @capture(x, @zero zerodef_)
            push!(zerodefs, zerodef)
            return @q begin end
        end
        if @capture(x, @expectation expectdef_)
            push!(expectdefs, expectdef)
            return @q begin end
        end
        return x
    end
    # Handle zero definition
    if length(zerodefs) > 1
        error("Only provide one zero implementation")
    end
    provided_zerodef, zerodef = if length(zerodefs) == 1
        zerodef = zerodefs[1]
        true, zerodef
    else
        false, Expr(:block)
    end
    # Handle expectation definition
    if length(expectdefs) > 1
        error("Only provide one expectation implementation")
    end
    provided_expectdef, expectdef = if length(expectdefs) == 1
        expectdef = expectdefs[1]
        expectdef = postwalk(prettify(expectdef)) do x
            @capture(x, constructor_Symbol(args__)) && constructor == scenarioname && return :(StochasticPrograms.ExpectedScenario($x))
            return x
        end
        true, expectdef
    else
        false, Expr(:block)
    end
    # Prepare automatic expectation definition
    combine = Expr(:call)
    push!(combine.args, scenarioname)
    push!(combine.args, :($(Expr(:parameters, :(probability = 1.0)))))
    for var in vars
        push!(combine.args, :(probability(s1)*s1.$var + probability(s2)*s2.$var))
    end
    # Prepare automatic scenariotext definition
    textdef = Expr(:block)
    for var in vars
        key = Meta.quot(var)
        push!(textdef.args, :(print(io, "\n  $($key): $(scenario.$var)")))
    end
    # Define scenario type
    code = @q begin
        if StochasticPrograms.supports_expected([$(vartypes...)], $provided_expectdef) || $provided_expectdef
            if StochasticPrograms.supports_zero([$(vartypes...)], $provided_zerodef) || $provided_zerodef
                struct $scenarioname <: AbstractScenario
                    probability::Probability
                    $def
                    function (::Type{$scenarioname})($(vardefs...); probability::AbstractFloat = 1.0)
                        return new(Probability(probability), $(vars...))
                    end
                end
                if $provided_zerodef
                    function Base.zero(::Type{$scenarioname})
                        $zerodef
                    end
                else
                    function Base.zero(::Type{$scenarioname})
                        return $scenarioname(zero.([$(vartypes...)])...; probability = 1.0)
                    end
                end
                if $provided_expectdef
                    function StochasticPrograms.expected(scenarios::Vector{$scenarioname})
                        isempty(scenarios) && return StochasticPrograms.ExpectedScenario(zero($scenarioname))
                        $expectdef
                    end
                else
                    function StochasticPrograms.expected(scenarios::Vector{$scenarioname})
                        isempty(scenarios) && return StochasticPrograms.ExpectedScenario(zero($scenarioname))
                        return StochasticPrograms.ExpectedScenario(reduce(scenarios) do s1, s2
                                                                   $combine
                                                                   end)
                    end
                end
                function StochasticPrograms.scenariotext(io::IO, scenario::$scenarioname)
                    $textdef
                    return io
                end
            else
                @warn "The scenario type $(string($(Meta.quot(scenarioname)))) was not defined. A user-provided implementation \n\n    function zero(::Type{{$(string($(Meta.quot(scenarioname))))})\n        ...\n    end\n\nis required."
            end
        else
            @warn "The scenario type $(string($(Meta.quot(scenarioname)))) was not defined. A user-provided implementation \n\n    function expected(scenarios::Vector{$(string($(Meta.quot(scenarioname))))})\n        ...\n    end\n\nis required."
        end
    end
    code = prettify(code)
    code = @q begin
        if (@__MODULE__) == $(esc(:Main))
            @everywhere begin
                $code
            end
        else
            $(esc(code))
        end
    end
    return prettify(code)
end
"""
    @zero(def)

Define the additive zero scenario inside a @scenario block using the syntax:
```julia
@zero begin
    ...
    return zero_scenario
end
```

## Examples

The following defines a zero scenario for the example scenario defined in [`@scenario`](@ref)

```julia
@zero begin
    return ExampleScenario(0.0)
end
```

See also [`@scenario`](@ref)
"""
macro zero(def) @warn "@zero should be used inside a @scenario block." end
"""
    @expectation(def)

Define how to form the expected scenario inside a @scenario block. The scenario collection is accessed through the reserved keyword `scenarios`.

```julia
@zero begin
    ...
    return zero_scenario
end
```

## Examples

The following defines expectation for the example scenario defined in [`@scenario`](@ref)

```julia
@expectation begin
    return ExampleScenario(sum([probability(s)*s.ξ for s in scenarios]))
end
```

See also [`@scenario`](@ref)
"""
macro expectation(def) @warn "@expectation should be used inside a @scenario block." end
"""
    @sampler(def)

Define a sampler type for some `scenariotype` compatible with StochasticPrograms using the syntax
```julia
@sampler samplername = begin
    ...internals...

    @sample scenariotype begin
        ...
        return scenario
    end
end
```
Any internal state required by the sampler, as well as any specialized constructor, are defined in the @sampler block as they would be in any Julia struct. Define the sample operation inside the [`@sample`](@ref) block and specify the `scenariotype` that the sampler returns. The defined sampler will be available on all Julia processes.

## Examples

The following defines a simple dummy sampler, with some internal weight value, for the scenario defined in [`@scenario`](@ref), and samples one scenario.

```jldoctest; setup = :(@scenario ExampleScenario = ξ::Float64), filter = r".*"
@sampler ExampleSampler = begin
    w::Float64

    ExampleSampler(w::AbstractFloat) = new(w)

    @sample ExampleScenario begin
        @parameters w
        return ExampleScenario(w*randn(), probability = rand())
    end
end
s = ExampleSampler(2.0)
s()

# output

ExampleScenario with probability 0.29
  ξ: 1.48

```

See also: [`@sample`](@ref), [`@scenario`](@ref)
"""
macro sampler(arg)
    @capture(prettify(arg), samplername_Symbol = samplerdef_) || error("Invalid syntax. Expected: sampler = begin samplerdef end")
    scenarioname = :undefined
    sampledefs = Vector{Expr}()
    def = postwalk(prettify(samplerdef)) do x
        if @capture(x, @parameters args__)
            code = Expr(:block)
            for param in args
                push!(code.args, :($param = sampler.$param))
            end
            return code
        elseif @capture(x, @sample sname_ sampledef_)
            scenarioname = sname
            push!(sampledefs, sampledef)
            return @q begin end
        else
            return x
        end
    end
    # Handle sample definition
    length(sampledefs) == 1 || error("Must provide exactly one @sample operation.")
    sampledef = sampledefs[1]
    # Define sampler type
    code = @q begin
        struct $samplername <: AbstractSampler{$scenarioname}
            $def
        end
        function (sampler::$samplername)()
            $sampledef
        end
    end
    code = prettify(code)
    code = @q begin
        if (@__MODULE__) == $(esc(:Main))
            @everywhere begin
                $code
            end
        else
            $(esc(code))
        end
    end
    return prettify(code)
end
"""
    @sample(def)

Define the sample operaton inside a @sampler block, using the syntax
```julia
@sample begin
    ...
    return sampled_scenario
end
```
The sampler object is referenced through the reserved keyword `sampler`, from which any internals can be accessed.
"""
macro sample(def) @warn "@sample should be used inside a @sampler block." end
"""
    @first_stage(def)

Add a first stage model generation recipe to `stochasticprogram` using the syntax
```julia
@first_stage stochasticprogram::StochasticProgram = begin
    ...
end [defer]
```
where JuMP syntax is used inside the block to define the first stage model. During definition, the first stage model is referenced through the reserved keyword `model`.

Optionally, give the keyword `defer` after the  to delay generation of the first stage model.

## Examples

The following defines the first stage model given by:

```math
  minimize 100x₁ + 150x₂
    s.t  x₁ + x₂ ≤ 120
         x₁ ≥ 40
         x₂ ≥ 20
```

```julia
@first_stage sp = begin
    @variable(model, x₁ >= 40)
    @variable(model, x₂ >= 20)
    @objective(model, Min, 100*x₁ + 150*x₂)
    @constraint(model, x₁ + x₂ <= 120)
end
```

See also: [`@second_stage`](@ref)
"""
macro first_stage(arg) esc(:(@first_stage $arg generate)) end
macro first_stage(arg, defer)
    return esc(:(@stage 1 $arg $defer))
end
"""
    @second_stage(def)

Add a second stage model generation recipe to `stochasticprogram` using the syntax
```julia
@second_stage stochasticprogram::StochasticProgram = begin
    @decision var1 var2 ...
    ...
end [defer]
```
where JuMP syntax is used inside the block to define the second stage model. During definition, the second stage model is referenced through the reserved keyword `model`.

Optionally, give the keyword `defer` after the  to delay generation of the first stage model.

## Examples

The following defines the second stage model given by:
```math
  minimize q₁(ξ)y₁ + q₂(ξ)y₂
    s.t  6y₁ + 10y₂ ≤ 60x₁
         8y₁ + 5y₂ ≤ 60x₂
         0 ≤ y₁ ≤ d₁(ξ)
         0 ≤ y₂ ≤ d₂(ξ)
```
where ``q₁(ξ), q₂(ξ), d₁(ξ), d₂(ξ)`` depend on the scenario ``ξ`` and ``x₁, x₂`` are first stage variables. Two scenarios are added so that two second stage models are generated.

```julia
@second_stage sp = begin
    @decision x₁ x₂
    @uncertain q₁ q₂ d₁ d₂
    @variable(model, 0 <= y₁ <= d₁)
    @variable(model, 0 <= y₂ <= d₂)
    @objective(model, Min, q₁*y₁ + q₂*y₂)
    @constraint(model, 6*y₁ + 10*y₂ <= 60*x₁)
    @constraint(model, 8*y₁ + 5*y₂ <= 80*x₂)
end
```

See also: [`@first_stage`](@ref)
"""
macro second_stage(arg) esc(:(@second_stage $arg generate)) end
macro second_stage(arg, defer)
    return esc(:(@stage 2 $arg $defer))
end
"""
    @parameters(def)

Define the problem parameters in a @stage block
```julia
@parameters param1, param2, ...
```
possibly with default values. Any defined parameter without a default value must be supplied as a keyword argument to [`instantiate`](@ref) or [`SAA`](@ref) when creating models.

## Examples

```julia
@parameters d

@parameters begin
    Crops = [:wheat, :corn, :beets]
    Cost = Dict(:wheat=>150, :corn=>230, :beets=>260)
    Budget = 500
end
```

See also [`@decision`](@ref), [`@uncertain`](@ref), [`@stage`](@ref)
"""
macro parameters(def) @warn "@parameters should be used inside a @stage block." end
"""
    @decision(def)

In a @stage block, annotate each decision taken in the previous stage using the syntax
```julia
@decision var1, var2, ...
```

## Examples

```julia
@decision x₁, x₂
```

See also [`@parameters`](@ref), [`@uncertain`](@ref), [`@stage`](@ref)
"""
macro decision(def) @warn "@decision should be used inside a @stage block." end
"""
    @uncertain(def)

In a @stage block, annotate each uncertain variable using the syntax
```julia
@uncertain var1, var2, ...
```
This assumes that the [`Scenario`] type is used. Alternatively, user-defined scenarios can be specified by annotating the type. Also, inside a @stochastic_model block, user-defined scenarios can be created during the @uncertain annotation, following [`@scenario`](@ref).

## Examples

```julia
@uncertain q₁ q₂ d₁ d₂

@scenario Simple = begin
    q₁::Float64
    q₂::Float64
    d₁::Float64
    d₂::Float64
end
@uncertain ξ::SimpleScenario

@stochastic_model begin
    ...
    @uncertain ξ::SimpleScenario = begin
        q₁::Float64
        q₂::Float64
        d₁::Float64
        d₂::Float64
    end
    ...
end
```

See also [`@scenario`](@ref), [`@parameters`](@ref), [`@decision`](@ref), [`@stage`](@ref)
"""
macro uncertain(def) @warn "@uncertain should be used inside a @stage block." end
"""
    @stage(def)

Add a stage model generation recipe to `stochasticprogram` using the syntax
```julia
@stage stage stochasticprogram::StochasticProgram = begin
    @parameters param1 param2 ...
    @decision var1 var2 ...
    @uncertain ξ
    ... JuMPdef ...
    ...
end [defer]
```
where JuMP syntax is used inside the block to define the stage model. During definition, the second stage model is referenced through the reserved keyword `model`.

Optionally, give the keyword `defer` after the  to delay generation of the stage model.

## Examples

The following defines the first stage model given by:
```math
  minimize 100x₁ + 150x₂
    s.t  x₁ + x₂ ≤ 120
         x₁ ≥ 40
         x₂ ≥ 20
```
and the second-stage model given by:
```math
  minimize q₁(ξ)y₁ + q₂(ξ)y₂
    s.t  6y₁ + 10y₂ ≤ 60x₁
         8y₁ + 5y₂ ≤ 60x₂
         0 ≤ y₁ ≤ d₁(ξ)
         0 ≤ y₂ ≤ d₂(ξ)
```
where ``q₁(ξ), q₂(ξ), d₁(ξ), d₂(ξ)`` depend on the scenario ``ξ`` and ``x₁, x₂`` are first stage variables. Two scenarios are added so that two second stage models are generated.

```jldoctest
ξ₁ = Scenario(q₁ = -24.0, q₂ = -28.0, d₁ = 500.0, d₂ = 100.0, probability = 0.4)
ξ₂ = Scenario(q₁ = -28.0, q₂ = -32.0, d₁ = 300.0, d₂ = 300.0, probability = 0.6)

sp = StochasticProgram([ξ₁, ξ₂])

@stage 1 sp = begin
    @variable(model, x₁ >= 40)
    @variable(model, x₂ >= 20)
    @objective(model, Min, 100*x₁ + 150*x₂)
    @constraint(model, x₁ + x₂ <= 120)
end

@stage 2 sp = begin
    @decision x₁ x₂
    @uncertain q₁ q₂ d₁ d₂
    @variable(model, 0 <= y₁ <= d₁)
    @variable(model, 0 <= y₂ <= d₂)
    @objective(model, Min, q₁*y₁ + q₂*y₂)
    @constraint(model, 6*y₁ + 10*y₂ <= 60*x₁)
    @constraint(model, 8*y₁ + 5*y₂ <= 80*x₂)
end

# output

Stochastic program with:
 * 2 decision variables
 * 2 recourse variables
 * 2 scenarios of type Scenario
Solver is default solver

```

See also: [`@parameters`](@ref), [`@decision`](@ref), [`@uncertain`](@ref)
"""
macro stage(stage, arg) esc(:(@stage $stage $arg generate)) end
macro stage(stage, args, defer)
    generate = if defer == :defer
        :false
    elseif defer == :generate
        :true
    else
        error("Unknown option $defer")
    end
    @capture(args, sp_Symbol = def_) || error("Invalid syntax. Expected stage, multistage = begin JuMPdef end")
    # Save variable definitions separately
    vardefs = Expr(:block)
    for line in block(def).args
        if @capture(line, @constraint(m_Symbol, constdef__)) || @capture(line, @objective(m_Symbol, objdef__)) || @capture(line, @parameters args__) || @capture(line, @decision args__) || @capture(line, @uncertain args__)
            continue
        else
            push!(vardefs.args, line)
        end
    end
    # Handle parameters
    def = postwalk(def) do x
        if @capture(x, @parameters arg_)
            code = Expr(:block)
            for paramdef in prettify(arg).args
                if @capture(paramdef, key_Symbol = val_) || @capture(paramdef, key_Symbol)
                    push!(code.args, :($key = stage.$key))
                else
                    error("Incorrect @parameters syntax. Specify parameter names, possibly with a default value.")
                end
            end
            pushfirst!(vardefs.args, code)
            return code
        elseif @capture(x, @parameters args__)
            code = Expr(:block)
            for param in args
                push!(code.args, :($param = stage.$param))
            end
            pushfirst!(vardefs.args, code)
            return code
        elseif @capture(x, @decision args__)
            stage == 1 && error("@decision links cannot be used in the first stage.")
            code = Expr(:block)
            for var in args
                varkey = Meta.quot(var)
                decisiondef = @q begin
                    haskey(parent.objDict, $varkey) || error("Decision $($varkey) not present in parent model.")
                    $var = parent.objDict[$varkey]
                end
                push!(code.args, decisiondef)
            end
            return code
        elseif @capture(x, @uncertain var_Symbol::t_Symbol)
            stage == 1 && error("@uncertain declarations cannot be used in the first stage.")
            return :($var::$t = scenario)
        elseif @capture(x, @uncertain vars__ from scenvar_Symbol::t_Symbol)
            stage == 1 && error("@uncertain declarations cannot be used in the first stage.")
            code = @q begin
                $scenvar::$t = scenario
            end
            for var in vars
                varkey = Meta.quot(var)
                push!(code.args, :($var = getproperty($scenvar, $varkey)))
            end
            return code
        elseif @capture(x, @uncertain vars__ from t_Symbol)
            stage == 1 && error("@uncertain declarations cannot be used in the first stage.")
            code = @q begin
                scenvar::$t = scenario
            end
            for var in vars
                varkey = Meta.quot(var)
                push!(code.args, :($var = getproperty(scenvar, $varkey)))
            end
            return code
        elseif @capture(x, @uncertain var_Symbol::t_Symbol = scenariodef_)
            stage == 1 && error("@uncertain declarations cannot be used in the first stage.")
            return @q begin
                $var::$t = scenario
            end
        elseif @capture(x, @uncertain vars__ from scenvar_Symbol::t_Symbol = scenariodef_)
            stage == 1 && error("@uncertain declarations cannot be used in the first stage.")
            code = @q begin
                $scenvar::$t = scenario
            end
            for var in vars
                varkey = Meta.quot(var)
                push!(code.args, :($var = getproperty($scenvar, $varkey)))
            end
            return code
        elseif @capture(x, @uncertain vars__ from t_Symbol = scenariodef_)
            stage == 1 && error("@uncertain declarations cannot be used in the first stage.")
            code = @q begin
                scenvar::$t = scenario
            end
            for var in vars
                varkey = Meta.quot(var)
                push!(code.args, :($var = getproperty(scenvar, $varkey)))
            end
            return code
        elseif @capture(x, @uncertain args__)
            stage == 1 && error("@uncertain declarations cannot be used in the first stage.")
            code = Expr(:block)
            for var in args
                varkey = Meta.quot(var)
                push!(code.args, :($var = scenario[$varkey]))
            end
            return code
        else
            return x
        end
    end
    # The first stage is treated different from the subsequent stages
    generatordefs = if stage == 1
        @q begin
            $(esc(sp)).generator[:stage_1] = ($(esc(:model))::JuMP.Model, $(esc(:stage))) -> begin
                $(esc(def))
	        return $(esc(:model))
            end
        end
    else
        @q begin
            $(esc(sp)).generator[Symbol(:stage_,$stage)] = ($(esc(:model))::JuMP.Model, $(esc(:stage)), $(esc(:scenario))::AbstractScenario, $(esc(:parent))::JuMP.Model) -> begin
                $(esc(def))
	        return $(esc(:model))
            end
        end
    end
    # Create definition code
    code = @q begin
        isa($(esc(sp)), StochasticProgram) || error("Given object is not a stochastic program.")
        $(esc(sp)).generator[Symbol(:stage_,$stage,:_vars)] = ($(esc(:model))::JuMP.Model, $(esc(:stage))) -> begin
            $(esc(vardefs))
	    return $(esc(:model))
        end
        $generatordefs
        if $generate
            generate_stage!($(esc(sp)), $stage)
        end
        $(esc(sp))
    end
    # Return code
    return prettify(code)
end
# ========================== #


text:
# An `Environment` object stores information about the environment in which a
# suite of benchmarks were executed.
#
# Fields:
#
#     uuid::String: A random UUID that uniquely identifies each run.
#
#     timestamp::String: The time when we began executing benchmarks.
#
#     julia_sha1::String: The SHA1 for the Julia Git revision we're working
#         from.
#
#     package_sha1::Nullable{String}: The SHA1 for the current repo's Git
#         revision (if any). This field is null when the code was not executed
#         inside of a Git repo.
#
#     os::String: The OS we're running on.
#
#     cpu_cores::Int: The number of CPU cores available.
#
#     arch::String: The architecture we're running on.
#
#     machine::String: The machine type we're running on.
#
#     use_blas64::Bool: Was BLAS configured to use 64-bits?
#
#     word_size::Int: The word size of the host machine.

immutable Environment
    uuid::Compat.String
    timestamp::Compat.String
    julia_sha1::Compat.String
    package_sha1::Nullable{Compat.String}
    os::Compat.String
    cpu_cores::Int
    arch::Compat.String
    machine::Compat.String
    use_blas64::Bool
    word_size::Int

    function Environment()
        uuid = string(Base.Random.uuid4())
        timestamp = Libc.strftime("%Y-%m-%d %H:%M:%S", round(Int, time()))
        julia_sha1 = Base.GIT_VERSION_INFO.commit
        package_sha1 = Nullable{Compat.String}()
        try
            sha1 = readchomp(pipeline(`git rev-parse HEAD`, stderr=Base.DevNull))
            package_sha1 = Nullable{Compat.String}(Compat.String(sha1))
        end
        os = string(Compat.KERNEL === :NT ? :Windows : Compat.KERNEL)
        cpu_cores = Sys.CPU_CORES
        arch = string(Sys.ARCH)
        machine = Base.MACHINE
        use_blas64 = Base.USE_BLAS64
        word_size = Sys.WORD_SIZE

        new(
            uuid,
            timestamp,
            julia_sha1,
            package_sha1,
            os,
            cpu_cores,
            arch,
            machine,
            use_blas64,
            word_size,
        )
    end
end

# Pretty-print information about the environment in which benchmarks are being
# executed
#
# Arguments:
#
#     io::IO: An `IO` object to be written to.
#
#     e::Environment: The `Environment` object that we want to print to `io`.

function Base.show(io::IO, e::Environment)
    names = Compat.String[
        "UUID",
        "Time",
        "Julia SHA1",
        "Package SHA1",
        "Machine kind",
        "CPU architecture",
        "CPU cores",
        "OS",
        "Word size",
        "64-bit BLAS",
    ]
    values = Any[
        e.uuid,
        e.timestamp,
        e.julia_sha1,
        get(e.package_sha1, "NULL"),
        e.machine,
        e.arch,
        e.cpu_cores,
        e.os,
        e.word_size,
        e.use_blas64,
    ]

    @printf(io, "================== Benchmark Environment ==================\n")
    max_length = maximum([length(n) for n in names]) +  1
    for (n, v) in zip(names, values)
        @printf(io, "%s: %s\n", lpad(n, max_length), v)
    end
end

# Log information about the environment in which benchmarks are being executed
# to a CSV file.
#
# Arguments:
#
#     filename::AbstractString: The name of a file to which we'll write information
#         about the `Environment` object, `e`.
#
#     e::Environment: The `Environment` that we want to write to disk.
#
#     append::Bool: Should we append to an existing file or create a new one?
#         Defaults to false.

function Base.writecsv(filename::AbstractString, e::Environment, append::Bool = false)
    names = Compat.String[
        "uuid",
        "timestamp",
        "julia_sha1",
        "package_sha1",
        "os",
        "cpu_cores",
        "arch",
        "machine",
        "use_blas64",
        "word_size",
    ]
    values = Any[
        e.uuid,
        e.timestamp,
        e.julia_sha1,
        get(e.package_sha1, "NULL"),
        e.os,
        string(e.cpu_cores),
        e.arch,
        e.machine,
        string(e.use_blas64),
        string(e.word_size),
    ]

    if append
        io = open(filename, "a")
    else
        io = open(filename, "w")
    end
    println(io, join(names, ','))
    println(io, join(values, ','))
    close(io)
end


text:
using Distributions, Logging

# There is some runtime code generation, so the easiest way
# to make sure everything we need is precompiled is just to
# run a super (super) quick fit during precompilation.
let

    logger = SimpleLogger(stdout, Logging.Error)
    with_logger(logger) do

        planet = DirectDetections.Planet(
            Priors(
                a = Uniform(1, 50),
                e = Beta(1.2, 5),
                τ = Uniform(0, 1),
                ω = Uniform(0, 2pi),
                i = Uniform(0, pi),
                Ω = Uniform(0, pi),
                mass = Uniform(0, 1),
            ),
            Astrometry(
                (epoch=1234.0, ra=123., dec=123., σ_ra=12., σ_dec=34.),
            ),
            name=:test
        )

        pma = ProperMotionAnom(
            (;
                ra_epoch=1234.,
                dec_epoch=1234.,
                pm_ra=12.,
                σ_pm_ra=.34,
                pm_dec=12.,
                σ_pm_dec=.34,
            ),
        )
        system = System(
            Priors(
                μ = TruncatedNormal(1.0, 0.2, 0.1, Inf),
                plx = TruncatedNormal(12.0, 0.2, 0.1, Inf),
            ),
            pma,
            planet,
            name=:Test
        )

        # This can't actually start since AdvancedHMC uses
        # Revise to load in ForwardDiff. It's not available
        # during precompile.

        # DirectDetections.hmc(
        #     system, 0.85,
        #     adaptation = 5,
        #     iterations = 10,
        #     tree_depth = 5,
        #     initial_samples = 5,
        # )

        # let's trigger some things manually though
        initial_θ_0 = sample_priors(system)
        ln_prior_transformed = make_ln_prior_transformed(system)
        arr2nt = DirectDetections.make_arr2nt(system)

    end
    nothing
end


text:
# using LoopVectorization, Test
# T = Float32

@testset "broadcast" begin
    M, N = 37, 47
    # M = 77;
    for T ∈ (Float32, Float64, Int32, Int64)
        @show T, @__LINE__
        R = T <: Integer ? (T(-100):T(100)) : T
        a = rand(R,99,99,99);
        b = rand(R,99,99,1);
        bl = LowDimArray{(true,true,false)}(b);
        br = reshape(b, (99,99));
        c1 = a .+ b;
        c2 = @avx a .+ bl;
        @test c1 ≈ c2
        fill!(c2, 99999); @avx c2 .= a .+ br;
        @test c1 ≈ c2
        fill!(c2, 99999); @avx c2 .= a .+ b;
        @test c1 ≈ c2
        br = reshape(b, (99,1,99));
        bl = LowDimArray{(true,false,true)}(br);
        @. c1 = a + br;
        fill!(c2, 99999); @avx @. c2 = a + bl;
        @test c1 ≈ c2
        fill!(c2, 99999); @avx @. c2 = a + br;
        @test c1 ≈ c2
        br = reshape(b, (1,99,99));
        bl = LowDimArray{(false,)}(br);
        @. c1 = a + br;
        fill!(c2, 99999);
        @test c1 ≈ @avx @. c2 = a + bl
        # @test c1 ≈ c2
        br = reshape(rand(R,99), (1,99,1));
        bl = LowDimArray{(false,)}(br);
        @. c1 = a + br;
        fill!(c2, 99999);
        @avx @. c2 = a + bl;
        @test c1 ≈ c2

        if T <: Integer
            xs = rand(-T(100):T(100), M);
        else
            xs = rand(T, M);
        end
        max_ = maximum(xs, dims=1);
        @test (@avx exp.(xs .- LowDimArray{(false,)}(max_))) ≈ exp.(xs .- LowDimArray{(false,)}(max_))

        if T === Int32
            a = rand(T(1):T(100), 73, 1);
            @test sqrt.(Float32.(a)) ≈ @avx sqrt.(a)
        elseif T === Int64
            a = rand(T(1):T(100), 73, 1);
            @test sqrt.(a) ≈ @avx sqrt.(a)
        else
            a = rand(T, 73, 1);
            @test sqrt.(a) ≈ @avx sqrt.(a)
        end

        a = rand(R, M); B = rand(R, M, N); c = rand(R, N); c′ = c';
        d1 =      @. a + B * c′;
        d2 = @avx @. a + B * c′;
        @test d1 ≈ d2

        @.      d1 = a + B * c′;
        @avx @. d2 = a + B * c′;
        @test d1 ≈ d2

        d3 = a .+ B * c;
        d4 = @avx a .+ B *ˡ c;
        @test d3 ≈ d4

        fill!(d3, -1000.0);
        fill!(d4, 91000.0);

        d3 .= a .+ B * c;
        @avx d4 .= a .+ B *ˡ c;
        @test d3 ≈ d4

        fill!(d4, 91000.0);
        @avx @. d4 = a + B *ˡ c;
        @test d3 ≈ d4

        M, K, N = 77, 83, 57;
        A = rand(R,M,K); B = rand(R,K,N); C = rand(R,M,N);
        At = copy(A');
        D1 = C .+ A * B;
        D2 = @avx C .+ A .*ˡ B;
        @test D1 ≈ D2
        if RUN_SLOW_TESTS
            fill!(D2, -999999); D2 = @avx C .+ At' *ˡ B;
            @test D1 ≈ D2
            fill!(D2, -999999); @test A * B ≈ (@avx @. D2 = A *ˡ B)
            D1 .= view(C, 1, :)' .+ A * B;
            fill!(D2, -999999);
            @avx D2 .= view(C, 1, :)' .+ A .*ˡ B;
            @test D1 ≈ D2
            C3d = rand(R,3,M,N);
            D1 .= view(C3d, 1, :, :) .+ A * B;
            fill!(D2, -999999);
            @avx D2 .= view(C3d, 1, :, :) .+ A .*ˡ B;
            @test D1 ≈ D2
        end
        D1 .= 9999;
        @avx D2 .= 9999;
        @test D1 == D2
        D1 .= -99999;
        @avx D2' .= -99999;
        @test D1 == D2

        if VERSION > v"1.2"
            b = rand(T,K); x = rand(R,N);
            D1 .= C .+ A * (b .+ x');
            @avx @. D2 = C + A *ˡ (b + x');
            @test D1 ≈ D2
            D2 = @avx @. C + A *ˡ (b + x');
            @test D1 ≈ D2
        end

        if T <: Union{Float32,Float64}
            D3 = cos.(B');
            D4 = @avx cos.(B');
            @test D3 ≈ D4

            fill!(D3, -1e3); fill!(D4, 9e9);
            Bt = transpose(B);
            @. D3 = exp(Bt);
            @avx @. D4 = exp(Bt);
            @test D3 ≈ D4

            D1 = similar(B); D2 = similar(B);
            D1t = transpose(D1);
            D2t = transpose(D2);
            @. D1t = exp(Bt);
            @avx @. D2t = exp(Bt);
            @test D1t ≈ D2t

            fill!(D1, -1e3);
            fill!(D2, 9e9);
            @. D1' = exp(Bt);
            lset = @avx @. D2' = exp(Bt);

            @test D1 ≈ D2

            a = rand(137);
            b1 = @avx @. 3*a + sin(a) + sqrt(a);
            b2 =      @. 3*a + sin(a) + sqrt(a);
            @test b1 ≈ b2
            three = 3; fill!(b1, -9999);
            @avx @. b1 = three*a + sin(a) + sqrt(a);
            @test b1 ≈ b2

            C = rand(100,10,10);
            D1 = C .^ 0.3;
            D2 = @avx C .^ 0.3;
            @test D1 ≈ D2
            @. D1 = C ^ 2;
            @avx @. D2 = C ^ 2;
            @test D1 ≈ D2
        end
    end
end


text:
using Spike
using Test

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

max_stars_repo_name: catmanium/Spike

max_issues_repo_name: catmanium/Spike

max_forks_repo_name: catmanium/Spike


text:
using DataDepsGenerators
using Test

using ReferenceTests

@testset "TERN test" begin
    @test_reference "../references/TERN/Influenza A.txt" generate(TERN(), "https://dataone.tern.org.au/mn/v2/object/aekos.org.au/collection/nsw.gov.au/nsw_atlas/vis_flora_module/KAHDRAIN.20150515")
end


text:
@testset "anti/diagonal attack" begin
    occ = [0,0,0,0,0,0,0,0,
           0,0,1,0,0,0,1,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,1,0,
           0,1,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0]

    att = [0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,1,0,1,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,1,0,1,0,0,
           0,0,1,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0]

    edg = [0,0,0,0,0,0,0,0,
           0,0,1,0,0,0,1,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,1,0,
           0,1,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0]
    occ_uint = Bobby.cvt_to_uint(BitArray(occ))
    ui = Bobby.INT2UINT[29]
    m, e = Bobby.cross_attack(occ_uint, ui)
    @test all(Int.(Bobby.cvt_to_bitarray(m)) .== att)
    @test all(Int.(Bobby.cvt_to_bitarray(e)) .== edg)
end

@testset "bishop fen" begin
    @test test_fen("k7/8/8/8/8/8/8/B6K w Q - 0 1", 1, [10])
    @test test_fen("k7/8/8/8/8/1p6/1P6/B6K w Q - 0 1", 1, [3])
    @test test_fen("k7/8/8/8/8/1p6/1p6/B6K w Q - 0 1", 1, [4])
    @test test_fen("k7/8/8/8/7q/7B/8/7K w Q - 0 1", 1, [3])
    @test test_fen("1k6/8/8/8/8/7q/6B1/7K w - -", 1, [2])
    @test test_fen("k7/8/8/8/8/8/8/6BK w Q - 0 1", 1, [9])
    @test test_fen("kb6/8/8/8/8/8/8/7K b Q - 0 1", 1, [9])
    @test test_fen("k7/b7/Q7/8/8/8/8/7K b Q - 0 1", 1, [1])
    @test test_fen("k1b5/8/QR6/8/8/8/8/7K b Q - 0 1", 1, [1])
    @test test_fen("k7/1b6/QR6/8/8/8/8/6K1 b - -", 1, [2])
end


text:
module TestConfigurations

using BenchmarkConfigSweeps
using BenchmarkConfigSweeps.Internal: asconfig
using Test

const BCS = BenchmarkConfigSweeps

function test_asconfig()
    @test asconfig(BCS.nthreads(2)).nthreads == BCS.nthreads(2)
    @test asconfig((BCS.nthreads(3), BCS.nthreads(2))).nthreads == BCS.nthreads(2)
    @test asconfig((BCS.env("A" => 0), BCS.nthreads(2))).nthreads == BCS.nthreads(2)
    @test asconfig((BCS.env("A" => 0), BCS.nthreads(2))).env == BCS.env("A" => 0)
    @test asconfig((BCS.env("A" => 0, "B" => false), BCS.env("A" => 1))).env ==
          BCS.env("A" => 1, "B" => false)
end

function test_readme()
    nthreads_list = [1, 4, 8, 16]
    configs = Iterators.product(
        zip(
            BenchmarkConfigSweeps.nthreads.(nthreads_list),
            BenchmarkConfigSweeps.env.("OPENBLAS_NUM_THREADS" .=> nthreads_list),
        ),
        BenchmarkConfigSweeps.env.(
            "JULIA_PROJECT" .=> ["baseline", "target"],
            "JULIA_LOAD_PATH" => "@",  # not varied
        ),
        BenchmarkConfigSweeps.julia.(["julia1.5", "julia1.6"]),
    )

    configs = map(asconfig, configs)
    @test length(configs) == 4 * 2 * 2
    @test configs[1].nthreads == BCS.nthreads(1)
    @test configs[1].env == BCS.env(
        "OPENBLAS_NUM_THREADS" => 1,
        "JULIA_PROJECT" => "baseline",
        "JULIA_LOAD_PATH" => "@",
    )
    @test configs[1].julia == BCS.julia("julia1.5")
end

end  # module


text:
using Encoref
using Test, ParquetFiles, DataFrames

include("testutils.jl")

function load_result(dir::String)
    res_dir = joinpath(dir, "result")
    Dict(
        split(resf, ".")[1] => load(joinpath(res_dir, resf)) |> DataFrame |> df -> hcat(df.i1, df.i2)
        for resf in readdir(res_dir)
    )
end

max_level = 0

test_dirs = ["test_data/test_sys_$i" for i in 0:0]

@testset "Encoref.jl" begin
    @test 1 == 1
end

@testset "step-init" begin

end

@testset "e2e" begin

    for test_dir in test_dirs
        true_res = load_result(test_dir)
        @test true_res == Dict(k=>sortslices(v, dims=1) for (k,v) in coref(CorefSystem(test_dir)))
    end

end


text:
"""
$(TYPEDEF)

Defines a time dependent Hamiltonian object with dense Matrices.

# Fields

$(FIELDS)
"""
struct DenseHamiltonian{T<:Number} <: AbstractDenseHamiltonian{T}
    "List of time dependent functions"
    f::Vector
    "List of constant matrices"
    m::Vector
    "Internal cache"
    u_cache::AbstractMatrix
    "Size"
    size::Tuple
    "Eigen decomposition routine"
    EIGS::Any
end

"""
$(SIGNATURES)

Constructor of the `DenseHamiltonian` type. `funcs` and `mats` are lists of time-dependent functions and the corresponding matrices. The Hamiltonian can be represented as ``∑ᵢfuncs[i](s)×mats[i]``. `unit` specifies wether `:h` or `:ħ` is set to one when defining `funcs` and `mats`. The `mats` will be scaled by ``2π`` if unit is `:h`.

`EIGS` is the initializer for the eigendecomposition routine for the Hamiltonian. It should return a function of the signature: `(H, s, lvl) -> (w, v)`. The default method `EIGEN_DEFAULT` will use `LAPACK` routine.
"""
function DenseHamiltonian(funcs, mats; unit = :h, EIGS = EIGEN_DEFAULT)
    if any((x) -> size(x) != size(mats[1]), mats)
        throw(ArgumentError("Matrices in the list do not have the same size."))
    end
    if is_complex(funcs, mats)
        mats = complex.(mats)
    end
    hsize = size(mats[1])
    # use static array for size smaller than 100
    if hsize[1] <= 10
        mats = [SMatrix{hsize[1],hsize[2]}(unit_scale(unit) * m) for m in mats]
    else
        mats = unit_scale(unit) * mats
    end
    cache = similar(mats[1])
    EIGS = EIGS(cache)
    DenseHamiltonian{eltype(mats[1])}(funcs, mats, cache, hsize, EIGS)
end

"""
    function (h::DenseHamiltonian)(s::Real)

Calling the Hamiltonian returns the value ``2πH(s)``. The argument `s` is in the dimensionless time. The returned matrix is in the unit of angular frequency.
"""
function (h::DenseHamiltonian)(s::Real)
    fill!(h.u_cache, 0.0)
    for i = 1:length(h.f)
        @inbounds axpy!(h.f[i](s), h.m[i], h.u_cache)
    end
    h.u_cache
end

# The argument `p` is not essential for `DenseHamiltonian`
# It exists to keep the `update_cache!` interface consistent across
# all `AbstractHamiltonian` types
function update_cache!(cache, H::DenseHamiltonian, p, s::Real)
    fill!(cache, 0.0)
    for i = 1:length(H.m)
        @inbounds axpy!(-1.0im * H.f[i](s), H.m[i], cache)
    end
end

function update_vectorized_cache!(cache, H::DenseHamiltonian, p, s::Real)
    hmat = H(s)
    iden = one(hmat)
    cache .= 1.0im * (transpose(hmat) ⊗ iden - iden ⊗ hmat)
end

function (h::DenseHamiltonian)(du, u::AbstractMatrix, p, s::Real)
    fill!(du, 0.0 + 0.0im)
    H = h(s)
    gemm!('N', 'N', -1.0im, H, u, 1.0 + 0.0im, du)
    gemm!('N', 'N', 1.0im, u, H, 1.0 + 0.0im, du)
end

function Base.convert(S::Type{T}, H::DenseHamiltonian{M}) where {T<:Complex,M}
    mats = [convert.(S, x) for x in H.m]
    cache = similar(H.u_cache, complex{M})
    DenseHamiltonian{eltype(mats[1])}(H.f, mats, cache, size(H), H.EIGS)
end

function Base.convert(S::Type{T}, H::DenseHamiltonian{M}) where {T<:Real,M}
    f_val = sum((x) -> x(0.0), H.f)
    if !(typeof(f_val) <: Real)
        throw(TypeError(:convert, "H.f", Real, typeof(f_val)))
    end
    mats = [convert.(S, x) for x in H.m]
    cache = similar(H.u_cache, real(M))
    DenseHamiltonian{eltype(mats[1])}(H.f, mats, cache, size(H), H.EIGS)
end


text:
module SpectraFit

ENV["GKS_ENCODING"] = "utf-8"

include("Quadrupolar.jl")
include("ChemicalShift.jl")
include("utility_functions.jl")
include("optim_estimation.jl")
include("bayesian_estimation.jl")
include("genetic_algorithm.jl")
include("bboptim_estimation.jl")

export get_experimental,
       Quadrupolar,
       ChemicalShift,
       fit_quadrupolar,
       fit_chemical_shift,
       metropolis_hastings,
       compare_theoreticals,
       get_data,
       plot_chemical_shift,
       get_output_table,
       genetic_algorithm,
       genetic_algorithm_csa,
       fit_quadrupolar_bb,
       fit_chemicalshift_bb

end # module


text:
module Diagnostics

export
    HorizontalAverage, Timeseries, FieldMaximum,
    CFL, AdvectiveCFL, DiffusiveCFL, NaNChecker

using Oceananigans,
      Oceananigans.Operators

using Oceananigans: AbstractDiagnostic

include("diagnostics_kernels.jl")
include("nan_checker.jl")
include("horizontal_average.jl")
include("timeseries.jl")
include("field_maximum.jl")
include("cfl.jl")

end


text:
# Flavours
abstract type Flavour end
struct VanillaFlavour <: Flavour end
struct RockyRoadFlavour <: Flavour end

Base.show(io::IO, f::VanillaFlavour) = print(io, "Vanilla flavour")
Base.show(io::IO, f::RockyRoadFlavour) = print(io, "Rocky Road flavour")

Vanilla = VanillaFlavour()
RockyRoad = RockyRoadFlavour()

export Vanilla, RockyRoad

struct V2013R{F<:Flavour} <: Version
    flavour::F
end

Base.show(io::IO, v::V2013R) = print(io, "v2013R ($(v.flavour))")

"""
    v2013R([flavour])

Identifier for the 2013R version of the model.
Defaults to the `Vanilla` flavour, but also accepts `RockyRoad`.

# Examples
```jldoctest
julia> v2013R()
v2013R (Vanilla flavour)
```

```jldoctest
julia> v2013R(Vanilla)
v2013R (Vanilla flavour)
```

```jldoctest
julia> v2013R(RockyRoad)
v2013R (Rocky Road flavour)
```
"""
function v2013R(flavour::F = VanillaFlavour()) where {F<:Flavour}
    V2013R{F}(flavour)
end

export v2013R

@extend struct VariablesV2013 <: Variables
    Eind::Array{VariableRef,1} # Industrial emissions (GtCO2 per year)
    Ω::Array{VariableRef,1} # Damages as fraction of gross output
    Λ::Array{VariableRef,1} # Cost of emissions reductions  (trillions 2005 USD per year)
    CPRICE::Array{VariableRef,1} # Carbon price (2005$ per ton of CO2)
    CEMUTOTPER::Array{VariableRef,1} # Period utility
end

function model_vars(version::V2013R{F}, model::Model, N::Int64, cca_ubound::Float64, μ_ubound::Array{Float64,1}, Tat_ubound::Float64, cprice_ubound::Array{Float64,1}) where F <: Flavour
    # Variables #
    @variable(model, 0.0 <= μ[i=1:N] <= μ_ubound[i]); # Emission control rate GHGs
    @variable(model, FORC[1:N]); # Increase in radiative forcing (watts per m2 from 1900)
    @variable(model, 0.0 <= Tₐₜ[1:N] <= Tat_ubound); # Increase temperature of atmosphere (degrees C from 1900)
    @variable(model, -1.0 <= Tₗₒ[1:N] <= 20.0); # Increase temperatureof lower oceans (degrees C from 1900)
    @variable(model, Mₐₜ[1:N] >= 10.0); # Carbon concentration increase in atmosphere (GtC from 1750)
    @variable(model, Mᵤₚ[1:N] >= 100.0); # Carbon concentration increase in shallow oceans (GtC from 1750)
    @variable(model, Mₗₒ[1:N] >= 1000.0); # Carbon concentration increase in lower oceans (GtC from 1750)
    @variable(model, E[1:N]); # Total CO2 emissions (GtCO2 per year)
    @variable(model, Eind[1:N]); # Industrial emissions (GtCO2 per year)
    @variable(model, C[1:N] >= 2.0); # Consumption (trillions 2005 US dollars per year)
    @variable(model, K[1:N] >= 1.0); # Capital stock (trillions 2005 US dollars)
    @variable(model, CPC[1:N] >= 0.01); #  Per capita consumption (thousands 2005 USD per year)
    @variable(model, I[1:N] >= 0.0); # Investment (trillions 2005 USD per year)
    @variable(model, S[1:N]); # Gross savings rate as fraction of gross world product
    @variable(model, RI[1:N]); # Real interest rate (per annum)
    @variable(model, Y[1:N] >= 0.0); # Gross world product net of abatement and damages (trillions 2005 USD per year)
    @variable(model, YGROSS[1:N] >= 0.0); # Gross world product GROSS of abatement and damages (trillions 2005 USD per year)
    @variable(model, YNET[1:N]); # Output net of damages equation (trillions 2005 USD per year)
    @variable(model, DAMAGES[1:N]); # Damages (trillions 2005 USD per year)
    @variable(model, Ω[1:N]); # Damages as fraction of gross output
    @variable(model, Λ[1:N]); # Cost of emissions reductions  (trillions 2005 USD per year)
    @variable(model, MCABATE[1:N]); # Marginal cost of abatement (2005$ per ton CO2)
    @variable(model, CCA[1:N] <= cca_ubound); # Cumulative industrial carbon emissions (GTC)
    @variable(model, PERIODU[1:N]); # One period utility function
    @variable(model, CPRICE[i=1:N] <= cprice_ubound[i]); # Carbon price (2005$ per ton of CO2)
    @variable(model, CEMUTOTPER[1:N]); # Period utility
    @variable(model, UTILITY); # Welfare function
    VariablesV2013(μ,FORC,Tₐₜ,Tₗₒ,Mₐₜ,Mᵤₚ,Mₗₒ,E,C,K,CPC,I,S,RI,Y,YGROSS,YNET,DAMAGES,MCABATE,CCA,PERIODU,UTILITY,Eind,Ω,Λ,CPRICE,CEMUTOTPER)
end

include("Results2013R.jl")
include("2013RVanilla.jl")
include("2013RRockyRoad.jl")


text:
# InclusionRepetition model is a two-part model like the hurdle but instead of a
# PositivePoisson multiplied by the Bernulli, it uses a regular Poisson for
# repetition (counts in excess of 1), conditional on inclusion.

"Selection returned object"
mutable struct InclusionRepetition{Z<:RegressionModel,P<:RegressionModel} <: TwoPartModel{Z,P}
  mzero::Z                # model for zeros
  mpos::P                 # model for positive counts
  fittedzero::Bool        # whether the model for zeros was fitted
  fittedpos::Bool         # whether the model for positives was fitted
end

# predicted counts vector given expected inclusion (μzero) and repetition (μpos)
μtpm(m::InclusionRepetition, μzero::V, μpos::V) where {T, V<:AbstractArray{T}} = μzero .* (one(T) .+ μpos)

function excessy!(ypos::V, ::Type{InclusionRepetition}) where {T, V<:AbstractVector{T}}
  @simd for i = 1:length(ypos)
    @inbounds ypos[i] -= one(T)
  end
  ypos
end

minypos(::Type{InclusionRepetition}) = 0.0

"""
    fit(InclusionRepetition,M,X,y; Xpos=Xpos, <keyword arguments>)

Fit an Inclusion-Repetition model of count vector y on X with potentially another
covariates matrix Xpos used to model positive counts (repetitions).

# Example with GLM:
```julia
  m = fit(InclusionRepetition,GeneralizedLinearModel,X,y; Xpos=Xpos)
  yhat = predict(m, X; Xpos=Xpos)
```

# Example with Lasso regularization:
```julia
  m = fit(InclusionRepetition,GammaLassoPath,X,y; Xpos=Xpos)
  yhat = predict(m, X; Xpos=Xpos, select=MinAICc())
```

# Arguments
- `M::RegressionModel`
- `counts` n-by-d matrix of counts (usually sparse)
- `dzero::UnivariateDistribution = Binomial()` distribution for zeros (inclusion) model
- `dpos::UnivariateDistribution = Poisson()` distribution for positives (repetition) model
- `lzero::Link=canonicallink(dzero)` link function for zeros model
- `lpos::Link=canonicallink(dpos)` link function for positives model

# Keywords
- `Xpos::Union{AbstractMatrix{T},Nothing} = nothing` covariates matrix for positives
  model or nothing to use X for both parts
- `dofit::Bool = true` fit the model or just construct its shell
- `wts::V = ones(y)` observation weights
- `offsetzero::AbstractVector = similar(y, 0)` offsets for zeros model
- `offsetpos::AbstractVector = similar(y, 0)` offsets for positives model
- `offset::AbstractVector = similar(y, 0)` offsets for both model parts
- `verbose::Bool=true`
- `showwarnings::Bool=false`
- `fitargs...` additional keyword arguments passed along to fit(M,...)
"""
function StatsBase.fit(::Type{InclusionRepetition},::Type{M},
  X::AbstractMatrix{T}, y::V,
  dzero::UnivariateDistribution = Binomial(),
  dpos::UnivariateDistribution = Poisson(),
  lzero::Link = canonicallink(dzero),
  lpos::Link = canonicallink(dpos);
  Xpos::Union{AbstractMatrix{T},Nothing} = nothing,
  dofit::Bool = true,
  wts::V = fill(one(eltype(y)),size(y)),
  offsetzero::AbstractVector = similar(y, 0),
  offsetpos::AbstractVector = similar(y, 0),
  offset::AbstractVector = similar(y, 0),
  verbose::Bool = false,
  showwarnings::Bool = false,
  fitargs...) where {M<:RegressionModel,T<:FP,V<:FPVector}

  ixpos, Iy = getIy(y)

  offsetzero, offsetpos = setoffsets(y, ixpos, offset, offsetzero, offsetpos)

  mzero, fittedzero = fitzero(M, X, Iy, dzero, lzero, dofit, wts, offsetzero, verbose, showwarnings, fitargs...)

  # Xpos optional argument allows to specify a data matrix only for positive counts
  if Xpos == nothing
    # use X for Xpos too
    @inbounds Xpos = X[ixpos,:]
  elseif size(Xpos,1) == length(y)
    # Xpos has same dimensions as X, take only positive y ones
    @inbounds Xpos = Xpos[ixpos,:]
  end

  @inbounds ypos = y[ixpos]
  @inbounds wtspos = wts[ixpos]

  mpos, fittedpos = fitpos(InclusionRepetition, M, Xpos, ypos, dpos, lpos, dofit, wtspos, offsetpos, verbose, showwarnings, fitargs...)

  InclusionRepetition(mzero,mpos,fittedzero,fittedpos)
end

"""
    fit(InclusionRepetition,M,f,df; fpos=Xpos, <keyword arguments>)

Takes dataframe and two formulas, one for each model part. Otherwise same arguments
as [`fit(::InclusionRepetition)`](@ref)

# Example
```julia
  fit(InclusionRepetition,GeneralizedLinearModel,@formula(y ~ x1*x2), df; fpos=@formula(y ~ x1*x2+x3))
```
"""
function StatsBase.fit(::Type{InclusionRepetition},::Type{M},
                      f::FormulaTerm,
                      data,
                      dzero::UnivariateDistribution = Binomial(),
                      dpos::UnivariateDistribution = Poisson(),
                      lzero::Link = canonicallink(dzero),
                      lpos::Link = canonicallink(dpos);
                      fpos::FormulaTerm = f,
                      dofit::Bool = true,
                      wts = fill(1.0,size(data,1)),
                      offsetzero = Float64[],
                      offsetpos = Float64[],
                      offset = Float64[],
                      verbose::Bool = false,
                      showwarnings::Bool = false,
                      contrasts::Dict{Symbol,<:Any} = Dict{Symbol,Any}(),
                      fitargs...) where {M<:RegressionModel}

  Tables.istable(data) || throw(ArgumentError("expected data in a Table, got $(typeof(data))"))

  cols = columntable(data)

  mfzero = ModelFrame(f, cols, model=M, contrasts=contrasts)
  mmzero = ModelMatrix(mfzero)
  y = response(mfzero)

  ixpos, Iy = getIy(y)

  offsetzero, offsetpos = setoffsets(y, ixpos, offset, offsetzero, offsetpos)

  # fit zero model to entire sample
  # TODO: should be wrapped in TableRegressionModel but can't figure out right now why it complains it is not defined
  # mzero = TableRegressionModel(fit(M, mmzero.m, Iy, dzero, lzero; wts=wts, offset=offsetzero, fitargs...), mfzero, mmzero)
  mzero, fittedzero = fitzero(M, mmzero.m, Iy, dzero, lzero, dofit, wts, offsetzero, verbose, showwarnings, fitargs...)

  mfpos = (f===fpos) ? mfzero : ModelFrame(fpos, cols, model=M, contrasts=contrasts)
  mmpos = (f===fpos) ? mmzero : ModelMatrix(mfpos)
  mmpos.m = mmpos.m[ixpos,:]

  mpos, fittedpos = fitpos(InclusionRepetition, M, mmpos.m, y[ixpos], dpos, lpos, dofit, wts[ixpos], offsetpos, verbose, showwarnings, fitargs...)

  InclusionRepetition(mzero,mpos,fittedzero,fittedpos)
end


text:
#=
  Fixed Point Convergence
  Fixed Point Existance
  Fixed Point Uniqueness
  Newtons Method
  Bisection Method
=#

function testConvergence(f,p0,TOL,N0)
  # takes starting point, tolerance, maximum iterations
  i = 1
  while i <= N0
    p = g(p0)
    if abs(p-p0) < TOL
      return true
    else
      i += 1
      p0 = p
    end
  end
  return false
  println("Failure: N0 tries")

function testExistance(f,a,b)
  # interval [a,b]
  if f(a) * f(b) <0
    return true
    #does this make asusmptions
  else
    return false
  end

function testUniqueness(f,a,b)
  # interval [a,b]
  if derivative(f,a) * derivative(f,b) < 0
    return true
    #does this make assumptions
  else
    return false
  end

function newtonsMethod(f,p0,TOL,N0)
  if testConvergence(f) == true
    i = 1
    while i <= N0
        p = p0 - (derivative(f,p0) / f(p0)
        if abs(p-p0) < TOL
          return p
        else
          i += 1
          p0 = p
        end
    end
  else
    #return eror e
  end

function secantMethod()

function bisectionMethod(f,a,b,DOA,NO)
  # interval [a,b], DOA degree of accuracy
  if testConvergence(f) == true
    i = 1
    while i <= N0
      m = (1/2)(a + b)
      if f(m) = 0 | abs(f(m)) < DOA
        return m
      elseif f(m) * f(a) < 0
        a = m
      elseif f(m) * f(b) < 0
        b = m
      else
        #return eror e
      end
    end
    return m
  else
    #return error e
  end

max_stars_repo_name: K-P-S/Julia-Scientific-Computing

max_issues_repo_name: K-P-S/Julia-Scientific-Computing

max_forks_repo_name: K-P-S/Julia-Scientific-Computing