text:
# types from SymEngine to Julia
## CSetBasic
type CSetBasic
    ptr::Ptr{Void}
end

function CSetBasic()
    z = CSetBasic(ccall((:setbasic_new, libsymengine), Ptr{Void}, ()))
    finalizer(z, CSetBasic_free)
    z
end

function CSetBasic_free(x::CSetBasic)
    if x.ptr != C_NULL
        ccall((:setbasic_free, libsymengine), Void, (Ptr{Void},), x.ptr)
        x.ptr = C_NULL
    end
end

function Base.length(s::CSetBasic)
    ccall((:setbasic_size, libsymengine), UInt, (Ptr{Void},), s.ptr)
end

function Base.getindex(s::CSetBasic, n::UInt)
    result = Basic()
    ccall((:setbasic_get, libsymengine), Void, (Ptr{Void}, UInt, Ptr{Basic}), s.ptr, n, &result)
    result
end

function Base.convert(::Type{Vector}, x::CSetBasic)
    n = Base.length(x)
    [x[i-1] for i in 1:n]
end
Base.convert(::Type{Set}, x::CSetBasic) = Set(convert(Vector, x))

## VecBasic Need this for get_args...

type CVecBasic
    ptr::Ptr{Void}
end

function CVecBasic()
    z = CVecBasic(ccall((:vecbasic_new, libsymengine), Ptr{Void}, ()))
    finalizer(z, CVecBasic_free)
    z
end

function CVecBasic_free(x::CVecBasic)
    if x.ptr != C_NULL
        ccall((:vecbasic_free, libsymengine), Void, (Ptr{Void},), x.ptr)
        x.ptr = C_NULL
    end
end

function Base.length(s::CVecBasic)
    ccall((:vecbasic_size, libsymengine), UInt, (Ptr{Void},), s.ptr)
end

function Base.getindex(s::CVecBasic, n::UInt)
    result = Basic()
    ccall((:vecbasic_get, libsymengine), Void, (Ptr{Void}, UInt, Ptr{Basic}), s.ptr, n, &result)
    result
end

function Base.convert(::Type{Vector}, x::CVecBasic)
    n = Base.length(x)
    [x[i-1] for i in 1:n]
end


text:
struct BenchmarkResult{T,I<:Union{Int,NTuple{3,Int}}}
    libraries::Vector{Symbol}
    sizes::Vector{I}
    gflops::Array{Float64,3}
    times::Array{Float64,3}
    threaded::Bool
end
function BenchmarkResult{T}(libraries, sizes, gflops, times, threaded) where {T}
    gflopsperm = permutedims(gflops, (2,3,1))
    timesperm = permutedims(times, (2,3,1))
    I = eltype(sizes)
    BenchmarkResult{T,I}(libraries, convert(Vector{I},sizes), gflopsperm, timesperm, threaded)
end

"""
    benchmark_result_type(benchmark_result::BenchmarkResult)
"""
function benchmark_result_type(::BenchmarkResult{T}) where {T}
    return T
end

function get_measure_index(measure::Symbol)::Int
    j = findfirst(==(measure), (:minimum,:median,:mean,:maximum,:hmean))
    if j === nothing
        throw(ArgumentError("`measure` argument must be one of (:minimum,:median,:mean,:maximum,:hmean), but was $(repr(measure))."))
    end
    return j
end
function _benchmark_result_df(sizes, libraries, mat, measure)
    j = get_measure_index(measure)
    df = DataFrame(Size = sizes)
    for i ∈ eachindex(libraries)
        setproperty!(df, libraries[i], mat[:,i,j])
    end
    return df
end
function _benchmark_result_df(br::BenchmarkResult, s::Symbol = :gflops, measure = :minimum)
    _benchmark_result_df(br.sizes, br.libraries, getproperty(br, s), measure)
end

"""
    benchmark_result_df(benchmark_result::BenchmarkResult, `measure` = :minimum)

`measure` refers to the BenchmarkTools summary on times. Valid options are:
`:minimum`, `:medain`, `:mean`, `:maximum`, and `:hmean`.

 -  `:minimum` would yield the maximum `GFLOPS`, and would be the usual estimate used in Julia.
 - `:hmean`, the harmonic mean of the times, is usful if you want an average GFLOPS, instead of a GFLOPS computed with the average times.
"""
function benchmark_result_df(benchmark_result::BenchmarkResult, measure = :minimum)
    df = _benchmark_result_df(benchmark_result, :times, measure)
    df = stack(df, Not(:Size), variable_name = :Library, value_name = :Seconds)
    df.GFLOPS = @. 2e-9 * matmul_length(df.Size) ./ df.Seconds
    return df
end

"""
    benchmark_result_threaded(benchmark_result::BenchmarkResult)
"""
function benchmark_result_threaded(benchmark_result::BenchmarkResult)
    return benchmark_result.threaded
end

function Base.show(io::IO, br::BenchmarkResult{T}) where {T}
    println(io, "Benchmark Result of Matrix{$T}, threaded = $(br.threaded)")
    df = _benchmark_result_df(br)
    println(io, df)
end

function maybe_sleep(x)
    x > 1e-3 && sleep(x)
end

function benchmark_fun!(
    f!::F,
    C,
    A,
    B,
    sleep_time,
    discard_first,
    reference,
    comment::String, # `comment` is a place to put the library name, the dimensions of the matrices, etc.
) where {F}
  maybe_sleep(sleep_time)
  if discard_first
    @elapsed f!(C, A, B)
  end
  t0 = @elapsed f!(C, A, B)
  if (reference !== nothing) && (!(C ≈ reference))
    msg = "C is not approximately equal to reference"
    @error(msg, comment)
    throw(ErrorException(msg))
  end
  if 2t0 < BenchmarkTools.DEFAULT_PARAMETERS.seconds
    maybe_sleep(sleep_time)
    br = @benchmark $f!($C, $A, $B)
    tmin = min(1e-9minimum(br).time, t0)
    tmedian = 1e-9median(br).time
    tmean = 1e-9mean(br).time
    tmax = 1e-9maximum(br).time # We'll exclude the first for this...
    thmean⁻¹ = 1e9mean(inv, br.times)
  else
    maybe_sleep(sleep_time)
    t1 = @elapsed f!(C, A, B)
    maybe_sleep(sleep_time)
    t2 = @elapsed f!(C, A, B)
    if (t0+t1) < 4BenchmarkTools.DEFAULT_PARAMETERS.seconds
      maybe_sleep(sleep_time)
      t3 = @elapsed f!(C, A, B)
      tmin = minimum((t0, t1, t2, t3))
      tmedian = median((t0, t1, t2, t3))
      tmean = mean((t0, t1, t2, t3))
      tmax = maximum((t0, t1, t2, t3))
      thmean⁻¹ = mean(inv, (t0, t1, t2, t3))
    else
      tmin = minimum((t0, t1, t2))
      tmedian = median((t0, t1, t2))
      tmean = mean((t0, t1, t2))
      tmax = maximum((t0, t1, t2))
      thmean⁻¹ = mean(inv, (t0, t1, t2))
    end
  end
  return tmin, tmedian, tmean, tmax, thmean⁻¹
end
_mat_size(M, N, ::typeof(adjoint)) = (N, M)
_mat_size(M, N, ::typeof(transpose)) = (N, M)
_mat_size(M, N, ::typeof(identity)) = (M, N)
function alloc_mat(_M, _N, memory::Vector{T}, off, f = identity) where {T}
    M, N = _mat_size(_M, _N, f)
    A = f(reshape(view(memory, (off+1):(off+M*N)), (M, N)))
    A, off + align(M*N, T)
end
matmul_sizes(s::Integer) = (s,s,s)
matmul_sizes(mkn::Tuple{Vararg{Integer,3}}) = mkn
matmul_length(s) = prod(matmul_sizes(s))

junk(::Type{T}) where {T <: Integer} = typemax(T) >> 1
junk(::Type{T}) where {T} = T(NaN)

struct LogSpace
    r::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}
end
Base.IteratorSize(::Type{LogSpace}) = Base.HasShape{1}()
"""
  logspace(start, stop, length)

Defines a monotonically increasing range, log spaced when possible. Useful for defining a range of sizes for benchmarks.

```julia
julia> collect(logspace(1,100,3))
3-element Vector{Int64}:
   1
  10
 100

julia> collect(logspace(1,10,3))
3-element Vector{Int64}:
  1
  3
 10

julia> collect(logspace(1,5,3))
3-element Vector{Int64}:
 1
 2
 5

julia> collect(logspace(1,3,3))
3-element Vector{Int64}:
 1
 2
 3
```
"""
logspace(start, stop, length) = LogSpace(range(log(start),log(stop), length = length))
function Base.iterate(ls::LogSpace)
    i_s = iterate(ls.r)
    i_s === nothing && return nothing
    i, _s = i_s
    v = round(Int, exp(i))
    v, (_s, v)
end
function Base.iterate(ls::LogSpace, (s,l))
    i_s = iterate(ls.r, s)
    i_s === nothing && return nothing
    i, _s = i_s
    v = max(round(Int, exp(i)), l+1)
    v, (_s, v)
end
Base.length(ls::LogSpace) = length(ls.r)
Base.size(ls::LogSpace) = (length(ls.r),)
Base.axes(ls::LogSpace) = axes(ls.r)
Base.eltype(::LogSpace) = Int
Base.convert(::Type{Vector{Int}}, l::LogSpace) = collect(l)

"""
    all_libs()
"""
function all_libs()
    libs = Symbol[
        :BLIS,
        :Gaius,
        :Octavian,
        :OpenBLAS,
        :Tullio,
        :LoopVectorization
    ]
    Sys.ARCH === :x86_64 && push!(libs, :MKL)
    return libs
end

function _integer_libs()
    libs_to_exclude = Symbol[:BLIS, :MKL, :OpenBLAS]
    return sort(unique(setdiff(all_libs(), libs_to_exclude)))
end

"""
    default_libs(T)
"""
function default_libs(::Type{T}) where {T}
    if T <: Integer
        return _integer_libs()
    else
        return all_libs()
    end
end

function luflop(m, n=m; innerflop=2)
  sum(1:min(m, n)) do k
    invflop = 1
    scaleflop = isempty(k+1:m) ? 0 : sum(k+1:m)
    updateflop = isempty(k+1:n) ? 0 : sum(k+1:n) do j
      isempty(k+1:m) ? 0 : sum(k+1:m) do i
        innerflop
      end
    end
    invflop + scaleflop + updateflop
  end * 1e-9
end
gemmflop(m,n,k) = 2e-9m*n*k

"""
    runbench(T = Float64;
             libs = default_libs(T),
             sizes = logspace(2, 4000, 200),
             threaded::Bool = Threads.nthreads() > 1,
             A_transform = identity,
             B_transform = identity,
             sleep_time = 0.0)

 - T: The element type of the matrices.
 - libs: Libraries to benchmark.
 - sizes: Sizes of matrices to benchmark. Must be an iterable with either
          `eltype(sizes) === Int` or `eltype(sizes) === NTuple{3,Int}`.
          If the former, the matrices are square, with each dimension equal to the value.
          If `i::NTuple{3,Int}`, it benchmarks `C = A * B` where `A` is `i[1]` by `i[2]`,
            `B` is `i[2]` by `i[3]` and `C` is `i[1]` by `i[3]`.
 - threaded: Should it benchmark multithreaded implementations?
 - A_transform: a function to apply to `A`. Defaults to `identity`, but can be `adjoint`.
 - B_transofrm: a function to apply to `B`. Defaults to `identity`, but can be `adjoint`.
 - sleep_time: The use of this keyword argument is discouraged. If set, it will call `sleep`
       in between benchmarks, the idea being to help keep the CPU cool. This is an unreliable
       means of trying to get more reliable benchmarks. Instead, it's reccommended you disable
       your systems turbo. Disabling it -- and reenabling when you're done benchmarking --
       should be possible without requiring a reboot.

"""
function runbench(
  ::Type{T} = Float64;
  libs = default_libs(T),
  sizes = logspace(2, 4000, 200),
  threaded::Bool = Threads.nthreads() > 1,
  A_transform = identity,
  B_transform = identity,
  sleep_time = 0.0
) where {T}
  if threaded
    Sys.ARCH === :x86_64 && mkl_set_num_threads(num_cores())
    openblas_set_num_threads(num_cores())
    blis_set_num_threads(num_cores())
  else
    Sys.ARCH === :x86_64 && mkl_set_num_threads(1)
    openblas_set_num_threads(1)
    blis_set_num_threads(1)
  end
  benchtime = BenchmarkTools.DEFAULT_PARAMETERS.seconds
  BenchmarkTools.DEFAULT_PARAMETERS.seconds = 0.5
  funcs = getfuncs(libs, threaded)
  sizevec = collect(sizes)
  # Hack to workaround https://github.com/JuliaCI/BenchmarkTools.jl/issues/127
  # Use the same memory every time, to reduce accumulation
  max_matrix_sizes = maximum(sizevec) do s
    M, K, N = matmul_sizes(s)
    align(M * K, T) + align(K * N, T) + align(M * N, T) * 2
  end
  memory = Vector{T}(undef, max_matrix_sizes)
  library = reduce(vcat, (libs for _ ∈ eachindex(sizevec)))
  times = Array{Float64}(undef, 5, length(sizes), length(libs))
  gflop = similar(times);

  discard_first = true # force when compiling

  p = Progress(length(sizes))
  gflop_report_type = NamedTuple{(:MedianGFLOPS, :MaxGFLOPS), Tuple{Float64, Float64}}
  last_perfs = Vector{Tuple{Symbol,Union{gflop_report_type,NTuple{3,Int}}}}(undef, length(libs)+1)
  for _j in 0:length(sizevec)-1
    if iseven(_j)
      j = (_j >> 1) + 1
    else
      j = length(sizevec) - (_j >> 1)
    end
    s = sizevec[j]
    M, K, N = matmul_sizes(s)
    A,  off = alloc_mat(M, K, memory,   0, A_transform)
    B,  off = alloc_mat(K, N, memory, off, B_transform)
    rand!(A); rand!(B);
    C0, off = alloc_mat(M, N, memory, off)
    C1, off = alloc_mat(M, N, memory, off)
    last_perfs[1] = (:Size, (M,K,N) .% Int)
    for i ∈ eachindex(funcs)
      C, ref = i == 1 ? (C0, nothing) : (fill!(C1,junk(T)), C0)
      lib = library[i]
      comment = "lib=$(lib), M=$(M), K=$(K), N=$(N)"
      t = benchmark_fun!(
        funcs[i],
        C,
        A,
        B,
        sleep_time,
        discard_first,
        ref,
        comment,
      )
      gffactor = gemmflop(M,K,N)
      @inbounds for k ∈ 1:4
        times[k,j,i] = t[k]
        gflop[k,j,i] = gffactor / t[k]
      end
      times[5,j,i] = inv(t[5])
      gflop[5,j,i] = gffactor * t[5]
      gflops = round.((gflop[1,j,i], gflop[2,j,i]), sigdigits = 4)
      gflops = (
        MedianGFLOPS = round(gflop[2,j,i], sigdigits = 4),
        MaxGFLOPS = round(gflop[1,j,i], sigdigits = 4)
      )
      last_perfs[i+1] = (libs[i], gflops)
    end
    ProgressMeter.next!(p; showvalues = last_perfs)
    if isodd(_j)
      discard_first = false
    end
  end
  BenchmarkTools.DEFAULT_PARAMETERS.seconds = benchtime # restore previous state
  BenchmarkResult{T}(libs, sizes, gflop, times, threaded)
end

@static if Sys.ARCH === :x86_64
const LUFUNCS = Dict(:RecursiveFactorization => RecursiveFactorization.lu!, :MKL => lumkl!, :OpenBLAS => luopenblas!)
else
const LUFUNCS = Dict(:RecursiveFactorization => RecursiveFactorization.lu!, :OpenBLAS => luopenblas)
end
struct LUWrapperFunc{F}; f::F; end
(lu::LUWrapperFunc)(A,B,C) = lu.f(copyto!(A,B))
function runlubench(
  ::Type{T} = Float64;
  libs = [:RecursiveFactorization, :MKL, :OpenBLAS],
  sizes = logspace(2, 4000, 200),
  threaded::Bool = Threads.nthreads() > 1,
  A_transform = identity,
  B_transform = identity,
  sleep_time = 0.0
) where {T}
  funcs = LUWrapperFunc.(getindex.(Ref(LUFUNCS), libs))
  if threaded
    mkl_set_num_threads(num_cores())
    openblas_set_num_threads(num_cores())
  else
    mkl_set_num_threads(1)
    openblas_set_num_threads(1)
  end
  benchtime = BenchmarkTools.DEFAULT_PARAMETERS.seconds
  BenchmarkTools.DEFAULT_PARAMETERS.seconds = 0.5
  sizevec = collect(sizes)
  # Hack to workaround https://github.com/JuliaCI/BenchmarkTools.jl/issues/127
  # Use the same memory every time, to reduce accumulation
  max_matrix_sizes = 2maximum(sizevec)^2 + (256 ÷ sizeof(T))
  memory = Vector{T}(undef, max_matrix_sizes)
  library = reduce(vcat, (libs for _ ∈ eachindex(sizevec)))
  times = Array{Float64}(undef, 5, length(sizes), length(libs))
  gflop = similar(times);

  discard_first = true # force when compiling

  p = Progress(length(sizes))
  gflop_report_type = NamedTuple{(:MedianGFLOPS, :MaxGFLOPS), Tuple{Float64, Float64}}
  last_perfs = Vector{Tuple{Symbol,Union{gflop_report_type,NTuple{2,Int}}}}(undef, length(libs)+1)
  for _j in 0:length(sizevec)-1
    if iseven(_j)
      j = (_j >> 1) + 1
    else
      j = length(sizevec) - (_j >> 1)
    end
    N = sizevec[j]
    M = N
    A,  off = alloc_mat(M, N, memory,   0, A_transform)
    rand!(A); #rand!(B);
    @inbounds for n ∈ 1:N, m ∈ 1:M
      A[m,n] = (A[m,n] + (m == n))
    end
    B, off = alloc_mat(M, N, memory, off, B_transform)
    last_perfs[1] = (:Size, (M,N) .% Int)
    for i ∈ eachindex(funcs)
      lib = library[i]
      comment = "lib=$(lib), M=$(M), N=$(N)"
      t = benchmark_fun!(
        funcs[i],
        B,
        A,
        nothing,
        sleep_time,
        discard_first,
        nothing,
        comment,
      )
      gffactor = luflop(M,N)
      @inbounds for k ∈ 1:4
        times[k,j,i] = t[k]
        gflop[k,j,i] = gffactor / t[k]
      end
      times[5,j,i] = inv(t[5])
      gflop[5,j,i] = gffactor * t[5]
      gflops = round.((gflop[1,j,i], gflop[2,j,i]), sigdigits = 4)
      gflops = (
        MedianGFLOPS = round(gflop[2,j,i], sigdigits = 4),
        MaxGFLOPS = round(gflop[1,j,i], sigdigits = 4)
      )
      last_perfs[i+1] = (libs[i], gflops)
    end
    ProgressMeter.next!(p; showvalues = last_perfs)
    if isodd(_j)
      discard_first = false
    end
  end
  BenchmarkTools.DEFAULT_PARAMETERS.seconds = benchtime # restore previous state
  BenchmarkResult{T}(libs, sizes, gflop, times, threaded)
end

ext: jl


text:
function initialize_DCside!(
    device_states,
    device::PSY.DynamicInverter{C, O, IC, PSY.FixedDCSource, P, F},
) where {
    C <: PSY.Converter,
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
}

    #Update inner_vars
    get_inner_vars(device)[Vdc_var] = PSY.get_voltage(PSY.get_dc_source(device))
end


text:
(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module SequencesDemo

using IntegerSequences

println("\nProduct 1*2*3")
p = ∏([1, 2, 3])
println(p)

println("\nList of first 10000 Kolakoski numbers (timing):")
@time KolakoskiList(10000)

println("\nList of first 10000 partition numbers (timing):")
@time PartitionNumberList(10000)

println("\nList of first 10000 Clausen numbers (timing):")
@time ClausenNumberList(10000)

println("\nList of first 10000 Ramanujan tau numbers (timing):")
@time RamanujanTauList(10000)

println("\nRamanujan tau of n = 10000:")
r = RamanujanTau(10000)
println(r)

println("\nAll partitions of 1, 2, 3, 4 and 5.")
for n in 1:5 Partition(n); println() end

println("\nNumber of points in square lattice on the
circle of radius √n.")
SeqShow(L004018(12))

println("\nThe product of the prime numbers dividing n.")
println(typeof([Radical(n) for n in 0:9]))
SeqShow([Radical(n) for n in 0:9])

println("\nOr:")
A = ZArray(9, Radical)
SeqShow(A)

println("\n=============")

println("\nThe exponential transform of Pascal's triangle.")
T = T055883(8)
ShowAsΔ(T)

println("\nThe row sums of the triangle above.")
Println(RowSums(T))

println("\nThe factorial of 100.")
f = F!(100)
println(f)

println("\nThe rising factorial ↑(20, 80).") # (a.k.a. Pochhammer)
rf = 20 ↑ 80
println(rf)

println("\nThe falling factorial ↓(80, 60).")
ff = 80 ↓ 60
println(ff)

println("\nThe number of divisors of n = 25920.")
t = τ(25920)
println(t)

println("\nThe sum of divisors of n = 25920.")
s = σ(25920)
println(s)

println("\nNumber of acyclic orientations of the
Turán graph T(2n, n) for n = 20.")
g = V033815(20)
println(g)

println("\nThe binomial function defined for all fmpz.
The same way as Maple and Mathematica compute the binomial.
First the familiar Pascal case:")
for n in 0:8
    for k in 0:n print(lpad(Binomial(n, k), 4)) end
    println()
end

println("\nAn extended region on the 2-dim lattice:")
for n in -5:5
    for k in -5:5
        print(lpad(Binomial(n, k), 5))
    end
    println()
end

println("\nKnuth, Graham and Patashnik write in
Concrete Mathematics: Hear us, O mathematicians
of the world! Let us not wait any longer! We can
make many formulas clearer by defining a new notation
now! Let us agree to write m ⊥ n, and to say m is
prime to n, if m and n are relatively prime.\n")
for n in 1:6
    for m in 1:6 print(⊥(m, n), " " ) end
    println()
end

println("\nUsing conversion to integer this can also be written as:\n")
for n in 1:6
    for m in 1:6 print(Int(⊥(m, n)), " " ) end
    println()
end

# Generating the OLMS logo.
# ⍊(n, k) means k is strong prime to n.
# n is strong prime to k iff n is prime to k and n does not divide k-1.
# Insert the line sum(L) == n - 5 && println(n) if you want to see why this
# triangle is called the Save-Prime-Triangle.

println("\nThe Save-Prime-Triangle, the logo of the OLMS
(open library of mathematical sequences).\n")
for n in 5:23
    T = [k for k in 1:n if ⍊(n, k)]
    L = [Int(k ∈ T) for k in 3:(n - 2)]
    println(L)
end

println("\n... have fun with Sequences!")

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
using HMCExamples

HMCExamples.main_FVGQ_1_kalman(ARGS)

max_stars_repo_name: HighDimensionalEconLab/HMCExamples

max_issues_repo_name: HighDimensionalEconLab/HMCExamples

max_forks_repo_name: HighDimensionalEconLab/HMCExamples


text:
using Scheduler
using Test
using Dates

@testset "Initializers" begin
    # Write your tests here.
    HAC = Qualification("HAC", Date(2019,1,1), Date(2020,6,1))
    @test HAC.Name == "HAC"
    @test HAC.Granted == Date(2019,1,1)
    @test HAC.Expires == Date(2020,6,1)
    @test HAC.CanExpire == true

    ANI = Qualification("ANI")
    @test ANI.Name == "ANI"
    @test ANI.Granted == Date(0)
    @test ANI.Expires == Date(0)
    @test ANI.CanExpire == false

    FCP = Qualification("FCP",Date(2018,1,1))
    @test FCP.Name == "FCP"
    @test FCP.Granted == Date(2018,1,1)
    @test FCP.Expires == Date(0)
    @test FCP.CanExpire == false

    crew = Aircrew("Alex",[HAC,ANI,FCP])
    @test crew.Name == "Alex"
    @test crew.Qualifications == [HAC,ANI,FCP]

    crew = Aircrew("old")
    @test crew.Name == "old"
    @test crew.Qualifications == []
    crew.Name = "hi"
    @test crew.Name == "hi"
    @test crew.Qualifications == []

    pos1 = Position("Pilot",["HAC","ANI"])
    @test pos1.RequiredQualifications == ["HAC","ANI"]

    evt = FlightEvent(
        DateTime(2020,1,1,8,0,0,0),
        DateTime(2020,1,1,10,0,0,0),
        [Position("Pilot",["HAC","FCP"]),Position("Copilot")]
    )
    @test evt.Start == DateTime(2020,1,1,8,0,0,0)
    @test evt.End == DateTime(2020,1,1,10,0,0,0)
    @test length(evt.Positions) == 2
    @test evt.Positions[1].RequiredQualifications == ["HAC","FCP"]
    @test evt.Positions[2].RequiredQualifications == []
end

@testset "isValidPersonForPosition" begin
    Alex = Aircrew("Alex",[Qualification("HAC"),Qualification("ANI"),Qualification("FCP")])
    Drew = Aircrew("Drew",[Qualification("HAC",Date(0),Date(2020,1,1))])
    John = Aircrew("John",[Qualification("HAC",Date(0),Date(2021,2,2))])
    Eric = Aircrew("Eric")
    Bill = NonAircrew("Bill")
    pilot = Position("Pilot",["HAC"])
    copilot = Position("Copilot")
    evt = FlightEvent(
        DateTime(2020,1,1,8),
        DateTime(2020,1,1,10),
        [pilot,copilot]
    )
    @test true   == Scheduler.isValidPersonForPosition(evt, pilot,   Alex)
    @test true   == Scheduler.isValidPersonForPosition(evt, copilot, Alex)

    @test false  == Scheduler.isValidPersonForPosition(evt, pilot,   Drew)
    @test true   == Scheduler.isValidPersonForPosition(evt, copilot, Drew)

    @test true   == Scheduler.isValidPersonForPosition(evt, pilot,   John)
    @test true   == Scheduler.isValidPersonForPosition(evt, copilot, John)

    @test false  == Scheduler.isValidPersonForPosition(evt, pilot,   Eric)
    @test true   == Scheduler.isValidPersonForPosition(evt, copilot, Eric)

    @test false  == Scheduler.isValidPersonForPosition(evt, pilot,   Bill)
    @test false  == Scheduler.isValidPersonForPosition(evt, copilot, Bill)

end

@testset "Assigning Names" begin
    Alex = Aircrew("Alex",[Qualification("HAC"),Qualification("ANI"),Qualification("FCP")])
    Drew = Aircrew("Drew",[Qualification("HAC",Date(0),Date(2020,1,1))])
    John = Aircrew("John",[Qualification("HAC",Date(0),Date(2021,2,2))])
    Eric = Aircrew("Eric")
    Bill = NonAircrew("Bill")
    pilot = Position("Pilot",["HAC"])
    copilot = Position("Copilot")
    evt = FlightEvent(
        DateTime(2020,1,1,8),
        DateTime(2020,1,1,10),
        [pilot,copilot]
    )
    evt = FlightEvent(
        DateTime(2020,1,1,8),
        DateTime(2020,1,1,10),
        [pilot,copilot]
    )
    people = [Alex, Drew, John, Eric, Bill]
end


text:
using Ensemblearn
using CSV
using DataFrames

export
    get_iris_data

pkg_path = abspath(joinpath(dirname(pathof(Ensemblearn)), ".."))
path_to_iris_data = "$(pkg_path)/test/data/iris.csv"

dict_mapping = Dict("Iris-setosa"=> 0, "Iris-versicolor"=> 1, "Iris-virginica"=> 2)

function get_iris_data()
    # Read test data
    df = CSV.File(path_to_iris_data) |> DataFrame

    # Filter data to have only 2 classes
    df = filter(r -> r.Name != "Iris-virginica", df)

    # Get labels as integers
    y = map(r -> dict_mapping[r], df.Name)
    x = hcat(df.PetalLength, df.PetalWidth, df.SepalLength, df.SepalWidth)

    x, y
end


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

## array.jl: Dense arrays

## Type aliases for convenience ##
"""
    AbstractVector{T}

Supertype for one-dimensional arrays (or array-like types) with
elements of type `T`. Alias for [`AbstractArray{T,1}`](@ref).
"""
const AbstractVector{T} = AbstractArray{T,1}

"""
    AbstractMatrix{T}

Supertype for two-dimensional arrays (or array-like types) with
elements of type `T`. Alias for [`AbstractArray{T,2}`](@ref).
"""
const AbstractMatrix{T} = AbstractArray{T,2}
const AbstractVecOrMat{T} = Union{AbstractVector{T}, AbstractMatrix{T}}
const RangeIndex = Union{Int, Range{Int}, AbstractUnitRange{Int}}
const DimOrInd = Union{Integer, AbstractUnitRange}
const IntOrInd = Union{Int, AbstractUnitRange}
const DimsOrInds{N} = NTuple{N,DimOrInd}
const NeedsShaping = Union{Tuple{Integer,Vararg{Integer}}, Tuple{OneTo,Vararg{OneTo}}}

"""
    Array{T,N} <: AbstractArray{T,N}

`N`-dimensional dense array with elements of type `T`.
"""
Array

"""
    Vector{T} <: AbstractVector{T}

One-dimensional dense array with elements of type `T`, often used to represent
a mathematical vector. Alias for [`Array{T,1}`](@ref).
"""
const Vector{T} = Array{T,1}

"""
    Matrix{T} <: AbstractMatrix{T}

Two-dimensional dense array with elements of type `T`, often used to represent
a mathematical matrix. Alias for [`Array{T,2}`](@ref).
"""
const Matrix{T} = Array{T,2}
const VecOrMat{T} = Union{Vector{T}, Matrix{T}}

const DenseVector{T} = DenseArray{T,1}
const DenseMatrix{T} = DenseArray{T,2}
const DenseVecOrMat{T} = Union{DenseVector{T}, DenseMatrix{T}}

## Basic functions ##

"""
    eltype(type)

Determine the type of the elements generated by iterating a collection of the given `type`.
For associative collection types, this will be a `Pair{KeyType,ValType}`. The definition
`eltype(x) = eltype(typeof(x))` is provided for convenience so that instances can be passed
instead of types. However the form that accepts a type argument should be defined for new
types.

# Examples
```jldoctest
julia> eltype(ones(Float32,2,2))
Float32

julia> eltype(ones(Int8,2,2))
Int8
```
"""
eltype(::Type) = Any
eltype(::Type{Any}) = Any
eltype(::Type{Bottom}) = throw(ArgumentError("Union{} does not have elements"))
eltype(t::DataType) = eltype(supertype(t))
eltype(x) = eltype(typeof(x))

import Core: arraysize, arrayset, arrayref

vect() = Array{Any,1}(0)
vect(X::T...) where {T} = T[ X[i] for i = 1:length(X) ]

function vect(X...)
    T = promote_typeof(X...)
    #T[ X[i] for i=1:length(X) ]
    # TODO: this is currently much faster. should figure out why. not clear.
    return copy!(Array{T,1}(length(X)), X)
end

size(a::Array, d) = arraysize(a, d)
size(a::Vector) = (arraysize(a,1),)
size(a::Matrix) = (arraysize(a,1), arraysize(a,2))
size(a::Array{<:Any,N}) where {N} = (@_inline_meta; ntuple(M -> size(a, M), Val(N)))

asize_from(a::Array, n) = n > ndims(a) ? () : (arraysize(a,n), asize_from(a, n+1)...)

length(a::Array) = arraylen(a)
elsize(a::Array{T}) where {T} = isbits(T) ? sizeof(T) : sizeof(Ptr)
sizeof(a::Array) = Core.sizeof(a)

function isassigned(a::Array, i::Int...)
    @_inline_meta
    ii = (sub2ind(size(a), i...) % UInt) - 1
    @boundscheck ii < length(a) % UInt || return false
    ccall(:jl_array_isassigned, Cint, (Any, UInt), a, ii) == 1
end

## copy ##

function unsafe_copy!(dest::Ptr{T}, src::Ptr{T}, n) where T
    # Do not use this to copy data between pointer arrays.
    # It can't be made safe no matter how carefully you checked.
    ccall(:memmove, Ptr{Void}, (Ptr{Void}, Ptr{Void}, UInt),
          dest, src, n*sizeof(T))
    return dest
end

function unsafe_copy!(dest::Array{T}, doffs, src::Array{T}, soffs, n) where T
    if isbits(T)
        unsafe_copy!(pointer(dest, doffs), pointer(src, soffs), n)
    else
        ccall(:jl_array_ptr_copy, Void, (Any, Ptr{Void}, Any, Ptr{Void}, Int),
              dest, pointer(dest, doffs), src, pointer(src, soffs), n)
    end
    return dest
end

function copy!(dest::Array{T}, doffs::Integer, src::Array{T}, soffs::Integer, n::Integer) where T
    n == 0 && return dest
    n > 0 || throw(ArgumentError(string("tried to copy n=", n, " elements, but n should be nonnegative")))
    if soffs < 1 || doffs < 1 || soffs+n-1 > length(src) || doffs+n-1 > length(dest)
        throw(BoundsError())
    end
    unsafe_copy!(dest, doffs, src, soffs, n)
end

copy!(dest::Array{T}, src::Array{T}) where {T} = copy!(dest, 1, src, 1, length(src))

copy(a::T) where {T<:Array} = ccall(:jl_array_copy, Ref{T}, (Any,), a)

function reinterpret(::Type{T}, a::Array{S,1}) where T where S
    nel = Int(div(length(a)*sizeof(S),sizeof(T)))
    # TODO: maybe check that remainder is zero?
    return reinterpret(T, a, (nel,))
end

function reinterpret(::Type{T}, a::Array{S}) where T where S
    if sizeof(S) != sizeof(T)
        throw(ArgumentError("result shape not specified"))
    end
    reinterpret(T, a, size(a))
end

function reinterpret(::Type{T}, a::Array{S}, dims::NTuple{N,Int}) where T where S where N
    function throwbits(::Type{S}, ::Type{T}, ::Type{U}) where {S,T,U}
        @_noinline_meta
        throw(ArgumentError("cannot reinterpret Array{$(S)} to ::Type{Array{$(T)}}, type $(U) is not a bits type"))
    end
    isbits(T) || throwbits(S, T, T)
    isbits(S) || throwbits(S, T, S)
    nel = div(length(a)*sizeof(S),sizeof(T))
    if prod(dims) != nel
        _throw_dmrsa(dims, nel)
    end
    ccall(:jl_reshape_array, Array{T,N}, (Any, Any, Any), Array{T,N}, a, dims)
end

# reshaping to same # of dimensions
function reshape(a::Array{T,N}, dims::NTuple{N,Int}) where T where N
    if prod(dims) != length(a)
        _throw_dmrsa(dims, length(a))
    end
    if dims == size(a)
        return a
    end
    ccall(:jl_reshape_array, Array{T,N}, (Any, Any, Any), Array{T,N}, a, dims)
end

# reshaping to different # of dimensions
function reshape(a::Array{T}, dims::NTuple{N,Int}) where T where N
    if prod(dims) != length(a)
        _throw_dmrsa(dims, length(a))
    end
    ccall(:jl_reshape_array, Array{T,N}, (Any, Any, Any), Array{T,N}, a, dims)
end

function _throw_dmrsa(dims, len)
    @_noinline_meta
    throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $len"))
end

## Constructors ##

similar(a::Array{T,1}) where {T}                    = Array{T,1}(size(a,1))
similar(a::Array{T,2}) where {T}                    = Array{T,2}(size(a,1), size(a,2))
similar(a::Array{T,1}, S::Type) where {T}           = Array{S,1}(size(a,1))
similar(a::Array{T,2}, S::Type) where {T}           = Array{S,2}(size(a,1), size(a,2))
similar(a::Array{T}, m::Int) where {T}              = Array{T,1}(m)
similar(a::Array, T::Type, dims::Dims{N}) where {N} = Array{T,N}(dims)
similar(a::Array{T}, dims::Dims{N}) where {T,N}     = Array{T,N}(dims)

# T[x...] constructs Array{T,1}
"""
    getindex(type[, elements...])

Construct a 1-d array of the specified type. This is usually called with the syntax
`Type[]`. Element values can be specified using `Type[a,b,c,...]`.

# Examples
```jldoctest
julia> Int8[1, 2, 3]
3-element Array{Int8,1}:
 1
 2
 3

julia> getindex(Int8, 1, 2, 3)
3-element Array{Int8,1}:
 1
 2
 3
```
"""
function getindex(::Type{T}, vals...) where T
    a = Array{T,1}(length(vals))
    @inbounds for i = 1:length(vals)
        a[i] = vals[i]
    end
    return a
end

getindex(::Type{T}) where {T} = (@_inline_meta; Array{T,1}(0))
getindex(::Type{T}, x) where {T} = (@_inline_meta; a = Array{T,1}(1); @inbounds a[1] = x; a)
getindex(::Type{T}, x, y) where {T} = (@_inline_meta; a = Array{T,1}(2); @inbounds (a[1] = x; a[2] = y); a)
getindex(::Type{T}, x, y, z) where {T} = (@_inline_meta; a = Array{T,1}(3); @inbounds (a[1] = x; a[2] = y; a[3] = z); a)

function getindex(::Type{Any}, @nospecialize vals...)
    a = Array{Any,1}(length(vals))
    @inbounds for i = 1:length(vals)
        a[i] = vals[i]
    end
    return a
end
getindex(::Type{Any}) = Array{Any,1}(0)

function fill!(a::Union{Array{UInt8}, Array{Int8}}, x::Integer)
    ccall(:memset, Ptr{Void}, (Ptr{Void}, Cint, Csize_t), a, x, length(a))
    return a
end

function fill!(a::Array{T}, x) where T<:Union{Integer,AbstractFloat}
    xT = convert(T, x)
    for i in eachindex(a)
        @inbounds a[i] = xT
    end
    return a
end

"""
    fill(x, dims)

Create an array filled with the value `x`. For example, `fill(1.0, (5,5))` returns a 5×5
array of floats, with each element initialized to `1.0`.

# Examples
```jldoctest
julia> fill(1.0, (5,5))
5×5 Array{Float64,2}:
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
```

If `x` is an object reference, all elements will refer to the same object. `fill(Foo(),
dims)` will return an array filled with the result of evaluating `Foo()` once.
"""
fill(v, dims::Dims)       = fill!(Array{typeof(v)}(dims), v)
fill(v, dims::Integer...) = fill!(Array{typeof(v)}(dims...), v)

for (fname, felt) in ((:zeros,:zero), (:ones,:one))
    @eval begin
        # allow signature of similar
        $fname(a::AbstractArray, T::Type, dims::Tuple) = fill!(similar(a, T, dims), $felt(T))
        $fname(a::AbstractArray, T::Type, dims...) = fill!(similar(a,T,dims...), $felt(T))
        $fname(a::AbstractArray, T::Type=eltype(a)) = fill!(similar(a,T), $felt(T))

        $fname(T::Type, dims::Tuple) = fill!(Array{T}(Dims(dims)), $felt(T))
        $fname(dims::Tuple) = ($fname)(Float64, dims)
        $fname(T::Type, dims...) = $fname(T, dims)
        $fname(dims...) = $fname(dims)
    end
end

"""
    eye([T::Type=Float64,] m::Integer, n::Integer)

`m`-by-`n` identity matrix.
The default element type is [`Float64`](@ref).

# Examples
```jldoctest
julia> eye(3, 4)
3×4 Array{Float64,2}:
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0

julia> eye(2, 2)
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0

julia> eye(Int, 2, 2)
2×2 Array{Int64,2}:
 1  0
 0  1
```
"""
function eye(::Type{T}, m::Integer, n::Integer) where T
    a = zeros(T,m,n)
    for i = 1:min(m,n)
        a[i,i] = oneunit(T)
    end
    return a
end

"""
    eye(m, n)

`m`-by-`n` identity matrix.
"""
eye(m::Integer, n::Integer) = eye(Float64, m, n)
eye(::Type{T}, n::Integer) where {T} = eye(T, n, n)
"""
    eye([T::Type=Float64,] n::Integer)

`n`-by-`n` identity matrix.
The default element type is [`Float64`](@ref).

# Examples
```jldoctest
julia> eye(Int, 2)
2×2 Array{Int64,2}:
 1  0
 0  1

julia> eye(2)
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0
```
"""
eye(n::Integer) = eye(Float64, n)

"""
    eye(A)

Constructs an identity matrix of the same dimensions and type as `A`.

# Examples
```jldoctest
julia> A = [1 2 3; 4 5 6; 7 8 9]
3×3 Array{Int64,2}:
 1  2  3
 4  5  6
 7  8  9

julia> eye(A)
3×3 Array{Int64,2}:
 1  0  0
 0  1  0
 0  0  1
```

Note the difference from [`ones`](@ref).
"""
eye(x::AbstractMatrix{T}) where {T} = eye(typeof(one(T)), size(x, 1), size(x, 2))

function _one(unit::T, x::AbstractMatrix) where T
    m,n = size(x)
    m==n || throw(DimensionMismatch("multiplicative identity defined only for square matrices"))
    eye(T, m)
end

one(x::AbstractMatrix{T}) where {T} = _one(one(T), x)
oneunit(x::AbstractMatrix{T}) where {T} = _one(oneunit(T), x)

## Conversions ##

convert(::Type{Vector}, x::AbstractVector{T}) where {T} = convert(Vector{T}, x)
convert(::Type{Matrix}, x::AbstractMatrix{T}) where {T} = convert(Matrix{T}, x)

convert(::Type{Array{T}}, x::Array{T,n}) where {T,n} = x
convert(::Type{Array{T,n}}, x::Array{T,n}) where {T,n} = x

convert(::Type{Array{T}}, x::AbstractArray{S,n}) where {T,n,S} = convert(Array{T,n}, x)
convert(::Type{Array{T,n}}, x::AbstractArray{S,n}) where {T,n,S} = copy!(Array{T,n}(size(x)), x)

promote_rule(a::Type{Array{T,n}}, b::Type{Array{S,n}}) where {T,n,S} = el_same(promote_type(T,S), a, b)

## copying iterators to containers

"""
    collect(element_type, collection)

Return an `Array` with the given element type of all items in a collection or iterable.
The result has the same shape and number of dimensions as `collection`.

# Examples
```jldoctest
julia> collect(Float64, 1:2:5)
3-element Array{Float64,1}:
 1.0
 3.0
 5.0
```
"""
collect(::Type{T}, itr) where {T} = _collect(T, itr, iteratorsize(itr))

_collect(::Type{T}, itr, isz::HasLength) where {T} = copy!(Array{T,1}(Int(length(itr)::Integer)), itr)
_collect(::Type{T}, itr, isz::HasShape) where {T}  = copy!(similar(Array{T}, indices(itr)), itr)
function _collect(::Type{T}, itr, isz::SizeUnknown) where T
    a = Array{T,1}(0)
    for x in itr
        push!(a,x)
    end
    return a
end

# make a collection similar to `c` and appropriate for collecting `itr`
_similar_for(c::AbstractArray, T, itr, ::SizeUnknown) = similar(c, T, 0)
_similar_for(c::AbstractArray, T, itr, ::HasLength) = similar(c, T, Int(length(itr)::Integer))
_similar_for(c::AbstractArray, T, itr, ::HasShape) = similar(c, T, indices(itr))
_similar_for(c, T, itr, isz) = similar(c, T)

"""
    collect(collection)

Return an `Array` of all items in a collection or iterator. For associative collections, returns
`Pair{KeyType, ValType}`. If the argument is array-like or is an iterator with the `HasShape()`
trait, the result will have the same shape and number of dimensions as the argument.

# Examples
```jldoctest
julia> collect(1:2:13)
7-element Array{Int64,1}:
  1
  3
  5
  7
  9
 11
 13
```
"""
collect(itr) = _collect(1:1 #= Array =#, itr, iteratoreltype(itr), iteratorsize(itr))

collect(A::AbstractArray) = _collect_indices(indices(A), A)

collect_similar(cont, itr) = _collect(cont, itr, iteratoreltype(itr), iteratorsize(itr))

_collect(cont, itr, ::HasEltype, isz::Union{HasLength,HasShape}) =
    copy!(_similar_for(cont, eltype(itr), itr, isz), itr)

function _collect(cont, itr, ::HasEltype, isz::SizeUnknown)
    a = _similar_for(cont, eltype(itr), itr, isz)
    for x in itr
        push!(a,x)
    end
    return a
end

_collect_indices(::Tuple{}, A) = copy!(Vector{eltype(A)}(), A)
_collect_indices(indsA::Tuple{Vararg{OneTo}}, A) =
    copy!(Array{eltype(A)}(length.(indsA)), A)
function _collect_indices(indsA, A)
    B = Array{eltype(A)}(length.(indsA))
    copy!(B, CartesianRange(indices(B)), A, CartesianRange(indsA))
end

if isdefined(Core, :Inference)
    _default_eltype(@nospecialize itrt) = Core.Inference.return_type(first, Tuple{itrt})
else
    _default_eltype(@nospecialize itr) = Any
end

_array_for(::Type{T}, itr, ::HasLength) where {T} = Array{T,1}(Int(length(itr)::Integer))
_array_for(::Type{T}, itr, ::HasShape) where {T} = similar(Array{T}, indices(itr))

function collect(itr::Generator)
    isz = iteratorsize(itr.iter)
    et = _default_eltype(typeof(itr))
    if isa(isz, SizeUnknown)
        return grow_to!(Array{et,1}(0), itr)
    else
        st = start(itr)
        if done(itr,st)
            return _array_for(et, itr.iter, isz)
        end
        v1, st = next(itr, st)
        collect_to_with_first!(_array_for(typeof(v1), itr.iter, isz), v1, itr, st)
    end
end

_collect(c, itr, ::EltypeUnknown, isz::SizeUnknown) =
    grow_to!(_similar_for(c, _default_eltype(typeof(itr)), itr, isz), itr)

function _collect(c, itr, ::EltypeUnknown, isz::Union{HasLength,HasShape})
    st = start(itr)
    if done(itr,st)
        return _similar_for(c, _default_eltype(typeof(itr)), itr, isz)
    end
    v1, st = next(itr, st)
    collect_to_with_first!(_similar_for(c, typeof(v1), itr, isz), v1, itr, st)
end

function collect_to_with_first!(dest::AbstractArray, v1, itr, st)
    i1 = first(linearindices(dest))
    dest[i1] = v1
    return collect_to!(dest, itr, i1+1, st)
end

function collect_to_with_first!(dest, v1, itr, st)
    push!(dest, v1)
    return grow_to!(dest, itr, st)
end

function collect_to!(dest::AbstractArray{T}, itr, offs, st) where T
    # collect to dest array, checking the type of each result. if a result does not
    # match, widen the result type and re-dispatch.
    i = offs
    while !done(itr, st)
        el, st = next(itr, st)
        S = typeof(el)
        if S === T || S <: T
            @inbounds dest[i] = el::T
            i += 1
        else
            R = typejoin(T, S)
            new = similar(dest, R)
            copy!(new,1, dest,1, i-1)
            @inbounds new[i] = el
            return collect_to!(new, itr, i+1, st)
        end
    end
    return dest
end

function grow_to!(dest, itr)
    out = grow_to!(similar(dest,Union{}), itr, start(itr))
    return isempty(out) ? dest : out
end

function grow_to!(dest, itr, st)
    T = eltype(dest)
    while !done(itr, st)
        el, st = next(itr, st)
        S = typeof(el)
        if S === T || S <: T
            push!(dest, el::T)
        else
            new = similar(dest, typejoin(T, S))
            copy!(new, dest)
            push!(new, el)
            return grow_to!(new, itr, st)
        end
    end
    return dest
end

## Iteration ##
start(A::Array) = 1
next(a::Array,i) = (@_propagate_inbounds_meta; (a[i],i+1))
done(a::Array,i) = (@_inline_meta; i == length(a)+1)

## Indexing: getindex ##

# This is more complicated than it needs to be in order to get Win64 through bootstrap
getindex(A::Array, i1::Int) = arrayref(A, i1)
getindex(A::Array, i1::Int, i2::Int, I::Int...) = (@_inline_meta; arrayref(A, i1, i2, I...)) # TODO: REMOVE FOR #14770

# Faster contiguous indexing using copy! for UnitRange and Colon
function getindex(A::Array, I::UnitRange{Int})
    @_inline_meta
    @boundscheck checkbounds(A, I)
    lI = length(I)
    X = similar(A, lI)
    if lI > 0
        unsafe_copy!(X, 1, A, first(I), lI)
    end
    return X
end
function getindex(A::Array, c::Colon)
    lI = length(A)
    X = similar(A, lI)
    if lI > 0
        unsafe_copy!(X, 1, A, 1, lI)
    end
    return X
end

# This is redundant with the abstract fallbacks, but needed for bootstrap
function getindex(A::Array{S}, I::Range{Int}) where S
    return S[ A[i] for i in I ]
end

## Indexing: setindex! ##
setindex!(A::Array{T}, x, i1::Int) where {T} = arrayset(A, convert(T,x)::T, i1)
setindex!(A::Array{T}, x, i1::Int, i2::Int, I::Int...) where {T} = (@_inline_meta; arrayset(A, convert(T,x)::T, i1, i2, I...)) # TODO: REMOVE FOR #14770

# These are redundant with the abstract fallbacks but needed for bootstrap
function setindex!(A::Array, x, I::AbstractVector{Int})
    @_propagate_inbounds_meta
    A === I && (I = copy(I))
    for i in I
        A[i] = x
    end
    return A
end
function setindex!(A::Array, X::AbstractArray, I::AbstractVector{Int})
    @_propagate_inbounds_meta
    @boundscheck setindex_shape_check(X, length(I))
    count = 1
    if X === A
        X = copy(X)
        I===A && (I = X::typeof(I))
    elseif I === A
        I = copy(I)
    end
    for i in I
        @inbounds x = X[count]
        A[i] = x
        count += 1
    end
    return A
end

# Faster contiguous setindex! with copy!
function setindex!(A::Array{T}, X::Array{T}, I::UnitRange{Int}) where T
    @_inline_meta
    @boundscheck checkbounds(A, I)
    lI = length(I)
    @boundscheck setindex_shape_check(X, lI)
    if lI > 0
        unsafe_copy!(A, first(I), X, 1, lI)
    end
    return A
end
function setindex!(A::Array{T}, X::Array{T}, c::Colon) where T
    @_inline_meta
    lI = length(A)
    @boundscheck setindex_shape_check(X, lI)
    if lI > 0
        unsafe_copy!(A, 1, X, 1, lI)
    end
    return A
end

setindex!(A::Array, x::Number, ::Colon) = fill!(A, x)
setindex!(A::Array{T, N}, x::Number, ::Vararg{Colon, N}) where {T, N} = fill!(A, x)

# efficiently grow an array

_growbeg!(a::Vector, delta::Integer) =
    ccall(:jl_array_grow_beg, Void, (Any, UInt), a, delta)
_growend!(a::Vector, delta::Integer) =
    ccall(:jl_array_grow_end, Void, (Any, UInt), a, delta)
_growat!(a::Vector, i::Integer, delta::Integer) =
    ccall(:jl_array_grow_at, Void, (Any, Int, UInt), a, i - 1, delta)

# efficiently delete part of an array

_deletebeg!(a::Vector, delta::Integer) =
    ccall(:jl_array_del_beg, Void, (Any, UInt), a, delta)
_deleteend!(a::Vector, delta::Integer) =
    ccall(:jl_array_del_end, Void, (Any, UInt), a, delta)
_deleteat!(a::Vector, i::Integer, delta::Integer) =
    ccall(:jl_array_del_at, Void, (Any, Int, UInt), a, i - 1, delta)

## Dequeue functionality ##

function push!(a::Array{T,1}, item) where T
    # convert first so we don't grow the array if the assignment won't work
    itemT = convert(T, item)
    _growend!(a, 1)
    a[end] = itemT
    return a
end

function push!(a::Array{Any,1}, @nospecialize item)
    _growend!(a, 1)
    arrayset(a, item, length(a))
    return a
end

function append!(a::Array{<:Any,1}, items::AbstractVector)
    itemindices = eachindex(items)
    n = length(itemindices)
    _growend!(a, n)
    copy!(a, length(a)-n+1, items, first(itemindices), n)
    return a
end

append!(a::Vector, iter) = _append!(a, iteratorsize(iter), iter)
push!(a::Vector, iter...) = append!(a, iter)

function _append!(a, ::Union{HasLength,HasShape}, iter)
    n = length(a)
    resize!(a, n+length(iter))
    @inbounds for (i,item) in zip(n+1:length(a), iter)
        a[i] = item
    end
    a
end

function _append!(a, ::IteratorSize, iter)
    for item in iter
        push!(a, item)
    end
    a
end

"""
    prepend!(a::Vector, items) -> collection

Insert the elements of `items` to the beginning of `a`.

# Examples
```jldoctest
julia> prepend!([3],[1,2])
3-element Array{Int64,1}:
 1
 2
 3
```
"""
function prepend! end

function prepend!(a::Array{<:Any,1}, items::AbstractVector)
    itemindices = eachindex(items)
    n = length(itemindices)
    _growbeg!(a, n)
    if a === items
        copy!(a, 1, items, n+1, n)
    else
        copy!(a, 1, items, first(itemindices), n)
    end
    return a
end

prepend!(a::Vector, iter) = _prepend!(a, iteratorsize(iter), iter)
unshift!(a::Vector, iter...) = prepend!(a, iter)

function _prepend!(a, ::Union{HasLength,HasShape}, iter)
    n = length(iter)
    _growbeg!(a, n)
    i = 0
    for item in iter
        @inbounds a[i += 1] = item
    end
    a
end
function _prepend!(a, ::IteratorSize, iter)
    n = 0
    for item in iter
        n += 1
        unshift!(a, item)
    end
    reverse!(a, 1, n)
    a
end

"""
    resize!(a::Vector, n::Integer) -> Vector

Resize `a` to contain `n` elements. If `n` is smaller than the current collection
length, the first `n` elements will be retained. If `n` is larger, the new elements are not
guaranteed to be initialized.

# Examples
```jldoctest
julia> resize!([6, 5, 4, 3, 2, 1], 3)
3-element Array{Int64,1}:
 6
 5
 4

julia> a = resize!([6, 5, 4, 3, 2, 1], 8);

julia> length(a)
8

julia> a[1:6]
6-element Array{Int64,1}:
 6
 5
 4
 3
 2
 1
```
"""
function resize!(a::Vector, nl::Integer)
    l = length(a)
    if nl > l
        ccall(:jl_array_grow_end, Void, (Any, UInt), a, nl-l)
    else
        if nl < 0
            throw(ArgumentError("new length must be ≥ 0"))
        end
        _deleteend!(a, l-nl)
    end
    return a
end

function sizehint!(a::Vector, sz::Integer)
    ccall(:jl_array_sizehint, Void, (Any, UInt), a, sz)
    a
end

function pop!(a::Vector)
    if isempty(a)
        throw(ArgumentError("array must be non-empty"))
    end
    item = a[end]
    _deleteend!(a, 1)
    return item
end

"""
    unshift!(collection, items...) -> collection

Insert one or more `items` at the beginning of `collection`.

# Examples
```jldoctest
julia> unshift!([1, 2, 3, 4], 5, 6)
6-element Array{Int64,1}:
 5
 6
 1
 2
 3
 4
```
"""
function unshift!(a::Array{T,1}, item) where T
    item = convert(T, item)
    _growbeg!(a, 1)
    a[1] = item
    return a
end

function shift!(a::Vector)
    if isempty(a)
        throw(ArgumentError("array must be non-empty"))
    end
    item = a[1]
    _deletebeg!(a, 1)
    return item
end

"""
    insert!(a::Vector, index::Integer, item)

Insert an `item` into `a` at the given `index`. `index` is the index of `item` in
the resulting `a`.

# Examples
```jldoctest
julia> insert!([6, 5, 4, 2, 1], 4, 3)
6-element Array{Int64,1}:
 6
 5
 4
 3
 2
 1
```
"""
function insert!(a::Array{T,1}, i::Integer, item) where T
    # Throw convert error before changing the shape of the array
    _item = convert(T, item)
    _growat!(a, i, 1)
    # _growat! already did bound check
    @inbounds a[i] = _item
    return a
end

"""
    deleteat!(a::Vector, i::Integer)

Remove the item at the given `i` and return the modified `a`. Subsequent items
are shifted to fill the resulting gap.

# Examples
```jldoctest
julia> deleteat!([6, 5, 4, 3, 2, 1], 2)
5-element Array{Int64,1}:
 6
 4
 3
 2
 1
```
"""
deleteat!(a::Vector, i::Integer) = (_deleteat!(a, i, 1); a)

function deleteat!(a::Vector, r::UnitRange{<:Integer})
    n = length(a)
    isempty(r) || _deleteat!(a, first(r), length(r))
    return a
end

"""
    deleteat!(a::Vector, inds)

Remove the items at the indices given by `inds`, and return the modified `a`.
Subsequent items are shifted to fill the resulting gap.

`inds` can be either an iterator or a collection of sorted and unique integer indices,
or a boolean vector of the same length as `a` with `true` indicating entries to delete.

# Examples
```jldoctest
julia> deleteat!([6, 5, 4, 3, 2, 1], 1:2:5)
3-element Array{Int64,1}:
 5
 3
 1

julia> deleteat!([6, 5, 4, 3, 2, 1], [true, false, true, false, true, false])
3-element Array{Int64,1}:
 5
 3
 1

julia> deleteat!([6, 5, 4, 3, 2, 1], (2, 2))
ERROR: ArgumentError: indices must be unique and sorted
Stacktrace:
 [1] _deleteat!(::Array{Int64,1}, ::Tuple{Int64,Int64}) at ./array.jl:926
 [2] deleteat!(::Array{Int64,1}, ::Tuple{Int64,Int64}) at ./array.jl:913
```
"""
deleteat!(a::Vector, inds) = _deleteat!(a, inds)
deleteat!(a::Vector, inds::AbstractVector) = _deleteat!(a, to_indices(a, (inds,))[1])

function _deleteat!(a::Vector, inds)
    n = length(a)
    s = start(inds)
    done(inds, s) && return a
    (p, s) = next(inds, s)
    q = p+1
    while !done(inds, s)
        (i,s) = next(inds, s)
        if !(q <= i <= n)
            if i < q
                throw(ArgumentError("indices must be unique and sorted"))
            else
                throw(BoundsError())
            end
        end
        while q < i
            @inbounds a[p] = a[q]
            p += 1; q += 1
        end
        q = i+1
    end
    while q <= n
        @inbounds a[p] = a[q]
        p += 1; q += 1
    end
    _deleteend!(a, n-p+1)
    return a
end

# Simpler and more efficient version for logical indexing
function deleteat!(a::Vector, inds::AbstractVector{Bool})
    n = length(a)
    length(inds) == n || throw(BoundsError(a, inds))
    p = 1
    for (q, i) in enumerate(inds)
        @inbounds a[p] = a[q]
        p += !i
    end
    _deleteend!(a, n-p+1)
    return a
end

const _default_splice = []

"""
    splice!(a::Vector, index::Integer, [replacement]) -> item

Remove the item at the given index, and return the removed item.
Subsequent items are shifted left to fill the resulting gap.
If specified, replacement values from an ordered
collection will be spliced in place of the removed item.

# Examples
```jldoctest splice!
julia> A = [6, 5, 4, 3, 2, 1]; splice!(A, 5)
2

julia> A
5-element Array{Int64,1}:
 6
 5
 4
 3
 1

julia> splice!(A, 5, -1)
1

julia> A
5-element Array{Int64,1}:
  6
  5
  4
  3
 -1

julia> splice!(A, 1, [-1, -2, -3])
6

julia> A
7-element Array{Int64,1}:
 -1
 -2
 -3
  5
  4
  3
 -1
```

To insert `replacement` before an index `n` without removing any items, use
`splice!(collection, n:n-1, replacement)`.
"""
function splice!(a::Vector, i::Integer, ins=_default_splice)
    v = a[i]
    m = length(ins)
    if m == 0
        _deleteat!(a, i, 1)
    elseif m == 1
        a[i] = ins[1]
    else
        _growat!(a, i, m-1)
        k = 1
        for x in ins
            a[i+k-1] = x
            k += 1
        end
    end
    return v
end

"""
    splice!(a::Vector, range, [replacement]) -> items

Remove items in the specified index range, and return a collection containing
the removed items.
Subsequent items are shifted left to fill the resulting gap.
If specified, replacement values from an ordered collection will be spliced in
place of the removed items.

To insert `replacement` before an index `n` without removing any items, use
`splice!(collection, n:n-1, replacement)`.

# Examples
```jldoctest splice!
julia> splice!(A, 4:3, 2)
0-element Array{Int64,1}

julia> A
8-element Array{Int64,1}:
 -1
 -2
 -3
  2
  5
  4
  3
 -1
```
"""
function splice!(a::Vector, r::UnitRange{<:Integer}, ins=_default_splice)
    v = a[r]
    m = length(ins)
    if m == 0
        deleteat!(a, r)
        return v
    end

    n = length(a)
    f = first(r)
    l = last(r)
    d = length(r)

    if m < d
        delta = d - m
        _deleteat!(a, (f - 1 < n - l) ? f : (l - delta + 1), delta)
    elseif m > d
        _growat!(a, (f - 1 < n - l) ? f : (l + 1), m - d)
    end

    k = 1
    for x in ins
        a[f+k-1] = x
        k += 1
    end
    return v
end

function empty!(a::Vector)
    _deleteend!(a, length(a))
    return a
end

# use memcmp for lexcmp on byte arrays
function lexcmp(a::Array{UInt8,1}, b::Array{UInt8,1})
    c = ccall(:memcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}, UInt),
              a, b, min(length(a),length(b)))
    return c < 0 ? -1 : c > 0 ? +1 : cmp(length(a),length(b))
end

const BitIntegerArray{N} = Union{map(T->Array{T,N}, BitInteger_types)...} where N
# use memcmp for == on bit integer types
function ==(a::Arr, b::Arr) where Arr <: BitIntegerArray
    size(a) == size(b) && 0 == ccall(
        :memcmp, Int32, (Ptr{Void}, Ptr{Void}, UInt), a, b, sizeof(eltype(Arr)) * length(a))
end

# this is ~20% faster than the generic implementation above for very small arrays
function ==(a::Arr, b::Arr) where Arr <: BitIntegerArray{1}
    len = length(a)
    len == length(b) && 0 == ccall(
        :memcmp, Int32, (Ptr{Void}, Ptr{Void}, UInt), a, b, sizeof(eltype(Arr)) * len)
end

function reverse(A::AbstractVector, s=first(linearindices(A)), n=last(linearindices(A)))
    B = similar(A)
    for i = first(linearindices(A)):s-1
        B[i] = A[i]
    end
    for i = s:n
        B[i] = A[n+s-i]
    end
    for i = n+1:last(linearindices(A))
        B[i] = A[i]
    end
    return B
end
function reverseind(a::AbstractVector, i::Integer)
    li = linearindices(a)
    first(li) + last(li) - i
end

function reverse!(v::AbstractVector, s=first(linearindices(v)), n=last(linearindices(v)))
    liv = linearindices(v)
    if n <= s  # empty case; ok
    elseif !(first(liv) ≤ s ≤ last(liv))
        throw(BoundsError(v, s))
    elseif !(first(liv) ≤ n ≤ last(liv))
        throw(BoundsError(v, n))
    end
    r = n
    @inbounds for i in s:div(s+n-1, 2)
        v[i], v[r] = v[r], v[i]
        r -= 1
    end
    return v
end

# concatenations of homogeneous combinations of vectors, horizontal and vertical

vcat() = Array{Any,1}(0)
hcat() = Array{Any,1}(0)

function hcat(V::Vector{T}...) where T
    height = length(V[1])
    for j = 2:length(V)
        if length(V[j]) != height
            throw(DimensionMismatch("vectors must have same lengths"))
        end
    end
    return [ V[j][i]::T for i=1:length(V[1]), j=1:length(V) ]
end

function vcat(arrays::Vector{T}...) where T
    n = 0
    for a in arrays
        n += length(a)
    end
    arr = Array{T,1}(n)
    ptr = pointer(arr)
    if isbits(T)
        elsz = Core.sizeof(T)
    else
        elsz = Core.sizeof(Ptr{Void})
    end
    for a in arrays
        na = length(a)
        nba = na * elsz
        if isbits(T)
            ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, UInt),
                  ptr, a, nba)
        else
            ccall(:jl_array_ptr_copy, Void, (Any, Ptr{Void}, Any, Ptr{Void}, Int),
                  arr, ptr, a, pointer(a), na)
        end
        ptr += nba
    end
    return arr
end

cat(n::Integer, x::Integer...) = reshape([x...], (ntuple(x->1, n-1)..., length(x)))

## find ##

"""
    findnext(A, i::Integer)

Find the next linear index >= `i` of a non-zero element of `A`, or `0` if not found.

# Examples
```jldoctest
julia> A = [0 0; 1 0]
2×2 Array{Int64,2}:
 0  0
 1  0

julia> findnext(A,1)
2

julia> findnext(A,3)
0
```
"""
function findnext(A, start::Integer)
    for i = start:length(A)
        if A[i] != 0
            return i
        end
    end
    return 0
end

"""
    findfirst(A)

Return the linear index of the first non-zero value in `A` (determined by `A[i]!=0`).
Returns `0` if no such value is found.

# Examples
```jldoctest
julia> A = [0 0; 1 0]
2×2 Array{Int64,2}:
 0  0
 1  0

julia> findfirst(A)
2

julia> findfirst(zeros(3))
0
```
"""
findfirst(A) = findnext(A, 1)

"""
    findnext(A, v, i::Integer)

Find the next linear index >= `i` of an element of `A` equal to `v` (using `==`), or `0` if not found.

# Examples
```jldoctest
julia> A = [1 4; 2 2]
2×2 Array{Int64,2}:
 1  4
 2  2

julia> findnext(A,4,4)
0

julia> findnext(A,4,3)
3
```
"""
function findnext(A, v, start::Integer)
    for i = start:length(A)
        if A[i] == v
            return i
        end
    end
    return 0
end
"""
    findfirst(A, v)

Return the linear index of the first element equal to `v` in `A`.
Returns `0` if `v` is not found.

# Examples
```jldoctest
julia> A = [4 6; 2 2]
2×2 Array{Int64,2}:
 4  6
 2  2

julia> findfirst(A,2)
2

julia> findfirst(A,3)
0
```
"""
findfirst(A, v) = findnext(A, v, 1)

"""
    findnext(predicate::Function, A, i::Integer)

Find the next linear index >= `i` of an element of `A` for which `predicate` returns `true`, or `0` if not found.

# Examples
```jldoctest
julia> A = [1 4; 2 2]
2×2 Array{Int64,2}:
 1  4
 2  2

julia> findnext(isodd, A, 1)
1

julia> findnext(isodd, A, 2)
0
```
"""
function findnext(testf::Function, A, start::Integer)
    for i = start:length(A)
        if testf(A[i])
            return i
        end
    end
    return 0
end

"""
    findfirst(predicate::Function, A)

Return the linear index of the first element of `A` for which `predicate` returns `true`.
Returns `0` if there is no such element.

# Examples
```jldoctest
julia> A = [1 4; 2 2]
2×2 Array{Int64,2}:
 1  4
 2  2

julia> findfirst(iseven, A)
2

julia> findfirst(x -> x>10, A)
0
```
"""
findfirst(testf::Function, A) = findnext(testf, A, 1)

"""
    findprev(A, i::Integer)

Find the previous linear index <= `i` of a non-zero element of `A`, or `0` if not found.

# Examples
```jldoctest
julia> A = [0 0; 1 2]
2×2 Array{Int64,2}:
 0  0
 1  2

julia> findprev(A,2)
2

julia> findprev(A,1)
0
```
"""
function findprev(A, start::Integer)
    for i = start:-1:1
        A[i] != 0 && return i
    end
    return 0
end

"""
    findlast(A)

Return the linear index of the last non-zero value in `A` (determined by `A[i]!=0`).
Returns `0` if there is no non-zero value in `A`.

# Examples
```jldoctest
julia> A = [1 0; 1 0]
2×2 Array{Int64,2}:
 1  0
 1  0

julia> findlast(A)
2

julia> A = zeros(2,2)
2×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0

julia> findlast(A)
0
```
"""
findlast(A) = findprev(A, length(A))

"""
    findprev(A, v, i::Integer)

Find the previous linear index <= `i` of an element of `A` equal to `v` (using `==`), or `0` if not found.

# Examples
```jldoctest
julia> A = [0 0; 1 2]
2×2 Array{Int64,2}:
 0  0
 1  2

julia> findprev(A, 1, 4)
2

julia> findprev(A, 1, 1)
0
```
"""
function findprev(A, v, start::Integer)
    for i = start:-1:1
        A[i] == v && return i
    end
    return 0
end

"""
    findlast(A, v)

Return the linear index of the last element equal to `v` in `A`.
Returns `0` if there is no element of `A` equal to `v`.

# Examples
```jldoctest
julia> A = [1 2; 2 1]
2×2 Array{Int64,2}:
 1  2
 2  1

julia> findlast(A,1)
4

julia> findlast(A,2)
3

julia> findlast(A,3)
0
```
"""
findlast(A, v) = findprev(A, v, length(A))

"""
    findprev(predicate::Function, A, i::Integer)

Find the previous linear index <= `i` of an element of `A` for which `predicate` returns `true`, or
`0` if not found.

# Examples
```jldoctest
julia> A = [4 6; 1 2]
2×2 Array{Int64,2}:
 4  6
 1  2

julia> findprev(isodd, A, 1)
0

julia> findprev(isodd, A, 3)
2
```
"""
function findprev(testf::Function, A, start::Integer)
    for i = start:-1:1
        testf(A[i]) && return i
    end
    return 0
end

"""
    findlast(predicate::Function, A)

Return the linear index of the last element of `A` for which `predicate` returns `true`.
Returns `0` if there is no such element.

# Examples
```jldoctest
julia> A = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> findlast(isodd, A)
2

julia> findlast(x -> x > 5, A)
0
```
"""
findlast(testf::Function, A) = findprev(testf, A, length(A))

"""
    find(f::Function, A)

Return a vector `I` of the linear indexes of `A` where `f(A[I])` returns `true`.
If there are no such elements of `A`, find returns an empty array.

# Examples
```jldoctest
julia> A = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> find(isodd,A)
2-element Array{Int64,1}:
 1
 2

julia> find(isodd, [2, 4])
0-element Array{Int64,1}
```
"""
function find(testf::Function, A)
    # use a dynamic-length array to store the indexes, then copy to a non-padded
    # array for the return
    tmpI = Array{Int,1}(0)
    inds = _index_remapper(A)
    for (i,a) = enumerate(A)
        if testf(a)
            push!(tmpI, inds[i])
        end
    end
    I = Array{Int,1}(length(tmpI))
    copy!(I, tmpI)
    return I
end
_index_remapper(A::AbstractArray) = linearindices(A)
_index_remapper(iter) = OneTo(typemax(Int))  # safe for objects that don't implement length

"""
    find(A)

Return a vector of the linear indexes of the non-zeros in `A` (determined by `A[i]!=0`). A
common use of this is to convert a boolean array to an array of indexes of the `true`
elements. If there are no non-zero elements of `A`, `find` returns an empty array.

# Examples
```jldoctest
julia> A = [true false; false true]
2×2 Array{Bool,2}:
  true  false
 false   true

julia> find(A)
2-element Array{Int64,1}:
 1
 4

julia> find(zeros(3))
0-element Array{Int64,1}
```
"""
function find(A)
    nnzA = countnz(A)
    I = Vector{Int}(nnzA)
    count = 1
    inds = _index_remapper(A)
    for (i,a) in enumerate(A)
        if a != 0
            I[count] = inds[i]
            count += 1
        end
    end
    return I
end

find(x::Number) = x == 0 ? Array{Int,1}(0) : [1]
find(testf::Function, x::Number) = !testf(x) ? Array{Int,1}(0) : [1]

findn(A::AbstractVector) = find(A)

"""
    findn(A)

Return a vector of indexes for each dimension giving the locations of the non-zeros in `A`
(determined by `A[i]!=0`).
If there are no non-zero elements of `A`, `findn` returns a 2-tuple of empty arrays.

# Examples
```jldoctest
julia> A = [1 2 0; 0 0 3; 0 4 0]
3×3 Array{Int64,2}:
 1  2  0
 0  0  3
 0  4  0

julia> findn(A)
([1, 1, 3, 2], [1, 2, 2, 3])

julia> A = zeros(2,2)
2×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0

julia> findn(A)
(Int64[], Int64[])
```
"""
function findn(A::AbstractMatrix)
    nnzA = countnz(A)
    I = similar(A, Int, nnzA)
    J = similar(A, Int, nnzA)
    count = 1
    for j=indices(A,2), i=indices(A,1)
        if A[i,j] != 0
            I[count] = i
            J[count] = j
            count += 1
        end
    end
    return (I, J)
end

"""
    findnz(A)

Return a tuple `(I, J, V)` where `I` and `J` are the row and column indexes of the non-zero
values in matrix `A`, and `V` is a vector of the non-zero values.

# Examples
```jldoctest
julia> A = [1 2 0; 0 0 3; 0 4 0]
3×3 Array{Int64,2}:
 1  2  0
 0  0  3
 0  4  0

julia> findnz(A)
([1, 1, 3, 2], [1, 2, 2, 3], [1, 2, 4, 3])
```
"""
function findnz(A::AbstractMatrix{T}) where T
    nnzA = countnz(A)
    I = zeros(Int, nnzA)
    J = zeros(Int, nnzA)
    NZs = Array{T,1}(nnzA)
    count = 1
    if nnzA > 0
        for j=indices(A,2), i=indices(A,1)
            Aij = A[i,j]
            if Aij != 0
                I[count] = i
                J[count] = j
                NZs[count] = Aij
                count += 1
            end
        end
    end
    return (I, J, NZs)
end

"""
    findmax(itr) -> (x, index)

Returns the maximum element of the collection `itr` and its index. If there are multiple
maximal elements, then the first one will be returned. `NaN` values are ignored, unless
all elements are `NaN`. Other than the treatment of `NaN`, the result is in line with `max`.

The collection must not be empty.

# Examples
```jldoctest
julia> findmax([8,0.1,-9,pi])
(8.0, 1)

julia> findmax([1,7,7,6])
(7, 2)

julia> findmax([1,7,7,NaN])
(7.0, 2)
```
"""
function findmax(a)
    if isempty(a)
        throw(ArgumentError("collection must be non-empty"))
    end
    s = start(a)
    mi = i = 1
    m, s = next(a, s)
    while !done(a, s)
        ai, s = next(a, s)
        i += 1
        ai != ai && continue # assume x != x => x is a NaN
        if m != m || isless(m, ai)
            m = ai
            mi = i
        end
    end
    return (m, mi)
end

"""
    findmin(itr) -> (x, index)

Returns the minimum element of the collection `itr` and its index. If there are multiple
minimal elements, then the first one will be returned. `NaN` values are ignored, unless
all elements are `NaN`. Other than the treatment of `NaN`, the result is in line with `min`.

The collection must not be empty.

# Examples
```jldoctest
julia> findmin([8,0.1,-9,pi])
(-9.0, 3)

julia> findmin([7,1,1,6])
(1, 2)

julia> findmin([7,1,1,NaN])
(1.0, 2)
```
"""
function findmin(a)
    if isempty(a)
        throw(ArgumentError("collection must be non-empty"))
    end
    s = start(a)
    mi = i = 1
    m, s = next(a, s)
    while !done(a, s)
        ai, s = next(a, s)
        i += 1
        ai != ai && continue
        if m != m || isless(ai, m)
            m = ai
            mi = i
        end
    end
    return (m, mi)
end

"""
    indmax(itr) -> Integer

Returns the index of the maximum element in a collection. If there are multiple maximal
elements, then the first one will be returned. `NaN` values are ignored, unless all
elements are `NaN`.

The collection must not be empty.

# Examples
```jldoctest
julia> indmax([8,0.1,-9,pi])
1

julia> indmax([1,7,7,6])
2

julia> indmax([1,7,7,NaN])
2
```
"""
indmax(a) = findmax(a)[2]

"""
    indmin(itr) -> Integer

Returns the index of the minimum element in a collection. If there are multiple minimal
elements, then the first one will be returned. `NaN` values are ignored, unless all
elements are `NaN`.

The collection must not be empty.

# Examples
```jldoctest
julia> indmin([8,0.1,-9,pi])
3

julia> indmin([7,1,1,6])
2

julia> indmin([7,1,1,NaN])
2
```
"""
indmin(a) = findmin(a)[2]

# similar to Matlab's ismember
"""
    indexin(a, b)

Returns a vector containing the highest index in `b` for
each value in `a` that is a member of `b` . The output
vector contains 0 wherever `a` is not a member of `b`.

# Examples
```jldoctest
julia> a = ['a', 'b', 'c', 'b', 'd', 'a'];

julia> b = ['a','b','c'];

julia> indexin(a,b)
6-element Array{Int64,1}:
 1
 2
 3
 2
 0
 1

julia> indexin(b,a)
3-element Array{Int64,1}:
 6
 4
 3
```
"""
function indexin(a::AbstractArray, b::AbstractArray)
    bdict = Dict(zip(b, 1:length(b)))
    [get(bdict, i, 0) for i in a]
end

function _findin(a, b)
    ind  = Int[]
    bset = Set(b)
    @inbounds for (i,ai) in enumerate(a)
        ai in bset && push!(ind, i)
    end
    ind
end

# If two collections are already sorted, findin can be computed with
# a single traversal of the two collections. This is much faster than
# using a hash table (although it has the same complexity).
function _sortedfindin(v, w)
    viter, witer = eachindex(v), eachindex(w)
    out  = eltype(viter)[]
    i, j = start(viter), start(witer)
    if done(viter, i) || done(witer, j)
        return out
    end
    viteri, i = next(viter, i)
    witerj, j = next(witer, j)
    @inbounds begin
        vi, wj = v[viteri], w[witerj]
        while true
            if isless(vi, wj)
                if done(viter, i)
                    break
                end
                viteri, i = next(viter, i)
                vi        = v[viteri]
            elseif isless(wj, vi)
                if done(witer, j)
                    break
                end
                witerj, j = next(witer, j)
                wj        = w[witerj]
            else
                push!(out, viteri)
                if done(viter, i)
                    break
                end
                # We only increment the v iterator because v can have
                # repeated matches to a single value in w
                viteri, i = next(viter, i)
                vi        = v[viteri]
            end
        end
    end
    return out
end

"""
    findin(a, b)

Returns the indices of elements in collection `a` that appear in collection `b`.

# Examples
```jldoctest
julia> a = collect(1:3:15)
5-element Array{Int64,1}:
  1
  4
  7
 10
 13

julia> b = collect(2:4:10)
3-element Array{Int64,1}:
  2
  6
 10

julia> findin(a,b) # 10 is the only common element
1-element Array{Int64,1}:
 4
```
"""
function findin(a::Array{<:Real}, b::Union{Array{<:Real},Real})
    if issorted(a, Sort.Forward) && issorted(b, Sort.Forward)
        return _sortedfindin(a, b)
    else
        return _findin(a, b)
    end
end
# issorted fails for some element types so the method above has to be restricted
# to element with isless/< defined.
findin(a, b) = _findin(a, b)

# Copying subregions
# TODO: DEPRECATE FOR #14770
function indcopy(sz::Dims, I::Vector)
    n = length(I)
    s = sz[n]
    for i = n+1:length(sz)
        s *= sz[i]
    end
    dst = eltype(I)[findin(I[i], i < n ? (1:sz[i]) : (1:s)) for i = 1:n]
    src = eltype(I)[I[i][findin(I[i], i < n ? (1:sz[i]) : (1:s))] for i = 1:n]
    dst, src
end

function indcopy(sz::Dims, I::Tuple{Vararg{RangeIndex}})
    n = length(I)
    s = sz[n]
    for i = n+1:length(sz)
        s *= sz[i]
    end
    dst::typeof(I) = ntuple(i-> findin(I[i], i < n ? (1:sz[i]) : (1:s)), n)::typeof(I)
    src::typeof(I) = ntuple(i-> I[i][findin(I[i], i < n ? (1:sz[i]) : (1:s))], n)::typeof(I)
    dst, src
end

## Filter ##

"""
    filter(f, a::AbstractArray)

Return a copy of `a`, removing elements for which `f` is `false`.
The function `f` is passed one argument.

# Examples
```jldoctest
julia> a = 1:10
1:10

julia> filter(isodd, a)
5-element Array{Int64,1}:
 1
 3
 5
 7
 9
```
"""
filter(f, As::AbstractArray) = As[map(f, As)::AbstractArray{Bool}]

"""
    filter!(f, a::AbstractVector)

Update `a`, removing elements for which `f` is `false`.
The function `f` is passed one argument.

# Examples
```jldoctest
julia> filter!(isodd, collect(1:10))
5-element Array{Int64,1}:
 1
 3
 5
 7
 9
```
"""
function filter!(f, a::AbstractVector)
    isempty(a) && return a

    idx = eachindex(a)
    state = start(idx)
    i, state = next(idx, state)

    for acurr in a
        if f(acurr)
            a[i] = acurr
            i, state = next(idx, state)
        end
    end

    deleteat!(a, i:last(idx))

    return a
end

function filter(f, a::Vector)
    r = Vector{eltype(a)}(0)
    for ai in a
        if f(ai)
            push!(r, ai)
        end
    end
    return r
end

# set-like operators for vectors
# These are moderately efficient, preserve order, and remove dupes.

function intersect(v1, vs...)
    ret = Vector{promote_eltype(v1, vs...)}(0)
    for v_elem in v1
        inall = true
        for vsi in vs
            if !in(v_elem, vsi)
                inall=false; break
            end
        end
        if inall
            push!(ret, v_elem)
        end
    end
    ret
end

function union(vs...)
    ret = Vector{promote_eltype(vs...)}(0)
    seen = Set()
    for v in vs
        for v_elem in v
            if !in(v_elem, seen)
                push!(ret, v_elem)
                push!(seen, v_elem)
            end
        end
    end
    ret
end
# setdiff only accepts two args

"""
    setdiff(a, b)

Construct the set of elements in `a` but not `b`. Maintains order with arrays. Note that
both arguments must be collections, and both will be iterated over. In particular,
`setdiff(set,element)` where `element` is a potential member of `set`, will not work in
general.

# Examples
```jldoctest
julia> setdiff([1,2,3],[3,4,5])
2-element Array{Int64,1}:
 1
 2
```
"""
function setdiff(a, b)
    args_type = promote_type(eltype(a), eltype(b))
    bset = Set(b)
    ret = Array{args_type,1}(0)
    seen = Set{eltype(a)}()
    for a_elem in a
        if !in(a_elem, seen) && !in(a_elem, bset)
            push!(ret, a_elem)
            push!(seen, a_elem)
        end
    end
    ret
end
# symdiff is associative, so a relatively clean
# way to implement this is by using setdiff and union, and
# recursing. Has the advantage of keeping order, too, but
# not as fast as other methods that make a single pass and
# store counts with a Dict.
symdiff(a) = a
symdiff(a, b) = union(setdiff(a,b), setdiff(b,a))
"""
    symdiff(a, b, rest...)

Construct the symmetric difference of elements in the passed in sets or arrays.
Maintains order with arrays.

# Examples
```jldoctest
julia> symdiff([1,2,3],[3,4,5],[4,5,6])
3-element Array{Int64,1}:
 1
 2
 6
```
"""
symdiff(a, b, rest...) = symdiff(a, symdiff(b, rest...))

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
using BerghAlgorithms
using Polymake
using Test

 @testset "StackyFan Helper Functions" begin
    @test stackyWeights(makeStackyFan([1 0; 1 1; 1 2],[[0,1],[1,2]],[2,2,2]))==[ 2 ,  2 ,  2 ]

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 1 1; 1 2],INPUT_CONES=[[0,1],[1,2]])
    @test stackyWeights(addStackStructure(X,[2,2,2]))==[ 2 ,  2 ,  2 ]

    @test encode([1,0,2,5])=="1,0,2,5"

    F=makeStackyFan([1 0; 1 1; 1 2; 1 3],[[0,1],[1,2],[2,3]],[1,2,3,4]);
    @test stackyWeights(F)==[ 1 ,  2 ,  3 ,  4 ]

    F=makeStackyFan([1 0; 1 1; 1 2; 1 3],[[0,1],[1,2],[2,3]],[1,2,3,4]);
    @test getRayStack(F,[1,2])==3

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 1 0;1 0 1;1 1 1],INPUT_CONES=[[0,1,2],[1,2,3]]);
    SX = addStackStructure(X, [2,3,5,7]);
    @test stackyWeights(rootConstruction(SX, [1, 4, 2, 1]))==[ 2 ,  12 ,  10 ,  7 ]

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 1 0;1 0 1;1 1 1],INPUT_CONES=[[0,1,2],[1,2,3]]);
    SX = StackyFan(X, [2,3,5,7]);
    @test stackyWeights(rootConstructionDistinguishedIndices(SX, [0,1,1,0], [4, 2, 1, 3]))==[ 2 ,  6 ,  5 ,  7 ]

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 1 0;1 0 1;1 1 1],INPUT_CONES=[[0,1,2],[1,2,3]]);
    F = StackyFan(X, [2,3,5,7]);
    distinguished = X.RAYS[[2,3],:];
    @test stackyWeights(rootConstructionDistinguished(F, distinguished, [4, 2]))==[ 2 ,  12 ,  10 ,  7 ]

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0;1 1; 0 1],INPUT_CONES=[[0,1],[1,2]]);
    s=[1,2];
    @test convert(Array{Int64,1},vec(Array(findBarycenter(s,X))))==[2,1]

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 1 2],INPUT_CONES=[[0,1]]);
    F=addStackStructure(X,[2,3]);
    @test convert(Array{Int64,1},Array(findStackyBarycenter([1,2],F)))==[5,6]

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 2 3; 1 5],INPUT_CONES=[[0,1],[1,2]]);
    F=addStackStructure(X,[1,2,3]);
    @test findStackyRayMatrix(F)==[1 0;4 6;3 15]

    B=[true,true,false,true]
    @test convertBool(B)==[0, 1, 3]

    v=[1,2]
    M=[0 1; 1 1; 1 0]
    @test getConeRank(v,M)==2

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 2 0;0 0 1;0 1 0; 1 1 1],INPUT_CONES=[[0,1,2],[0,2,3,4]])
    @test getDimension(X)==3

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0; 1 1 0; 1 0 1; 1 1 1],INPUT_CONES=[[0,1,2,3]])
    @test getConeFaces(X,[1,2,3,4],Array(X.RAYS))==[[ 1 ,  2 ], [ 1 ,  3 ], [ 3 ,  4 ], [ 2 ,  4 ]]

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 1 0;1 0 1;1 1 1],INPUT_CONES=[[0,1,2,3]])
    @test X.SIMPLICIAL==false
    @test makeSimplicial(X).SIMPLICIAL==true
    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0 0;0 1 0 0;0 0 1 0;1 -1 1 0; 1 0 -2 0],INPUT_CONES=[[0,1,2,3],[0,4]])
    @test X.SIMPLICIAL==false
    @test makeSimplicial(X).SIMPLICIAL==true

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 1 2],INPUT_CONES=[[0,1]]);
    F=addStackStructure(X,[2,3]);
    @test stackyWeights(stackyBlowup(F,[0,1],[1,1]))==[2,1,3]

    F=makeStackyFan([1 0; 1 1; 1 2; 1 3],[[0,1],[1,2],[2,3]],[1,1,1,1])
    @test getConesPolymake(F)[1].RAYS==[1 0; 1 1]

    A=[1 2; 3 4];
    @test slicematrix(A)==[[1,2],[3,4]]

    A=[1 2 3;4 5 6; 7 8 9];
    S=Set([1,3]);
    @test rowMinors(A,S)==[1 2 3; 7 8 9]

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0;1 1; 0 1],INPUT_CONES=[[0,1],[1,2]]);
    M=X.MAXIMAL_CONES;
    @test convertIncidenceMatrix(M)==[[1,2],[2,3]]

    C=Polymake.polytope.Cone(INPUT_RAYS=[1 0; 1 2]);
    @test coneMultiplicity(C)==2

    @test getMultiplicity([1,2],[1 0; 1 2; 1 3])==2

    @test typeof(coneConvert([1, 2, 4],[1 0 0; 0 1 0; 0 0 1; 1 1 1]))==Polymake.BigObjectAllocated

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0; 1 1 0; 1 0 1; 1 1 1],INPUT_CONES=[[0,1,2,3]]);
    @test Set(getCones(X))==Set([[1,2,3,4],[1,2],[1,3],[3,4],[2,4],[1],[2],[3],[4]])

    C=Polymake.polytope.Cone(INPUT_RAYS=[1 0 0; 1 0 1; 1 1 0; 1 1 1]);
    @test convert(Array{Int64,2},Array(findFaceContainingRay(C,[1,1,1]).RAYS))==[1 0 1; 1 1 1]

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0 0; 1 0 1 0; 1 1 0 0; 1 1 1 0; 0 0 0 1],INPUT_CONES=[[0,1,2,3,4]]);
    @test findMinimalCone(X,[2,1,0,0])==[1,3]

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0 0; 1 0 1 0; 1 1 0 0; 1 1 1 0; 0 0 0 1],INPUT_CONES=[[0,1,2,3,4]]);
    @test convert(Array{Int64,2},Array(Polymake.common.primitive(starSubdivision(X,[2,1,0,0]).RAYS)))==[1 0 0 0;1 0 1 0;0 0 0 1;2 1 0 0; 1 1 0 0; 1 1 1 0]

    @test distinguishedAndIntPoint([1,2,4],[1 0 0; 1 2 0;2 1 3; 1 0 3],[1,0,0,0])==true

    @test convertToIncidence([2,3,5],6)==[0,1,1,0,1,0]

    @test compareCones([1,2],[2,3],[1 0 0; 0 1 0; 0 0 1],[1,1,0])==-1
    @test compareCones([1,2],[1,3],[1 0;1 2;1 -1],[1,1,1])==1

    @test extremalCones([[1,2],[2,3],[3,4]],[1 0;1 2; 1 5; 1 8],[0,1,1,0])==[[3,4]]

    C=Polymake.polytope.Cone(INPUT_RAYS=[1 2; 2 1]);
    @test interiorPoints(C)==[[1,1]]

    A=[[1,1,1],[2,1,3],[0,5,4]];
    @test minimalByLex(A)==[0,5,4]

    @test minimalByDist([[0,1,5,7],[3,3,2,2],[8,5,3,6],[2,1,1,10]],[0,1,1,0])==[2,1,1,10]

    @test coneRayDecomposition([1,2,3],[3 5 7; 8 16 9;2 1 3;1 1 1],[2,2,3],[1,1,1,1])==[6,5,52,0]

end

@testset "Visualization" begin
    plot3dSimpCone([[1,0,0],[0,1,0],[0,0,1]]);

    plot2dCone([[1,0],[0,1]]);

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 2 5],INPUT_CONES=[[0,1]]);
    showSimpFan(X);

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 2 5],INPUT_CONES=[[0,1]]);
    F=addStackStructure(X,[1,2]);
    showSimpStackyFan(F);

    C=Polymake.polytope.Cone(INPUT_RAYS=[1 0 0; 1 0 1; 1 1 0; 1 1 1]);
    @test coneVectorOrder(C)==[[1,1,0],[1,1,1],[1,0,1],[1,0,0]]

    plot3dCone([[1,1,0],[1,1,1],[1,0,1],[1,0,0]]);

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0; 1 0 1; 1 1 0; 1 1 1],INPUT_CONES=[[0,1,2,3]]);
    showFan(X);

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0; 1 0 1; 1 1 0; 1 1 1],INPUT_CONES=[[0,1,2,3]]);
    F=addStackStructure(X,[1,2,2,5]);
    showStackyFan(F);
end

 @testset "BerghA.jl" begin
     @testset "test stackyWeights" begin
         X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 2 5],INPUT_CONES=[[0,1]]);
         F=addStackStructure(X,[1,1]);
         @test BerghA(F,[1,0]).stacks["1,2"] == 5
         @test stackyWeights(BerghA(F,[1,1])) == [ 5 ,  2 ,  5 ,  10 ]
         @test stackyWeights(BerghA(F,[1,0])) == [ 5 ,  5 ,  1 ,  2 ,  5 ,  10 ]
         @test stackyWeights(BerghA(F,[0,1])) == [ 1 ,  5 ,  5 ,  2 ,  5 ,  1 ,  5 ,  10 ]

         X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 3; 4 5 6; 2 3 1],INPUT_CONES=[[0,1,2]])
         F=addStackStructure(X,[1,1,1])
         @test stackyWeights(BerghA(F,[1,1,1])) == [ 28 ,  21 ,  84 ,  28 ,  84 ,  84 ,  42 ,  84 ,  168 ]

         X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 2; 2 1 1; 5 3 9],INPUT_CONES=[[0,1,2]])
         F=addStackStructure(X,[1,1,1])
         @test stackyWeights(BerghA(F,[1,1,1])) == [ 4 ,  4 ,  8 ,  4 ,  4 ,  8 ]

         X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 4 3; 1 5],INPUT_CONES=[[0,1],[1,2]])
         F=addStackStructure(X,[1,1,1])
         @test stackyWeights(BerghA(F,[1,1,1])) == [ 4 ,  34 ,  6 ,  68 ,  34 ,  68 ,  34 ,  6 ,  12 ]

         X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[323 179; 44 135],INPUT_CONES=[[0,1]])
         F=addStackStructure(X,[1,1])
         @test stackyWeights(BerghA(F,[1,1])) == [ 491602456800 ,  49160245680 ,  294961474080 ,  468192816 ,  12173013216 ,  73038079296 ,  12173013216 ,  97384105728 ,  73038079296 ,  196640982720 ,  131093988480 ,  156064272 ,  936385632 ,  468192816 ,  196640982720 ,  49160245680 ,  9832049136 ,  292152317184 ,  5899229481600 ,  3932819654400 ,  589922948160 ,  393281965440 ,  12173013216 ,  8115342144 ,  12173013216 ,  5899229481600 ,  589922948160 ,  737403685200 ,  51126655507200 ,  292152317184 ,  51126655507200 ,  5899229481600 ,  5899229481600 ,  589922948160 ,  589922948160 ,  196640982720 ,  2949614740800 ,  5899229481600 ,  294961474080 ,  589922948160 ,  196640982720 ,  589922948160 ,  393281965440 ,  936385632 ,  24346026432 ,  24346026432 ,  11798458963200 ,  1179845896320 ,  737403685200 ,  1474807370400 ,  2949614740800 ,  5899229481600 ,  589922948160 ,  11798458963200 ,  1179845896320 ,  2949614740800 ,  1474807370400 ,  5899229481600 ,  102253311014400 ]

    end
 end

# # @testset "slow test" begin
# #          X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[4 1; 7 9],INPUT_CONES=[[0,1]])
# #         F=addStackStructure(X,[1,1])
# #         @test stackyWeights(BerghA(F,[1,0])) == [ 609 ,  29 ,  1 ,  174 ,  29 ,  1740 ,  1218 ,  58 ,  145 ,  1044 ,  1044 ,  290 ,  290 ,  145 ,  406 ,  348 ,  261 ,  14616 ,  14616 ,  609 ,  3480 ,  870 ,  609 ,  174 ,  609 ,  9744 ,  14616 ,  1218 ,  58 ,  145 ,  435 ,  725 ,  1305 ,  1740 ,  6960 ,  3480 ,  870 ,  3480 ,  58464 ,  6090 ,  3480 ,  1392 ,  696 ,  1044 ,  2088 ,  261 ,  174 ,  261 ,  609 ,  406 ,  609 ,  609 ,  406 ,  609 ,  1218 ,  812 ,  1218 ,  2088 ,  261 ,  1044 ,  1160 ,  1740 ,  1740 ,  870 ,  1305 ,  1305 ,  14616 ,  10440 ,  145 ,  435 ,  609 ,  116 ,  580 ,  290 ,  580 ,  1740 ,  3480 ,  3480 ,  261 ,  522 ,  261 ,  522 ,  522 ,  1218 ,  1218 ,  1218 ,  1218 ,  2436 ,  2436 ,  4176 ,  2088 ,  3480 ,  2610 ,  29232 ,  6090 ,  1218 ,  290 ,  145 ,  290 ,  12180 ,  261 ,  522 ]
# # end

@testset "makeSmooth tests" begin
    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[4 -1; 0 1],INPUT_CONES=[[0, 1]])
    @test X.SMOOTH_FAN == false
    @test makeSmooth(X).SMOOTH_FAN == true
    @test makeSmooth(X).INPUT_RAYS == [1 -1/4; 0 1; 1 0]

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 1 0;1 0 1;1 1 1],INPUT_CONES=[[0,1,2,3]])
    @test X.SIMPLICIAL == false
    @test X.SMOOTH_FAN == false
    @test makeSmooth(X).SMOOTH_FAN == true

    #This test fails, the last ray is given as 4 2 2
    @test makeSmooth(X).INPUT_RAYS == [1 0 0;1 1 0;1 0 1;1 1 1;2 1 1]
    @test makeSmooth(X).SIMPLICIAL == true

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0 0;0 1 0 0;0 0 1 0;1 -1 1 0;1 0 -2 0],INPUT_CONES=[[0,1,2,3],[0,4]])
    @test X.SIMPLICIAL == false
    @test X.SMOOTH_FAN == false
    @test makeSmooth(X).SMOOTH_FAN == true
    @test makeSmooth(X).SIMPLICIAL == true

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 2 3;-1 1 1;1 -1 1;-1 -1 1;1 1 -1;-1 1 -1;1 -1 -1;-1 -1 -1],INPUT_CONES=[[0,1,2,3],[0,1,4,5],[0,2,4,6],[1,3,5,7],[2,3,6,7],[4,5,6,7]])
    @test X.SIMPLICIAL == false
    @test X.SMOOTH_FAN == false
    @test makeSmooth(X).SMOOTH_FAN == true
    @test makeSmooth(X).SIMPLICIAL == true

    X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;0 1 0;0 0 1;0 -1 -1; -1 0 -1; -2 -1 0],INPUT_CONES=[[0,1,2],[0,1,3],[1,3,4],[1,2,4],[2,4,5],[0,2,5],[0,3,5],[3,4,5]])
    @test X.SIMPLICIAL == true
    @test X.SMOOTH_FAN == false
    @test makeSmooth(X).SMOOTH_FAN == true
    @test makeSmooth(X).SIMPLICIAL == true

end

@testset "toric_blowup tests" begin
       X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 1 0;1 0 1;1 1 1],INPUT_CONES=[[0,1,2,3]])
       @test X.SIMPLICIAL == false
       A = toric_blowup([0,1,2,3], X, nothing)
       @test A.INPUT_RAYS == [1 0 0; 1 1 0; 1 0 1; 1 1 1;2 1 1]
       @test Set(convertIncidenceMatrix(A.INPUT_CONES)) == Set([[1, 2, 5], [1, 3, 5], [3, 4, 5], [2, 4, 5]])
       @test A.SIMPLICIAL == true
       @test A.SMOOTH_FAN == true

       A = toric_blowup([0,1], X, nothing)
       @test A.INPUT_RAYS == [1 0 0; 1 1 0; 1 0 1;1 1 1;2 1 0]
       @test Set(convertIncidenceMatrix(A.INPUT_CONES)) == Set([[1, 3, 5], [2, 4, 5], [3, 4, 5]])
       @test A.SIMPLICIAL == true

       A = toric_blowup([0,1,2,3], X, [5,3,4])
       @test A.INPUT_RAYS == [1 0 0;1 1 0;1 0 1;1 1 1;5 3 4]
       @test Set(convertIncidenceMatrix(A.INPUT_CONES)) == Set([[1, 2, 5], [1, 3, 5], [2, 4, 5], [3, 4, 5]])
       @test A.SIMPLICIAL == true

        A = toric_blowup([0], X, nothing)
        @test A.INPUT_RAYS == [1 0 0;1 1 0;1 0 1;1 1 1]
        @test Set(convertIncidenceMatrix(A.INPUT_CONES)) == Set([[1, 2, 4], [1, 3, 4]])
        @test A.SIMPLICIAL == true

       X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 1 2],INPUT_CONES=[[0,1]])
       B=toric_blowup([0,1],X,[1,1]);
       @test B.INPUT_RAYS == [1 0;1 2; 1 1]
       @test Set(convertIncidenceMatrix(B.INPUT_CONES)) == Set([[1, 3], [2, 3]])
       @test B.SIMPLICIAL == true
       @test B.SMOOTH_FAN == true

end

@testset "AlgC.jl" begin
    A=[1,2,3,4];
    remove!(A,1);
    @test A==[2,3,4]

    @test getIndex([0,1,0],[1 0 0; 0 1 0; 0 0 1])==2

    @test isIndependent(3,[1,2,3],[1 0 0; 0 1 0; 1 2 3])==false
    @test isIndependent(3,[1,2,3],[1 0 0; 0 1 0; 1 1 1])==true

    @test independencyIndex([1,2,3],[1 0 0 ; 1 2 0; 2 0 3; 0 0 5])==3

    F=makeStackyFan([1 0 0; 1 2 0; 0 0 1],[[0,1,2]],[1,1,2]);
    @test isRelevant([1,2,0],[1,2,3],F)==true
    F=makeStackyFan([1 0 0; 0 1 0; 0 0 1],[[0,1,2]],[1,1,2]);
    @test isRelevant([0,1,0],[1,2,3],F)==false

    F=makeStackyFan([1 0 0; 0 1 0; 1 0 2],[[0,1,2]],[1,2,1]);
    div=Dict([1,0,0]=>0,[0,1,0]=>0,[1,0,2]=>0);
    @test toroidalIndex([1,2,3],F,div)==3
    div=Dict([1,0,0]=>0,[0,1,0]=>0,[1,0,2]=>1);
    @test toroidalIndex([1,2,3],F,div)==2

    F=makeStackyFan([1 0 0; 0 1 0; 1 0 2],[[0,1,2]],[1,2,1]);
    div=Dict([1,0,0]=>0,[0,1,0]=>0,[1,0,2]=>0);
    @test divisorialIndex([1,2,3],F,div)==3
    div=Dict([1,0,0]=>0,[0,1,0]=>0,[1,0,2]=>1);
    @test divisorialIndex([1,2,3],F,div)==1
    H, div = BerghC(F,[0,0,0]);
    @test convert(Array{Int64,2},Polymake.common.primitive(H.fan.RAYS))== [1 0 0; 0 1 0; 1 1 1; 1 0 2; 1 0 1]
    @test div == Dict([1, 1, 1] => 1, [1, 0, 0] => 0, [1, 0, 1] => 1, [0, 1, 0] => 1, [1, 0, 2] => 0)

    @test coneContains([1,2,3],[1,2,3,4])==true
    @test coneContains([1,2,5],[1,2,3,4])==false

    F=makeStackyFan([1 2 0;1 3 0; 3 0 1],[[0,1,2]],[1,1,5]);
    div=Dict([1,2,0]=>0,[1,3,0]=>0,[3,0,1]=>0);
    @test minMaxDivisorial(F,div)==[[3]]
    F=makeStackyFan([1 1 0;1 3 0; 3 0 1],[[0,1,2]],[1,1,5]);
    div=Dict([1,1,0]=>0,[1,3,0]=>0,[3,0,1]=>0);
    @test minMaxDivisorial(F,div)==[[1,2,3]]

    F=makeStackyFan([1 1 0;1 3 0; 0 0 1],[[0,1,2]],[1,1,5]);
    H, div = BerghC(F,[0,0,0]);
    @test convert(Array{Int64,2},Polymake.common.primitive(H.fan.RAYS))==[1 1 0; 0 0 1; 2 4 5; 1 3 0; 1 2 0]
    @test div==Dict([1,1,0]=>0,[0,0,1]=>1,[2,4,5]=>1,[1,2,0]=>1,[1,3,0]=>0)

    F=makeStackyFan([1 0;1 3; 5 17],[[0,1],[1,2]],[1,1,5]);
    H, div = BerghC(F,[0,0,0]);
    @test convert(Array{Int64,2},Polymake.common.primitive(H.fan.RAYS))==[1 0; 2 3; 1 3; 13 44; 5 17]
    @test div==Dict([1,0]=>0,[13,44]=>1,[1,3]=>0,[2,3]=>1,[5,17]=>1)

    F=makeStackyFan([1 0;1 3],[[0,1]],[1,1]);
    H, div = BerghC(F,[0,0,0]);
    @test convert(Array{Int64,2},Polymake.common.primitive(H.fan.RAYS))==[1 0; 2 3; 1 3]
    @test div==Dict([1,0]=>0,[1,3]=>0,[2,3]=>1)

    H, div = BerghC(F, [1,1,1]);
    @test convert(Array{Int64,2},Polymake.common.primitive(H.fan.RAYS))==[1 0; 1 3]
    @test div==Dict([1,0]=>1,[1,3]=>1)

    F=makeStackyFan([1 0;1 3],[[0,1]],[3,3]);
    H, div = BerghC(F,[0,0,0]);
    @test convert(Array{Int64,2},Polymake.common.primitive(H.fan.RAYS))==[1 0; 2 3; 1 3]
    @test div==Dict([1, 0] => 1, [1, 3] => 1, [2, 3] => 1)

    F=makeStackyFan([0 0 1;1 0 0; 1 3 0],[[0,1,2]],[1,1,1]);
    H, div = BerghC(F,[0,0,0]);
    @test convert(Array{Int64,2},Polymake.common.primitive(H.fan.RAYS))==[0 0 1; 1 0 0; 2 3 0; 1 3 0]
    @test div==Dict([0, 0, 1] => 0, [2, 3, 0] => 1, [1, 0, 0] => 0, [1, 3, 0] => 0)

    F=makeStackyFan([0 0 1;1 0 0; 1 3 0],[[0,1,2]],[1,2,5]);
    H, div = BerghC(F,[0,0,0]);
    @test convert(Array{Int64,2},Polymake.common.primitive(H.fan.RAYS))==[0 0 1; 1 0 0; 7 15 0; 1 3 0]
    @test div==Dict([7, 15, 0] => 1, [0, 0, 1] => 0, [1, 0, 0] => 1, [1, 3, 0] => 1)

    F = makeStackyFan([1 0; 1 1; 1 2],[[0,1],[1,2]],[2,2,2]);
    H, div = BerghC(F,[0,0,0]);
    @test convert(Array{Int64,2},Polymake.common.primitive(H.fan.RAYS)) ==[1 0; 2 1; 1 1; 2 3; 1 2]
    @test div==Dict([1, 0] => 1, [1, 2] => 1, [2, 1] => 1, [1, 1] => 1, [2, 3] => 1)

end

max_stars_repo_name: nbenabla/destackification

max_issues_repo_name: nbenabla/destackification

max_forks_repo_name: nbenabla/destackification


text:
# Extend the graph utility functions to support Javis `JGraph`
for extend in [
    :is_directed,
    :edgetype,
    :ne,
    :nv,
    :outneighbors,
    :inneighbors,
    :has_vertex,
    :has_edge,
]
    eval(
        quote
            local func = $extend
            Graphs.$extend(jg::JGraph, args...) = $extend(jg.graph.adjacency_graph, args...)
        end
    )
end

"""
    add_vertex!()

Add a vertex to the canvas. The layout is regenerated depending upon the `mode` of the graph selected.

The function syntax chosen is similar to `Graphs.add_vertex!`.
"""
function Graphs.add_vertex!(jg::JGraph)
end

"""
    rem_vertex!()

Remove a vertex from the canvas. The layout is regenerated depending upon the `mode` of the graph selected.

The function syntax chosen is similar to `Graphs.rem_vertex!`.
"""
function Graphs.rem_vertex!(jg::JGraph, vertex::T) where {T}
end

"""
    add_edge!()

Add an edge to the canvas.

The function syntax chosen is similar to `Graphs.add_edge!`.
"""
function Graphs.add_edge!(jg::JGraph)
end

"""
    rem_edge!()

Remove an edge from the canvas.

The function syntax chosen is similar to `Graphs.rem_edge!`.
"""
function Graphs.rem_edge!(jg::JGraph, from_vertex::T, to_vertex::T) where {T}
end