text:
# broadcasting operations

export AbstractGPUArrayStyle

using Base.Broadcast

import Base.Broadcast: BroadcastStyle, Broadcasted, AbstractArrayStyle

const BroadcastGPUArray{T} = Union{AnyGPUArray{T},
                                   Base.RefValue{<:AbstractGPUArray{T}}}

"""
Abstract supertype for GPU array styles. The `N` parameter is the dimensionality.

Downstream implementations should provide a concrete array style type that inherits from
this supertype.
"""
abstract type AbstractGPUArrayStyle{N} <: AbstractArrayStyle{N} end

# Wrapper types otherwise forget that they are GPU compatible
# NOTE: don't directly use GPUArrayStyle here not to lose downstream customizations.
BroadcastStyle(W::Type{<:WrappedGPUArray})= BroadcastStyle(Adapt.parent(W){Adapt.eltype(W), Adapt.ndims(W)})
backend(W::Type{<:WrappedGPUArray}) = backend(Adapt.parent(W){Adapt.eltype(W), Adapt.ndims(W)})

# Ref is special: it's not a real wrapper, so not part of Adapt,
# but it is commonly used to bypass broadcasting of an argument
# so we need to preserve its dimensionless properties.
BroadcastStyle(::Type{Base.RefValue{AT}}) where {AT<:AbstractGPUArray} =
    typeof(BroadcastStyle(AT))(Val(0))
backend(::Type{Base.RefValue{AT}}) where {AT<:AbstractGPUArray} = backend(AT)
# but make sure we don't dispatch to the optimized copy method that directly indexes
function Broadcast.copy(bc::Broadcasted{<:AbstractGPUArrayStyle{0}})
    ElType = Broadcast.combine_eltypes(bc.f, bc.args)
    isbitstype(ElType) || error("Cannot broadcast function returning non-isbits $ElType.")
    dest = copyto!(similar(bc, ElType), bc)
    return @allowscalar dest[CartesianIndex()]  # 0D broadcast needs to unwrap results
end

# we need to override the outer copy method to make sure we never fall back to scalar
# iteration (see, e.g., CUDA.jl#145)
@inline function Broadcast.copy(bc::Broadcasted{<:AbstractGPUArrayStyle})
    ElType = Broadcast.combine_eltypes(bc.f, bc.args)
    if !Base.isconcretetype(ElType)
        error("""GPU broadcast resulted in non-concrete element type $ElType.
                 This probably means that the function you are broadcasting contains an error or type instability.""")
    end
    copyto!(similar(bc, ElType), bc)
end

@inline function Base.copyto!(dest::BroadcastGPUArray, bc::Broadcasted{Nothing})
    axes(dest) == axes(bc) || Broadcast.throwdm(axes(dest), axes(bc))
    isempty(dest) && return dest
    bc′ = Broadcast.preprocess(dest, bc)

    # grid-stride kernel
    function broadcast_kernel(ctx, dest, bc′, nelem)
        for i in 1:nelem
            I = @cartesianidx(dest, i)
            @inbounds dest[I] = bc′[I]
        end
        return
    end
    elements = length(dest)
    elements_per_thread = typemax(Int)
    heuristic = launch_heuristic(backend(dest), broadcast_kernel, dest, bc′, 1;
                                 elements, elements_per_thread)
    config = launch_configuration(backend(dest), heuristic;
                                  elements, elements_per_thread)
    gpu_call(broadcast_kernel, dest, bc′, config.elements_per_thread;
             threads=config.threads, blocks=config.blocks)

    return dest
end

# Base defines this method as a performance optimization, but we don't know how to do
# `fill!` in general for all `BroadcastGPUArray` so we just go straight to the fallback
@inline Base.copyto!(dest::BroadcastGPUArray, bc::Broadcasted{<:Broadcast.AbstractArrayStyle{0}}) =
    copyto!(dest, convert(Broadcasted{Nothing}, bc))

## map

allequal(x) = true
allequal(x, y, z...) = x == y && allequal(y, z...)

function Base.map(f, x::BroadcastGPUArray, xs::AbstractArray...)
    # if argument sizes match, their shape needs to be preserved
    xs = (x, xs...)
    if allequal(size.(xs)...)
         return f.(xs...)
    end

    # if not, treat them as iterators
    indices = LinearIndices.(xs)
    common_length = minimum(length.(indices))

    # construct a broadcast to figure out the destination container
    ElType = Broadcast.combine_eltypes(f, xs)
    isbitstype(ElType) || error("Cannot map function returning non-isbits $ElType.")
    dest = similar(x, ElType, common_length)

    return map!(f, dest, xs...)
end

function Base.map!(f, dest::BroadcastGPUArray, xs::AbstractArray...)
    # custom broadcast, ignoring the container size mismatches
    # (avoids the reshape + view that our mapreduce impl has to do)
    indices = LinearIndices.((dest, xs...))
    common_length = minimum(length.(indices))
    common_length==0 && return

    bc = Broadcast.instantiate(Broadcast.broadcasted(f, xs...))
    if bc isa Broadcast.Broadcasted
        bc = Broadcast.preprocess(dest, bc)
    end

    # grid-stride kernel
    function map_kernel(ctx, dest, bc, nelem)
        for i in 1:nelem
            j = linear_index(ctx, i)
            j > common_length && return

            J = CartesianIndices(axes(bc))[j]
            @inbounds dest[j] = bc[J]
        end
        return
    end
    elements = common_length
    elements_per_thread = typemax(Int)
    heuristic = launch_heuristic(backend(dest), map_kernel, dest, bc, 1;
                                 elements, elements_per_thread)
    config = launch_configuration(backend(dest), heuristic;
                                  elements, elements_per_thread)
    gpu_call(map_kernel, dest, bc, config.elements_per_thread;
             threads=config.threads, blocks=config.blocks)

    return dest
end


text:
using PtFEM

data = Dict(
  # Frame(nels, nn, ndim, np_types, nip, finite_element(nod, nodof))
  :struc_el => Frame(10, 6, 2, 1, 1, Line(2, 2)),
  :properties => [2.0e5;],
  :x_coords => [0.0, 4.0, 4.0, 8.0, 8.0, 12.0],
  :y_coords => [3.0, 0.0, 3.0, 3.0, 0.0, 0.0],
  :g_num => [
    1 1 3 3 3 2 2 5 4 5;
    2 3 4 5 2 4 5 4 6 6
  ],
  :support => [
    (1, [0 0]),
    (2, [1 0])
    ],
  :loaded_nodes => [
    (6, [0.0 -10.0])],
  :penalty => 1e19
)

data |> display
println()

@time fem, dis_df, fm_df = p42(data)
println()

display(dis_df)
println()
display(fm_df)

ext: jl


text:
"Exact solution of an ODE."
struct IntegratorExactODE{DT, TT, D, QT <: Function} <: ODEIntegrator{DT,TT}
    q::QT
    Δt::TT

    function IntegratorExactODE{DT,D}(solution::solType, Δt::TT) where {DT, TT, D, solType <: Function}
        new{DT,TT,D,solType}(solution, Δt)
    end
end

@inline Base.ndims(::IntegratorExactODE{DT,TT,D}) where {DT,TT,D} = D
timestep(int::IntegratorExactODE) = int.Δt

function integrate_step!(int::IntegratorExactODE{DT,TT}, sol::AtomicSolutionODE{DT,TT}) where {DT,TT}
    # reset atomic solution
    reset!(sol, timestep(int))

    # compute new solution
    int.q(sol.t̄ + timestep(int), sol.q̄, sol.q, timestep(int))
end


text:
# typealias Basis{T<:ScalarFunc3d{Float64},N} Array{T,N}

typealias Basis{T<:Func,N} AbstractArray{T,N}

_nindices(sym, n) = [symbol(sym, :_, i) for i = 1:n]  # e.g. i_1,...,i_N
_nref(sym, ind...) = Expr(:ref, sym, ind...)  # e.g. A[i_1,...,i_N]

function Base.super(T,S)
    @assert issubtype(T,S)
    while issubtype(super(T),S) && !is(T,Any)
        T = super(T)
    end
    T
end

for f in (:(Base.call), :laplacian)

    eval(quote

        @generated function $f{Ti,To,T<:Func,N,M}(b::Basis{T,N}, dst::AbstractArray{To}, c::AbstractArray{Ti,M})
            Ti_,To_ = super(T,Func).parameters
            @assert Ti == Ti_ && To == To_ "invalid types for $(super(T,Func))"
            i      = _nindices(:i,M)
            j      = _nindices(:j,N)
            c_i    = _nref(:c,i...)
            b_j    = _nref(:b,j...)
            dst_ij = _nref(:dst,i...,j...)
            quote
                Base.Cartesian.@nloops $N j b begin
                    Base.Cartesian.@nloops $M i c begin
                        $dst_ij = $($f)($b_j,$c_i)
                    end
                end
            end
        end

        @generated function $f{Ti,T<:Func,N,M}(b::Basis{T,N}, c::AbstractArray{Ti,M})
            Ti_,To = super(T,Func).parameters
            @assert Ti == Ti_ "invalid input type for $(super(T,Func))"
            quote
                dst = Array{$To}((size(c)...,size(b)...))
                $($f)(b,dst,c)
                dst
            end
        end

        @generated function $f{Ti,T<:Func,N}(b::Basis{T,N}, c::Ti)
            Ti_,To = super(T,Func).parameters
            @assert Ti == Ti_ "invalid input type for $(super(T,Func))"
            quote
                dst = Array{$To}(size(b))
                for i in eachindex(b)
                    dst[i] = $($f)(b[i],c)
                end
                dst
            end
        end

        @generated function $f{T<:Func,N,X}(b::Basis{T,N}, c::X...)
            Ti,To = super(T,Func).parameters
            @assert issubtype(Ti,Point) && X == first(Ti.parameters) "invalid input type for $(super(T,Func))"
            quote
                dst = Array{$To}(size(b))
                for i in eachindex(b)
                    dst[i] = $($f)(b[i],c...)
                end
                dst
            end
        end

    end)

end

# @generated function Base.call{T,N,M}(b::Basis{T,M}, dst::AbstractArray{Float64}, c::AbstractArray{Point3{Float64},N})
#     i      = _nindices(:i,N)
#     j      = _nindices(:j,M)
#     c_i    = _nref(:c,i...)
#     b_j    = _nref(:b,j...)
#     dst_ij = _nref(:dst,i...,j...)
#     quote
#         Base.Cartesian.@nloops $M j b begin
#             Base.Cartesian.@nloops $N i c begin
#                 $dst_ij = $b_j($c_i)
#             end
#         end
#     end
# end

# function Base.call(b::Basis, c::AbstractArray{Point3{Float64}})
#     dst = Array{Float64}((size(c)...,size(b)...))
#     b(dst,c)
#     dst
# end

# Base.call(b::Basis, x::Float64, y::Float64, z::Float64) = reshape(Float64[b[i](x,y,z) for i in eachindex(b)], size(b))
# Base.call(b::Basis, p::Point3{Float64}) = reshape(Float64[b[i](p) for i in eachindex(b)], size(b))

# @generated function laplacian{T,N,M}(b::Basis{T,M}, dst::AbstractArray{Float64}, c::AbstractArray{Point3{Float64},N})
#     i      = _nindices(:i,N)
#     j      = _nindices(:j,M)
#     c_i    = _nref(:c,i...)
#     b_j    = _nref(:b,j...)
#     dst_ij = _nref(:dst,i...,j...)
#     quote
#         Base.Cartesian.@nloops $M j b begin
#             Base.Cartesian.@nloops $N i c begin
#                 $dst_ij = laplacian($b_j,$c_i)
#             end
#         end
#     end
# end

# function laplacian(b::Basis, c::AbstractArray{Point3{Float64}})
#     dst = Array{Float64}((size(c)...,size(b)...))
#     laplacian(b,dst,c)
#     dst
# end

# @generated function add!{T,N,M}(b::Basis{T,M}, dst::AbstractArray{Float64}, c::AbstractArray{Point3{Float64},N})
#     i     = _nindices(:i,N)
#     j     = _nindices(:j,M)
#     c_i   = _nref(:c,i...)
#     b_j   = _nref(:b,j...)
#     dst_i = _nref(:dst,i...)
#     quote
#         @assert size(dst) == size(c)
#         Base.Cartesian.@nloops $N i c begin
#             Base.Cartesian.@nloops $M j b begin
#                 $dst_i += $b_j($c_i)
#             end
#         end
#     end
# end

# @generated function sum!{N,M}(b::Basis{M}, dst::AbstractArray{Float64}, c::AbstractArray{Point3{Float64},N})
#     i = _nindices(:i,N)
#     j = _nindices(:j,M)
#     c_i = _nref(:c,i...)
#     b_j = _nref(:b,j...)
#     dst_i = _nref(:dst,i...)
#     quote
#         @assert size(dst) == size(c)
#         Base.Cartesian.@nloops $N i c begin
#             $dst_i = 0.0
#             Base.Cartesian.@nloops $M j b begin
#                 $dst_i += laplacian($b_j,$c_i)
#             end
#         end
#     end
# end

# function sum(b::Basis, c::AbstractArray{Point3{Float64}})
#     dst = zeros(c)
#     add!(b,dst,c)
#     dst
# end

max_stars_repo_name: johncsnyder/Element


text:
```Implement the following sorting algorithms: Selection sort, Insertion sort,
 Merge sort, Quick sort, Stooge Sort. Check Wikipedia for descriptions.```

v = rand(-10:10, 1000)

 function selectionsort!(arr::Vector{<:Real})
    len = length(arr)
    if len < 2 return arr end
    for i in 1:len-1
        lmin, j = findmin(arr[i+1:end])
        if lmin < arr[i]
            arr[i+j] = arr[i]
            arr[i] = lmin
        end
    end
    return arr
end

@time  selectionsort!(v)

function insertionsort!(A::Array{T}) where T <: Number
    for i in 1:length(A)-1
        value = A[i+1]
        j = i
        while j > 0 && A[j] > value
            A[j+1] = A[j]
            j -= 1
        end
        A[j+1] = value
    end
    return A
end

x = rand(10)
@show x insertionsort!(x)

function mergesort(arr::Vector)
    if length(arr) ≤ 1 return arr end
    mid = length(arr) ÷ 2
    lpart = mergesort(arr[1:mid])
    rpart = mergesort(arr[mid+1:end])
    rst = similar(arr)
    i = ri = li = 1
    @inbounds while li ≤ length(lpart) && ri ≤ length(rpart)
        if lpart[li] ≤ rpart[ri]
            rst[i] = lpart[li]
            li += 1
        else
            rst[i] = rpart[ri]
            ri += 1
        end
        i += 1
    end
    if li ≤ length(lpart)
        copy!(rst, i, lpart, li)
    else
        copy!(rst, i, rpart, ri)
    end
    return rst
end

mergesort(v)

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
using StateSpaceReconstruction
using PerronFrobenius
using StaticArrays
using DelayEmbeddings

Vpts = [rand(4) for i = 1:50]
Dpts = Dataset([rand(4) for i = 1:50])
Spts = Dataset([rand(4) for i = 1:50])
Mpts = Dataset([rand(4) for i = 1:50])
pts = rand(50, 4)
Tpts = rand(4, 50)

@test transferoperator(pts, RectangularBinning(3)) isa TransferOperatorRectangularBinning
@test transferoperator(Dpts, RectangularBinning(3)) isa TransferOperatorRectangularBinning
@test transferoperator(Spts, RectangularBinning(3)) isa TransferOperatorRectangularBinning
@test transferoperator(Mpts, RectangularBinning(3)) isa TransferOperatorRectangularBinning
@test transferoperator(pts, RectangularBinning(3)) isa TransferOperatorRectangularBinning
@test transferoperator(Tpts, RectangularBinning(3)) isa TransferOperatorRectangularBinning

@testset "Transfer operator from rectangular binnings" begin
    @testset "Grid approach" begin
        @testset "Independent realization $i" for i=1:3
            points_2D = rand(2, 200)
            points_3D = rand(3, 400)
            points_4D = rand(4, 500)

            E_2D = invariantize(cembed(points_2D))
            E_3D = invariantize(cembed(points_3D))
            E_4D = invariantize(cembed(points_4D))

            points_2D = E_2D.points
            points_3D = E_3D.points
            points_4D = E_4D.points

            # The binning scheme
            ϵ = 3

            @testset "From precomputed bin visits" begin
                @testset "2D" begin
                    bins_visited_by_orbit = assign_bin_labels(E_2D, ϵ)
                    bininfo = get_binvisits(bins_visited_by_orbit)
                    TO = estimate_transferoperator_from_binvisits(bininfo)

                    @test typeof(TO) <: TransferOperatorRectangularBinning

                    # Last row might sum to zero, because the last point does not need to
                    # be contained in the last bin. However, the remaining row sums must
                    # be one.
                    if !is_markov(TO)
                        @warn "There were all-zero columns in the transfer matrix"
                        @warn "Removing first column and last row"
                        if is_markov(TO.transfermatrix[1:(end-1), 2:end])
                            info("That made the transfer matrix Markov")
                        else
                            @warn "That did NOT make the transfer matrix Markov"
                            #@save "from_binvisits_2D.jld2" TO, E_2D, bininfo
                        end
                        @test is_markov(TO.transfermatrix[1:(end-1), 2:end])
                    end
                end

                @testset "3D" begin
                    bins_visited_by_orbit = assign_bin_labels(E_3D, ϵ)
                    bininfo = get_binvisits(bins_visited_by_orbit)
                    TO = estimate_transferoperator_from_binvisits(bininfo)
                    @test typeof(TO) <: TransferOperatorRectangularBinning

                    # Last row might sum to zero, because the last point does not need to
            		# be contained in the last bin. However, the remaining row sums must
            		# be one.
                    if !is_markov(TO)
                        @warn "There were all-zero columns in the transfer matrix"
                        @warn "Removing first column and last row"
                        if is_markov(TO.transfermatrix[1:(end-1), 2:end])
                            info("That made the transfer matrix Markov")
                        else
                            @warn "That did NOT make the transfer matrix Markov"
                            #@save "from_binvisits_3D.jld2" TO, E_3D, bininfo
                        end
                        @test is_markov(TO.transfermatrix[1:(end-1), 2:end])
                    end
                end

                @testset "4D" begin
                    bins_visited_by_orbit = assign_bin_labels(E_4D, ϵ)
                    bininfo = get_binvisits(bins_visited_by_orbit)
                    TO = estimate_transferoperator_from_binvisits(bininfo)
                    @test typeof(TO) <: TransferOperatorRectangularBinning

                    # Last row might sum to zero, because the last point does not need to
            		# be contained in the last bin. However, the remaining row sums must
            		# be one.
                    if !is_markov(TO)
                        @warn "There were all-zero columns in the transfer matrix"
                        @warn "Removing first column and last row"
                        if is_markov(TO.transfermatrix[1:(end-1), 2:end])
                            info("That made the transfer matrix Markov")
                        else
                            @warn "That did NOT make the transfer matrix Markov"
                            #@save "from_binvisits_4D.jld2" TO, E_4D, bininfo
                        end
                        @test is_markov(TO.transfermatrix[1:(end-1), 2:end])
                    end
                end

            end

            @testset "From embedding" begin
                @testset "2D" begin
                    TO = TransferOperatorEstimatorRectangularBinning(E_2D, ϵ)
                    @test typeof(TO) <: TransferOperatorRectangularBinning
                    if !is_markov(TO)
                        @warn "There were all-zero columns in the transfer matrix"
                        @warn "Removing first column and last row"
                        if is_markov(TO.transfermatrix[1:(end-1), 2:end])
                            info("That made the transfer matrix Markov")
                        else
                            @warn "That did NOT make the transfer matrix Markov"
                            #@save "from_embedding_2D.jld2" TO, E_2D
                        end
                        @test is_markov(TO.transfermatrix[1:(end-1), 2:end])
                    end
                end

                @testset "3D" begin
                    TO = TransferOperatorEstimatorRectangularBinning(E_3D, ϵ)
                    @test typeof(TO) <: TransferOperatorRectangularBinning
                    if !is_markov(TO)
                        @warn "There were all-zero columns in the transfer matrix"
                        @warn "Removing first column and last row"
                        if is_markov(TO.transfermatrix[1:(end-1), 2:end])
                            info("That made the transfer matrix Markov")
                        else
                            @warn "That did NOT make the transfer matrix Markov"
                            #@save "from_embedding_3D.jld2" TO, E_3D
                        end
                        @test is_markov(TO.transfermatrix[1:(end-1), 2:end])
                    end
                end

                @testset "4D" begin
                    TO = TransferOperatorEstimatorRectangularBinning(E_4D, ϵ)
                    @test typeof(TO) <: TransferOperatorRectangularBinning
                    if !is_markov(TO)
                        @warn "There were all-zero columns in the transfer matrix"
                        @warn "Removing first column and last row"
                        if is_markov(TO.transfermatrix[1:(end-1), 2:end])
                            info("That made the transfer matrix Markov")
                        else
                            @warn "That did NOT make the transfer matrix Markov"
                            #@save "from_embedding_4D.jld2" TO, E_4D
                        end
                        @test is_markov(TO.transfermatrix[1:(end-1), 2:end])
                    end
                end
            end

            @testset "From points" begin
                @testset "2D" begin
                    TO = TransferOperatorEstimatorRectangularBinning(points_2D, ϵ)

                    @test typeof(TO) <: TransferOperatorRectangularBinning

                    if !is_markov(TO)
                        @warn "There were all-zero columns in the transfer matrix"
                        @warn "Removing first column and last row"
                        if is_markov(TO.transfermatrix[1:(end-1), 2:end])
                            info("That made the transfer matrix Markov")
                        else
                            @warn "That did NOT make the transfer matrix Markov"
                            #@save "from_points_2D.jld2" TO, points_2D
                        end
                        @test is_markov(TO.transfermatrix[1:(end-1), 2:end])
                    end
                end

                @testset "3D" begin
                    TO = TransferOperatorEstimatorRectangularBinning(points_3D, ϵ)

                    @test typeof(TO) <: TransferOperatorRectangularBinning

                    if !is_markov(TO)
                        @warn "There were all-zero columns in the transfer matrix"
                        @warn "Removing first column and last row"
                        if is_markov(TO.transfermatrix[1:(end-1), 2:end])
                            info("That made the transfer matrix Markov")
                        else
                            @warn "That did NOT make the transfer matrix Markov"
                            #@save "from_points_3D.jld2" TO, points_3D
                        end
                        @test is_markov(TO.transfermatrix[1:(end-1), 2:end])
                    end
                end

                @testset "4D" begin
                    TO = TransferOperatorEstimatorRectangularBinning(points_4D, ϵ)
                    @test typeof(TO) <: TransferOperatorRectangularBinning
                    if !is_markov(TO)
                        @warn "There were all-zero columns in the transfer matrix"
                        @warn "Removing first column and last row"
                        if is_markov(TO.transfermatrix[1:(end-1), 2:end])
                            info("That made the transfer matrix Markov")
                        else
                            @warn "That did NOT make the transfer matrix Markov"
                            #@save "from_points_4D.jld2" TO, points_4D
                        end
                        @test is_markov(TO.transfermatrix[1:(end-1), 2:end])
                    end
                end
            end
        end
    end
end


text:
function ask(q)
    println(q)
    readline()
end

function onerun(setups, wind_arduinos, led_arduino, camera)
    md = Dict()
    labels = setups.label
    @label start
    l = ask("Which setup?")
    i = findfirst(==(l), labels)
    if isnothing(i)
        @info "I cannot find your choice, $l, in the available setups" labels
        @goto start
    end
    md["setup"] = label_setup(setups[i,:])
    for _ in 1:3
        update_arena!(wind_arduinos, led_arduino, setups[i,:])
    end
    ask("Press Enter to start recording")
    md["recording_time"] = string(now())
    folder = datadir / md["recording_time"]
    mkdir(folder)
    @async record(camera, wind_arduinos, folder)
    ask("Press Enter to stop recording")
    close(camera)
    md["beetle_id"]  = ask("Beetle ID")
    md["comment"] = ask("Comments")
    res = ask("Save? [Y]/n")
    if  isempty(res) || occursin(r"^y"i, res)
        open(folder / "metadata.toml", "w") do io
            TOML.print(io, md)
        end
    else
        rm(folder, force=true, recursive=true)
    end
end

function main(; setup_file = HTTP.get(setupsurl).body, fan_ports = ["/dev/serial/by-id/usb-Arduino__www.arduino.cc__0043_957353530323510141D0-if00", "/dev/serial/by-id/usb-Arduino__www.arduino.cc__0043_95635333930351917172-if00", "/dev/serial/by-id/usb-Arduino__www.arduino.cc__0043_95735353032351010260-if00", "/dev/serial/by-id/usb-Arduino__www.arduino.cc__0043_55838323435351213041-if00", "/dev/serial/by-id/usb-Arduino__www.arduino.cc__0043_957353530323514121D0-if00"], led_port = "/dev/serial/by-id/usb-1a86_USB2.0-Serial-if00-port0", cam_url = "http://130.235.245.118:8080/stream/video.mjpeg")

    wind_arduinos = [FanArduino(id, port) for (id, port) in enumerate(fan_ports) if isconnected(port)]
    led_arduino = LEDArduino(led_port)
    camera = PiCam(cam_url)

    df = CSV.File(setup_file, header = 1:2) |> DataFrame
    setups = select(df, :setup_label => identity => :label, r"fan" => ByRow(parse2wind ∘ tuple) => :fans, r"star" => ByRow(parse2stars ∘ tuple) => :stars)

    while true
        try
            onerun(setups, wind_arduinos, led_arduino, camera)
        catch ex
            res = ask("Quit? [Y]/n")
            if  isempty(res) || occursin(r"^y"i, res)
                close(camera)
                close(led_arduino)
                close.(wind_arduinos)
                throw(ex)
                return nothing
            else
                continue
            end
        end
    end
end

function record(camera, wind_arduinos, folder)
    open(camera)

    fan_io = open(folder / "fans.csv", "w")
    println(fan_io, "time,", join([join(["fan$(a.id)_speed$j" for j in 1:3], ",") for a in wind_arduinos], ","))

    tmp = folder / "temp.stream"

    open(tmp, "w") do stream_io
        i = 0
        while isopen(camera)
            i += 1

            img = get_frame(camera)
            appendencode!(camera.encoder, stream_io, img, i)

            t, rpms = get_rpms(wind_arduinos)
            println(fan_io, t, ",",join(Iterators.flatten(rpms), ","))

        end
        finishencode!(camera.encoder, stream_io)
    end
    close(fan_io)
    camera.encoder = prepareencoder(camera.buff; framerate, AVCodecContextProperties, codec_name)
end

function update_arena!(wind_arduinos, led_arduino, setup)
    for a in wind_arduinos
        a.pwm[] = setup.fans[a.id].pwm
    end
    led_arduino.msg[] = parse2arduino(setup.stars)
end

function backup()
    todo = readpath(datadir)
    n = length(todo)
    done = Vector{SystemPath}(undef, n)
    @showprogress 1 "Uploading..." for (i, folder) in enumerate(todo)
        tmp = folder / "temp.stream"
        video = folder / "track.mp4"
        mux(tmp, video, framerate, silent = true)

        tb = Tar.create(string(folder))
        source = AbstractPath(tb)
        name = basename(folder)
        destination = S3Path(bucket, name * ".tar", config = s3config)
        mv(source, destination)
        @assert AWSS3.s3_exists(s3config, "dackebeetle", name * ".tar") "upload failed for $name"
        done[i] = folder
    end
    foreach(done) do folder
        rm(folder, recursive = true)
    end
end

label_setup(x) = string("fans=", Int[i.pwm for i in x.fans], "; stars=", [string(i.cardinality, " ", i.elevation, " ", i.intensity, " ", i.radius) for i in x.stars])


text:
using ImageMagick, Base.Test

img = load(joinpath("images", "bad_exif_orientation.jpg"))
@test size(img) == (512,512)


text:
export PrioritizedDQNLearner

using Setfield: @set
using Random: AbstractRNG, GLOBAL_RNG
import Functors
using LinearAlgebra: dot

Base.@kwdef mutable struct PrioritizedDQNLearner{A<:Approximator{<:TwinNetwork}} <: AbstractLearner
    approximator::A
    loss_func::Any  # !!! here the loss func must return the loss before reducing over the batch dimension
    n::Int = 1
    γ::Float32 = 0.99f0
    β_priority::Float32 = 0.5f0
    is_enable_double_DQN::Bool = true
    rng::AbstractRNG = GLOBAL_RNG
    # for logging
    loss::Float32 = 0.0f0
end

(L::PrioritizedDQNLearner)(s::AbstractArray) = L.approximator(s)

Functors.functor(x::PrioritizedDQNLearner) = (; approximator=x.approximator), y -> @set x.approximator = y.approximator

function RLBase.optimise!(
    learner::PrioritizedDQNLearner,
    batch::Union{
        NamedTuple{(:key, :priority, SS′ART...)},
        NamedTuple{(:key, :priority, SS′L′ART...)}
    }
)
    A = learner.approximator
    Q = A.model.source
    Qₜ = A.model.target
    γ = learner.γ
    β = learner.β_priority
    loss_func = learner.loss_func
    n = learner.n

    s, s′, a, r, t = map(x -> batch[x], SS′ART)
    batch_size = length(a)
    a = CartesianIndex.(a, 1:batch_size)
    k, p = batch.key, batch.priority
    p′ = similar(p)

    w = 1.0f0 ./ ((p .+ 1.0f-10) .^ β)
    w ./= maximum(w)

    q′ = learner.is_enable_double_DQN ? Q(s′) : Qₜ(s′)

    if haskey(batch, :next_legal_actions_mask)
        q′ .+= ifelse.(batch[:next_legal_actions_mask], 0.0f0, typemin(Float32))
    end

    q′ₐ = learner.is_enable_double_DQN ? Qₜ(s′)[dropdims(argmax(q′, dims=1), dims=1)] : dropdims(maximum(q′; dims=1), dims=1)

    G = r .+ γ^n .* (1 .- t) .* q′ₐ

    gs = gradient(params(A)) do
        qₐ = Q(s)[a]
        batch_losses = loss_func(G, qₐ)
        loss = dot(vec(w), vec(batch_losses)) * 1 // batch_size
        ignore() do
            p′ .= vec((batch_losses .+ 1.0f-10) .^ β)
            learner.loss = loss
        end
        loss
    end

    optimise!(A, gs)
    k => p′
end

function RLBase.optimise!(policy::QBasedPolicy{<:PrioritizedDQNLearner}, trajectory::Trajectory)
    for batch in trajectory
        k, p = optimise!(policy, batch) |> send_to_host
        trajectory[:priority, k] = p
    end
end


text:
"""
	meanconfid(x;sigma)
Compute the mean with confidence interval assuming normal
distribution (or large sample size => central limit theorem)

**Input**
* x: vector of samples (NaNs will be removed)
* sigma: optional input in case sigma (standard deviation, not variance) is known
* alpha: (1-alpha) gives the conficende interval for x_mean

**Output**
* x_mean: estimated mean
* x_confid: tuple with (lower,upper) confidence bounds defined by alpha

**Example**
```
# known standard deviation
x = [2.,3.,5.,6.,9.]
x_mean,x_confid = meanconfid(x,sigma=2.5)

# including NaN and Unknown sigma
push!(x,NaN);
x_mean,x_confid = meanconfid(x)
```
"""
function meanconfid(x::Vector{Float64};alpha::Float64=0.05,sigma::Float64=NaN)
	x_use = HypoTest.prepdata(x);
	x_mean = Distributions.mean(x_use);
	s = isnan(sigma) ? Distributions.std(x_use) : sigma;
	if isnan(sigma) && length(x_use) < 100
		d = Distributions.TDist(length(x_use)-1);
	else
		d =  Distributions.Normal();
	end
	q = Distributions.quantile(d,1-alpha/2);
	return x_mean,(x_mean-q*s/sqrt(length(x_use)),x_mean+q*s/sqrt(length(x_use)))
end

"""
	stdconfid(x;x_mean)
Compute the standard deviation with confidence interval assuming normal
distribution

**Input**
* x: vector of samples (NaNs will be removed)
* x_mean: optional input in case true mean is known (=>n instead n-1)
* alpha: (1-alpha) gives the conficende interval for sigma

**Output**
* sigma: estimated standard deviation (not variance)
* x_confid: tuple with (lower,upper) confidence bounds defined by alpha

**Example**
```
x = [1.95,1.96,1.98,1.99,2.00,2.00,2.01,2.01,2.01,2.01,
	2.02,2.02,2.03,2.03,2.04,2.04,2.05,2.06,2.08,2.10];
# unknown mean value
s,s_interval = stdconfid(x,alpha=0.05)

# known mean
x = [19.01,19.02,18.99,19.00,19.05,19.05,19.0,18.98,18.99,19.0];
s,s_interval = stdconfid(x,alpha=0.05,x_mean=19.)
```
"""
function stdconfid(x::Vector{Float64};alpha::Float64=0.05,x_mean::Float64=NaN)
	x_use = HypoTest.prepdata(x);
	sigma = Distributions.std(x_use, mean = isnan(x_mean) ? nothing : x_mean,
								corrected = isnan(x_mean) ? true : false);
	k = isnan(x_mean) ? length(x_use)-1 : length(x_use);
	d = Distributions.Chisq(k);
	qlow,qhigh = Distributions.quantile(d,1-alpha/2),Distributions.quantile(d,alpha/2);
	return sigma,(sqrt(k*sigma^2/qlow), sqrt(k*sigma^2/qhigh))
end

max_stars_repo_name: emenems/HypoTest

max_issues_repo_name: emenems/HypoTest

max_forks_repo_name: emenems/HypoTest


text:
using NMEA
using Base.Test

nmeas = NMEAData()
f = open("testdata.txt", "r")
for line = readlines(f)
    mtype = parse_msg!(nmeas, line)
    if (mtype == "GGA")
        @test get(nmeas.last_GGA.valid)
    elseif (mtype == "RMC")
        @test get(nmeas.last_RMC.valid)
    elseif (mtype == "GSA")
        @test get(nmeas.last_GSA.valid)
    elseif (mtype == "GSV")
        @test get(nmeas.last_GSV.valid)
    elseif (mtype == "GBS")
        @test get(nmeas.last_GBS.valid)
    elseif (mtype == "VTG")
        @test get(nmeas.last_VTG.valid)
    elseif (mtype == "GLL")
        @test get(nmeas.last_GLL.valid)
    elseif (mtype == "ZDA")
        @test get(nmeas.last_ZDA.valid)
    elseif (mtype == "DTM")
        @test get(nmeas.last_DTM.valid)
    else
        continue
    end
end


text:
using Base.Test
using DPMeansClustering

@testset "DPMeansClustering" begin
	include("test_dp.jl")
end


text:
n = 30000 # sample size
x = sum(randn(n,3).^2,dims=2) # generates 10000 draws from χ²(3)
using StatPlots, Distributions
density(x, label="nonparametric density") # creates a nonparametric density plot from the sample
plot!(Chisq(3),label="χ²(3) density") # plots the true density on top


text:
using ChainRulesTestUtils, DifferenceEquations, Distributions, LinearAlgebra, Test, Zygote
using DelimitedFiles
using DiffEqBase
using FiniteDiff: finite_difference_gradient

# joint case
function joint_likelihood_2(A_0, A_1, A_2, B, C_0, C_1, C_2, u0, noise, observables, D; kwargs...)
    problem = QuadraticStateSpaceProblem(A_0, A_1, A_2, B, u0, (0, size(observables, 2)); C_0, C_1,
                                         C_2, observables_noise = D, noise, observables, kwargs...)
    return solve(problem).logpdf
end

# Matrices from RBC
A_0_rbc = [-7.824904812740593e-5, 0.0]
A_1_rbc = [0.9568351489231076 6.209371005755285; 3.0153731819288737e-18 0.20000000000000007]
A_2_rbc = cat([-0.00019761505863889124 0.03375055315837927; 0.0 0.0],
              [0.03375055315837913 3.128758481817603; 0.0 0.0]; dims = 3)
B_2_rbc = reshape([0.0; -0.01], 2, 1) # make sure B is a matrix
C_0_rbc = [7.824904812740593e-5, 0.0]
C_1_rbc = [0.09579643002426148 0.6746869652592109; 1.0 0.0]
C_2_rbc = cat([-0.00018554166974717046 0.0025652363153049716; 0.0 0.0],
              [0.002565236315304951 0.3132705036896446; 0.0 0.0]; dims = 3)
D_2_rbc = MvNormal(Diagonal(abs2.([0.1, 0.1])))
u0_2_rbc = zeros(2)

observables_2_rbc = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/RBC_observables.csv"),
                            ',')' |> collect
noise_2_rbc = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/RBC_noise.csv"), ',')' |>
              collect

# Data and Noise
T = 5
observables_2_rbc = observables_2_rbc[:, 1:T]
noise_2_rbc = noise_2_rbc[:, 1:T]

@testset "basic inference" begin
    prob = QuadraticStateSpaceProblem(A_0_rbc, A_1_rbc, A_2_rbc, B_2_rbc, u0_2_rbc,
                                      (0, size(observables_2_rbc, 2)); C_0 = C_0_rbc, C_1 = C_1_rbc,
                                      C_2 = C_2_rbc, observables_noise = D_2_rbc,
                                      noise = noise_2_rbc, observables = observables_2_rbc)
    @inferred QuadraticStateSpaceProblem(A_0_rbc, A_1_rbc, A_2_rbc, B_2_rbc, u0_2_rbc,
                                         (0, size(observables_2_rbc, 2)); C_0 = C_0_rbc,
                                         C_1 = C_1_rbc, C_2 = C_2_rbc, observables_noise = D_2_rbc,
                                         noise = noise_2_rbc, observables = observables_2_rbc)

    DiffEqBase.get_concrete_problem(prob, false)
    @inferred DiffEqBase.get_concrete_problem(prob, false)

    sol = solve(prob)
    @inferred solve(prob)
end

@testset "quadratic rbc joint likelihood" begin
    @test joint_likelihood_2(A_0_rbc, A_1_rbc, A_2_rbc, B_2_rbc, C_0_rbc, C_1_rbc, C_2_rbc,
                             u0_2_rbc, noise_2_rbc, observables_2_rbc, D_2_rbc) ≈ -690.81094364573
    @inferred joint_likelihood_2(A_0_rbc, A_1_rbc, A_2_rbc, B_2_rbc, C_0_rbc, C_1_rbc, C_2_rbc,
                                 u0_2_rbc, noise_2_rbc, observables_2_rbc, D_2_rbc) # would this catch inference problems in the solve?
    gradient((args...) -> joint_likelihood_2(args..., observables_2_rbc, D_2_rbc), A_0_rbc, A_1_rbc,
             A_2_rbc, B_2_rbc, C_0_rbc, C_1_rbc, C_2_rbc, u0_2_rbc, noise_2_rbc)
    test_rrule(Zygote.ZygoteRuleConfig(),
               (args...) -> joint_likelihood_2(args..., observables_2_rbc, D_2_rbc), A_0_rbc,
               A_1_rbc, A_2_rbc, B_2_rbc, C_0_rbc, C_1_rbc, C_2_rbc, u0_2_rbc, noise_2_rbc;
               rrule_f = rrule_via_ad, check_inferred = false)
end

# Load FVGQ data for checks
A_0_raw = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_A_0.csv"), ',')
A_0_FVGQ = vec(A_0_raw)
A_1_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_A_1.csv"), ',')
A_2_raw = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_A_2.csv"), ',')
A_2_FVGQ = reshape(A_2_raw, length(A_0_FVGQ), length(A_0_FVGQ), length(A_0_FVGQ))
B_2_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_B.csv"), ',')
C_0_raw = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_C_0.csv"), ',')
C_0_FVGQ = vec(C_0_raw)
C_1_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_C_1.csv"), ',')
C_2_raw = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_C_2.csv"), ',')
C_2_FVGQ = reshape(C_2_raw, length(C_0_FVGQ), length(A_0_FVGQ), length(A_0_FVGQ))
# D_raw = readdlm(joinpath(pkgdir(DifferenceEquations), "FVGQ_D.csv"); header = false)))
D_2_FVGQ = MvNormal(Diagonal(abs2.(ones(6) * 1e-3)))
u0_2_FVGQ = zeros(size(A_1_FVGQ, 1))

observables_2_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations),
                                      "test/data/FVGQ20_observables.csv"), ',')' |> collect

noise_2_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_noise.csv"), ',')' |>
               collect

@testset "quadratic FVGQ joint likelihood" begin
    @test joint_likelihood_2(A_0_FVGQ, A_1_FVGQ, A_2_FVGQ, B_2_FVGQ, C_0_FVGQ, C_1_FVGQ, C_2_FVGQ,
                             u0_2_FVGQ, noise_2_FVGQ, observables_2_FVGQ, D_2_FVGQ) ≈
          -1.473244794713955e10
    @inferred joint_likelihood_2(A_0_FVGQ, A_1_FVGQ, A_2_FVGQ, B_2_FVGQ, C_0_FVGQ, C_1_FVGQ,
                                 C_2_FVGQ, u0_2_FVGQ, noise_2_FVGQ, observables_2_FVGQ, D_2_FVGQ)
    gradient((args...) -> joint_likelihood_2(args..., observables_FVGQ, D_FVGQ), A_0_FVGQ, A_1_FVGQ,
             A_2_FVGQ, B_2_FVGQ, C_0_FVGQ, C_1_FVGQ, C_2_FVGQ, u0_2_FVGQ, noise_2_FVGQ)
    #test_rrule(Zygote.ZygoteRuleConfig(), (args...) -> joint_likelihood_2(args..., observables, D),
    #    A_0, A_1, A_2, B, C_0, C_1, C_2, u0, noise; rrule_f = rrule_via_ad,
    #    check_inferred = false)
end


text:
using Test
using ImageIO
using FileIO: File, DataFormat, Stream, @format_str

@testset "ImageIO" begin

    @testset "PNGs" begin
        img = rand(UInt8, 10, 10)
        f = File{DataFormat{:PNG}}("test_fpath.png")
        ImageIO.save(f, img)
        img_saveload = ImageIO.load(f)
        @test all(img .== reinterpret(UInt8, img_saveload))

        open(io->ImageIO.save(Stream(format"PNG", io), img), "test_io.png", "w")
        img_saveload = open(io->ImageIO.load(Stream(format"PNG", io)), "test_io.png")
        @test all(img .== reinterpret(UInt8, img_saveload))
    end

end

rm("test_fpath.png", force=true)
rm("test_io.png", force=true)


text:
# `permutedims` seems to work better with AD (cf. KernelFunctions.jl)
aT_b(a::AbstractVector{<:Real}, b::AbstractMatrix{<:Real}) = permutedims(a) * b
# `permutedims` can't be used here since scalar output is desired
aT_b(a::AbstractVector{<:Real}, b::AbstractVector{<:Real}) = dot(a, b)

# flatten arrays with fallback for scalars
_vec(x::AbstractArray{<:Real}) = vec(x)
_vec(x::Real) = x


text:
using EasyDataAugmentation
using Documenter

DocMeta.setdocmeta!(EasyDataAugmentation, :DocTestSetup, :(using EasyDataAugmentation); recursive=true)

makedocs(;
    modules=[EasyDataAugmentation],
    authors="lilianabs <lilianabsmath@google.com> and contributors",
    repo="https://github.com/lilianabs/EasyDataAugmentation.jl/blob/{commit}{path}#{line}",
    sitename="EasyDataAugmentation.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://lilianabs.github.io/EasyDataAugmentation.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/lilianabs/EasyDataAugmentation.jl",
    devbranch="main",
)


text:
module s
using Gtk.ShortNames

import eve
using DataFrames

immutable WinPoint
  x :: Clong
  y :: Clong
end
immutable WinMsg
  hwnd :: Ptr{GObject}
  msg :: Cuint
  wparam :: Ptr{Cuint}
  lparam :: Ptr{Void}
  time :: Culong
  pt :: WinPoint
end

function DataFrameView(df)
  store = @ListStore(eltypes(df)...)
  for r in eachrow(df)
    push!(store, (Array(r)...))
  end
  view = @TreeView(TreeModel(store))
  renderer = @CellRendererText()
  for i in eachindex(names(df))
    push!(view, @TreeViewColumn(string(names(df)[i]), renderer, Dict("text"=>i-1)))
  end
  view
end

function MarginToolView(data :: eve.ItemSummery)
  grid = @Grid()
  setproperty!(grid, "orientation", 1)
  push!(grid, @Label(string(data.itemID)))
  push!(grid, @Label("Max Buy: $(data.maxbuy)"))
  push!(grid, @Label("Min Sell: $(data.minsell)"))
  push!(grid, @Label("Margin: $(100data.margin)%"))
  grid
end
function HotkeyFilter(xevt :: Ptr{WinMsg}, gevt :: Ptr{Gtk.GdkEventKey}, data :: Ptr{Void})
  evt = unsafe_load(xevt)
  if evt.msg == 0x0312
    return Cint(2)
    GAccessor.
  else
    return Cint(0)
  end
end
function UITest()
  t = eve.importOrders()
  win = @Window("test", 1920, 1080)
  (sells, buys) = groupby(t, :bid)
  listViewSells = DataFrameView(sells)
  listViewBuys = DataFrameView(buys)

  f = @Frame("Orders")
  panes = @Paned(true)
  panes[1] = listViewSells
  panes[2] = listViewBuys
  push!(f, panes)
  margin = @Grid()
  dv = MarginDataView()
  push!(margin, dv)
  importOrdersBtn = @Button("Import Orders")
  push!(margin, importOrdersBtn)
  signal_connect(importOrdersBtn, "clicked") do widget
    UpdateMarginToolView(dv, eve.summerizeItem(eve.importMarketLog()))
  end
  tabs = @Notebook()
  push!(tabs, f, "Orders")
  push!(tabs, margin, "Margin Tool")

  push!(win, tabs)
  gdkWin = Gtk.gdk_window(win)
  ccall( (:gdk_window_add_filter, Gtk.libgdk), Void, (Ptr{Void}, Ptr{TreeView}, Ptr{Void}),
    gdkWin,
    cfunction(HotkeyFilter, Cint, (Ptr{WinMsg}, Ptr{Void}, Ptr{Void})),
    pointer_from_objref(listViewSells))
  hwnd = ccall( (:gdk_win32_window_get_impl_hwnd, Gtk.libgdk),  Ptr{Void}, (Ptr{Void},), gdkWin)
  result = ccall( (:RegisterHotKey, "User32"), stdcall, Cint, (Ptr{Void}, Cint, Cuint, Cuint),
    hwnd, 1, 0, 0x28)
  #sel = getproperty(listViewSells, :selection, TreeSelection)
  sel = Gtk.GAccessor.selection(listViewSells)
  signal_connect(sel, "changed") do selection
    print(selected(selection))
  end
  showall(win)
end
end


text:
using Enzyme
using Test

@testset "ABI & Calling convention" begin

    f(x) = x

    # GhostType -> Nothing
    res = autodiff(Reverse, f, Const, Const(nothing))
    @test res === ()

    @test () === autodiff(Forward, f, Const, Const(nothing))

    res = autodiff(Reverse, f, Const(nothing))
    @test res === ()

    @test () === autodiff(Forward, f, Const(nothing))

    res = Enzyme.autodiff_deferred(f, Const(nothing))
    @test res === ()
    @test () === Enzyme.fwddiff_deferred(f, Const(nothing))

    # ConstType -> Type{Int}
    res = autodiff(Reverse, f, Const, Const(Int))
    @test res === ()
    @test () === autodiff(Forward, f, Const, Const(Int))

    res = autodiff(Reverse, f, Const(Int))
    @test res === ()
    @test () === autodiff(Forward, f, Const(Int))

    res = Enzyme.autodiff_deferred(f, Const(Int))
    @test res === ()
    @test () === Enzyme.fwddiff_deferred(f, Const(Int))

    # Complex numbers
    cres,  = Enzyme.autodiff(Reverse, f, Active, Active(1.5 + 0.7im))
    @test cres ≈ 1.0 + 0.0im
    cres,  = Enzyme.autodiff(Forward, f, DuplicatedNoNeed, Duplicated(1.5 + 0.7im, 1.0 + 0im))
    @test cres ≈ 1.0 + 0.0im

    cres,  = Enzyme.autodiff(Reverse, f, Active(1.5 + 0.7im))
    @test cres ≈ 1.0 + 0.0im
    cres,  = Enzyme.autodiff(Forward, f, Duplicated(1.5 + 0.7im, 1.0+0im))
    @test cres ≈ 1.0 + 0.0im

    cres, = Enzyme.autodiff_deferred(f, Active(1.5 + 0.7im))
    @test cres ≈ 1.0 + 0.0im
    cres,  = Enzyme.fwddiff_deferred(f, Duplicated(1.5 + 0.7im, 1.0+0im))
    @test cres ≈ 1.0 + 0.0im

    # Unused singleton argument
    unused(_, y) = y
    res0, = autodiff(Reverse, unused, Active, Const(nothing), Active(2.0))
    @test res0 ≈ 1.0
    res0, = autodiff(Forward, unused, DuplicatedNoNeed, Const(nothing), Duplicated(2.0, 1.0))
    @test res0 ≈ 1.0
    res0, = autodiff(Forward, unused, DuplicatedNoNeed, Const(nothing), DuplicatedNoNeed(2.0, 1.0))
    @test res0 ≈ 1.0

    res0, = autodiff(Reverse, unused, Const(nothing), Active(2.0))
    @test res0 ≈ 1.0
    res0, = autodiff(Forward, unused, Const(nothing), Duplicated(2.9, 1.0))
    @test res0 ≈ 1.0

    # returning an Array, with Const
    function squareRetArray(x)
        x[1] *= 2
        x
    end
    x = [0.0]
    dx = [1.2]
    autodiff(Reverse, squareRetArray, Const, Duplicated(x, dx))
    @test dx[1] ≈ 2.4

    dx = [1.2]
    @test () === autodiff(Forward, squareRetArray, Const, Duplicated(x, dx))
    @test dx[1] ≈ 2.4

    x = [0.0]
    dx = [1.2]
    Enzyme.autodiff_deferred(squareRetArray, Const, Duplicated(x, dx))

    dx = [1.2]
    @test () === autodiff(Forward, squareRetArray, Const, Duplicated(x, dx))
    @test dx[1] ≈ 2.4

    # Multi arg => sret
    mul(x, y) = x * y
    pair = autodiff(Reverse, mul, Active, Active(2.0), Active(3.0))
    @test pair[1] ≈ 3.0
    @test pair[2] ≈ 2.0

    pair = autodiff(Reverse, mul, Active(2.0), Active(3.0))
    @test pair[1] ≈ 3.0
    @test pair[2] ≈ 2.0

    pair = Enzyme.autodiff_deferred(mul, Active(2.0), Active(3.0))
    @test pair[1] ≈ 3.0
    @test pair[2] ≈ 2.0

    # Multi output
    # TODO broken arg convention?
    # tup(x) = (x, x*2)
    # pair = first(autodiff(Forward, tup, DuplicatedNoNeed, Duplicated(3.14, 1.0)))
    # @test pair[1] ≈ 1.0
    # @test pair[2] ≈ 2.0
    # pair = first(autodiff(Forward, tup, Duplicated(3.14, 1.0)))
    # @test pair[1] ≈ 1.0
    # @test pair[2] ≈ 2.0
    # pair = first(Enzyme.fwddiff_deferred(tup, Duplicated(3.14, 1.0)))
    # @test pair[1] ≈ 1.0
    # @test pair[2] ≈ 2.0

    # SequentialType
    struct Foo
        baz::Int
        qux::Float64
    end

    g(x) = x.qux
    res2,  = autodiff(Reverse, g, Active, Active(Foo(3, 1.2)))
    @test res2.qux ≈ 1.0

    @test 1.0≈ first(autodiff(Forward, g, DuplicatedNoNeed, Duplicated(Foo(3, 1.2), Foo(0, 1.0))))

    res2,  = autodiff(Reverse, g, Active(Foo(3, 1.2)))
    @test res2.qux ≈ 1.0

    @test 1.0≈ first(autodiff(Forward, g, Duplicated(Foo(3, 1.2), Foo(0, 1.0))))

    unused2(_, y) = y.qux
    resF, = autodiff(Reverse, unused2, Active, Const(nothing), Active(Foo(3, 2.0)))
    @test resF.qux ≈ 1.0

    @test 1.0≈ first(autodiff(Forward, unused2, DuplicatedNoNeed, Const(nothing), Duplicated(Foo(3, 1.2), Foo(0, 1.0))))

    resF, = autodiff(Reverse, unused2, Const(nothing), Active(Foo(3, 2.0)))
    @test resF.qux ≈ 1.0

    @test 1.0≈ first(autodiff(Forward, unused2, Const(nothing), Duplicated(Foo(3, 1.2), Foo(0, 1.0))))

    h(x, y) = x.qux * y.qux
    res3 = autodiff(Reverse, h, Active, Active(Foo(3, 1.2)), Active(Foo(5, 3.4)))
    @test res3[1].qux ≈ 3.4
    @test res3[2].qux ≈ 1.2

    @test 7*3.4 + 9 * 1.2 ≈ first(autodiff(Forward, h, DuplicatedNoNeed, Duplicated(Foo(3, 1.2), Foo(0, 7.0)), Duplicated(Foo(5, 3.4), Foo(0, 9.0))))

    res3 = autodiff(Reverse, h, Active(Foo(3, 1.2)), Active(Foo(5, 3.4)))
    @test res3[1].qux ≈ 3.4
    @test res3[2].qux ≈ 1.2

    @test 7*3.4 + 9 * 1.2 ≈ first(autodiff(Forward, h, Duplicated(Foo(3, 1.2), Foo(0, 7.0)), Duplicated(Foo(5, 3.4), Foo(0, 9.0))))

    caller(f, x) = f(x)
    res4, = autodiff(Reverse, caller, Active, (x)->x, Active(3.0))
    @test res4 ≈ 1.0

    res4, = autodiff(Forward, caller, DuplicatedNoNeed, (x)->x, Duplicated(3.0, 1.0))
    @test res4 ≈ 1.0

    res4, = autodiff(Reverse, caller, (x)->x, Active(3.0))
    @test res4 ≈ 1.0

    res4, = autodiff(Forward, caller, (x)->x, Duplicated(3.0, 1.0))
    @test res4 ≈ 1.0

    struct LList
        next::Union{LList,Nothing}
        val::Float64
    end

    function sumlist(n::LList)
        sum = 0.0
        while n !== nothing
            sum += n.val
            n = n.next
        end
        sum
    end

    regular = LList(LList(nothing, 1.0), 2.0)
    shadow  = LList(LList(nothing, 0.0), 0.0)
    ad = autodiff(Reverse, sumlist, Active, Duplicated(regular, shadow))
    @test ad === ()
    @test shadow.val ≈ 1.0 && shadow.next.val ≈ 1.0

    @test 2.0 ≈ first(autodiff(Forward, sumlist, DuplicatedNoNeed, Duplicated(regular, shadow)))

    mulr(x, y) = x[] * y[]
    x = Ref(2.0)
    y = Ref(3.0)
    dx = Ref(0.0)
    dy = Ref(0.0)
    n = autodiff(Reverse, mulr, Active, Duplicated(x, dx), Duplicated(y, dy))
    @test n === ()
    @test dx[] ≈ 3.0
    @test dy[] ≈ 2.0

    x = Ref(2.0)
    y = Ref(3.0)
    dx = Ref(5.0)
    dy = Ref(7.0)
    @test 5.0*3.0 + 2.0*7.0≈ first(autodiff(Forward, mulr, DuplicatedNoNeed, Duplicated(x, dx), Duplicated(y, dy)))

    mid, = Enzyme.autodiff(Reverse, (fs, x) -> fs[1](x), Active, (x->x*x,), Active(2.0))
    @test mid ≈ 4.0

    mid, = Enzyme.autodiff(Reverse, (fs, x) -> fs[1](x), Active, [x->x*x], Active(2.0))
    @test mid ≈ 4.0

    mid, = Enzyme.autodiff(Forward, (fs, x) -> fs[1](x), DuplicatedNoNeed, (x->x*x,), Duplicated(2.0, 1.0))
    @test mid ≈ 4.0

    mid, = Enzyme.autodiff(Forward, (fs, x) -> fs[1](x), DuplicatedNoNeed, [x->x*x], Duplicated(2.0, 1.0))
    @test mid ≈ 4.0

    # deserves_argbox yes and no
    struct Bar
        r::Ref{Int}
    end

    # ConstType

    # primitive type Int128, Float64, Float128

    # returns: sret, const/ghost, !deserve_retbox
end

@testset "Mutable Struct ABI" begin
    mutable struct MStruct
        val::Float32
    end

    function sqMStruct(domain::Vector{MStruct}, x::Float32)
       @inbounds domain[1] = MStruct(x*x)
       return nothing
    end

    orig   = [MStruct(0.0)]
    shadow = [MStruct(17.0)]
    Enzyme.autodiff(Forward, sqMStruct, Duplicated(orig, shadow), Duplicated(Float32(3.14), Float32(1.0)))
     @test 2.0*3.14 ≈ shadow[1].val

end

@testset "Thread helpers" begin
    function timesID(x)
        x * Threads.threadid()
    end
    function timesNum(x)
        x * Threads.nthreads()
    end
    @test Threads.threadid() ≈ Enzyme.autodiff(Reverse, timesID, Active(2.0))[1]
    @test Threads.nthreads() ≈ Enzyme.autodiff(Reverse, timesNum, Active(2.0))[1]
end

@testset "Closure ABI" begin
    function clo2(x)
        V = [x]
        y -> V[1] * y
    end
    f = clo2(2.0)
    @test 2.0 ≈ Enzyme.autodiff(Reverse, f, Active(3.0))[1]

    @test 2.0 ≈ Enzyme.autodiff(Forward, f, Duplicated(3.0, 1.0))[1]

    df = clo2(0.0)
    @test 2.0 ≈ Enzyme.autodiff(Reverse, Duplicated(f, df), Active(3.0))[1]
    @test 3.0 ≈ df.V[1]

    @test 2.0 * 7.0 + 3.0 * 5.0 ≈ first(Enzyme.autodiff(Forward, Duplicated(f, df), Duplicated(5.0, 7.0)))
end

@testset "Callable ABI" begin
    function method(f, x)
        return f(x)
    end

    struct AFoo
       x::Float64
    end

    function (f::AFoo)(x::Float64)
       return f.x * x
    end

    @test Enzyme.autodiff(Reverse, method, Active, AFoo(2.0), Active(3.0))[1]≈ 2.0
    @test Enzyme.autodiff(Reverse, AFoo(2.0), Active, Active(3.0))[1]≈ 2.0

    @test Enzyme.autodiff(Forward, method, DuplicatedNoNeed, AFoo(2.0), Duplicated(3.0, 1.0))[1]≈ 2.0
    @test Enzyme.autodiff(Forward, AFoo(2.0), DuplicatedNoNeed, Duplicated(3.0, 1.0))[1]≈ 2.0

    struct ABar
    end

    function (f::ABar)(x::Float64)
       return 2.0 * x
    end

    @test Enzyme.autodiff(Reverse, method, Active, ABar(), Active(3.0))[1]≈ 2.0
    @test Enzyme.autodiff(Reverse, ABar(), Active, Active(3.0))[1]≈ 2.0

    @test Enzyme.autodiff(Forward, method, DuplicatedNoNeed, ABar(), Duplicated(3.0, 1.0))[1]≈ 2.0
    @test Enzyme.autodiff(Forward, ABar(), DuplicatedNoNeed, Duplicated(3.0, 1.0))[1]≈ 2.0
end


text:
using LogCompose, Test, Logging

function test()
    config = joinpath(@__DIR__, "testapp.toml")
    simple_logfile = "simple.log"
    rm(simple_logfile; force=true)

    let logger = LogCompose.logger(config, "simple"; section="loggers")
        with_logger(logger) do
            @info("testsimple")
        end
        flush(logger.stream)

        log_file_contents = readlines(simple_logfile)
        @test findfirst("testsimple", log_file_contents[1]) !== nothing
    end

    let pipe = Pipe()
        Base.link_pipe!(pipe)
        writer = Base.pipe_writer(pipe)
        redirect_stdout(writer) do
            let logger = LogCompose.logger(config, "console"; section="loggers")
                with_logger(logger) do
                    @info("testconsole")
                end
                flush(logger.stream)
            end
        end
        close(writer)
        log_file_contents = readlines(Base.pipe_reader(pipe))
        @test findfirst("testconsole", log_file_contents[1]) !== nothing
        close(pipe)
    end

    let logger = LogCompose.logger(config, "null"; section="loggers")
        with_logger(logger) do
            @info("testnull")
        end
    end

    invalid_config = Dict{String,Any}(
        "invalid" => Dict{String,Any}(
            "type" => "UnknownModule.UnknownLogger"
        )
    )
    @test_throws ErrorException LogCompose.logger(invalid_config, "invalid")
    @test_throws ErrorException LogCompose.logcompose(String, invalid_config, invalid_config)

    try
        rm(simple_logfile; force=true)
    catch ex
        # ignore (occasionally fails with resource busy exception on Windows, because logger has not been gc'd yet)
    end
end

test()


text:
# To be executed from CAGD.jl

todisplay = VERSION <= VersionNumber("1.2") ? true : false

using CAGD
using SparseArrays
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
if todisplay  using ViewerGL; GL = ViewerGL  end

include("./CAGD.jl/paper/examples/cellular_complex.jl");
include("./CAGD.jl/paper/examples/CSG_example.jl");
include("./CAGD.jl/paper/examples/cycles_boundaries.jl");
include("./CAGD.jl/paper/examples/planar_algebra_2D.jl");
include("./CAGD.jl/paper/examples/solid_algebra_3D.jl");
include("./CAGD.jl/paper/examples/space_arrangement_2d.jl");
include("./CAGD.jl/paper/examples/space_arrangement.jl");
include("./CAGD.jl/paper/examples/subdivision_of_2_cells.jl");


text:
function HMMviterbi(v,phghm,ph1,pvgh)
    #HMMVITERBI Viterbi most likely joint hidden state of a HMM
    # [maxstate logprob]=HMMviterbi(v,phghm,ph1,pvgh)
    #
    # Inputs:
    # v : visible (obervation) sequence being a vector v=[2 1 3 3 1 ...]
    # phghm : homogeneous transition distribution phghm(i,j)=p(h(t)=i|h(t-1)=j)
    # ph1 : initial distribution
    # pvgh : time-independent emission distribution pvgh(i,j)=p(v(t)=i|h(t)=j)
    #
    # Outputs:
    # maxstate : most likely joint hidden (latent) state sequence
    # logprob : associated log probability of the most likely hidden sequence

    T=length(v); H=size(phghm,1)
    mu=zeros(H,T)
    mu[:,T]=ones(H,1)
    hs=zeros(Int64,1,T)
    pvghtrans=pvgh'
    for t=T:-1:2
        tmp = repmat(pvghtrans[:,v[t]].*mu[:,t],1,H).*phghm
        tmp = maximum(tmp,1)'
        mu[:,t-1]= tmp./sum(tmp) # normalise to avoid underflow
    end

    # backtrack:
    val, hs[1]=findmax(ph1.*pvghtrans[:,v[1]].*mu[:,1])
    for t=2:T
        tmp = pvghtrans[:,v[t]].*phghm[:,hs[t-1]]
        val, hs[t]=findmax(tmp.*mu[:,t])
    end

    logprob=log(ph1[hs[1]])+log(pvgh[v[1],hs[1]])
    for t=2:T
        logprob=logprob+log(phghm[hs[t],hs[t-1]])+log(pvgh[v[t],hs[t]])
    end
    return hs, logprob
end


text:
function  createcgcells(cgelcon,tlocal,ne)

nce, nve = size(tlocal);
cells = zeros(Int,ne*nce,nve);
for el=1:ne
    ind = (nce*(el-1)+1):nce*el;
    cells[ind,:] = reshape(cgelcon[tlocal[:],el],(nce,nve));
end
cells = cells .- 1;

return cells

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
###############
# API methods #
###############

"""
    ForwardDiff.gradient(f, x::AbstractArray, cfg::GradientConfig = GradientConfig(f, x), check=Val{true}())

Return `∇f` evaluated at `x`, assuming `f` is called as `f(x)`.

This method assumes that `isa(f(x), Real)`.

Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care.
"""
function gradient(f, x::AbstractArray, cfg::GradientConfig{T} = GradientConfig(f, x), ::Val{CHK}=Val{true}()) where {T, CHK}
    CHK && checktag(T, f, x)
    if chunksize(cfg) == length(x)
        return vector_mode_gradient(f, x, cfg)
    else
        return chunk_mode_gradient(f, x, cfg)
    end
end

"""
    ForwardDiff.gradient!(result::Union{AbstractArray,DiffResult}, f, x::AbstractArray, cfg::GradientConfig = GradientConfig(f, x), check=Val{true}())

Compute `∇f` evaluated at `x` and store the result(s) in `result`, assuming `f` is called as
`f(x)`.

This method assumes that `isa(f(x), Real)`.

"""
function gradient!(result::Union{AbstractArray,DiffResult}, f::F, x::AbstractArray, cfg::GradientConfig{T} = GradientConfig(f, x), ::Val{CHK}=Val{true}()) where {T, CHK, F}
    CHK && checktag(T, f, x)
    if chunksize(cfg) == length(x)
        vector_mode_gradient!(result, f, x, cfg)
    else
        chunk_mode_gradient!(result, f, x, cfg)
    end
    return result
end

@inline gradient(f, x::SArray)                      = vector_mode_gradient(f, x)
@inline gradient(f, x::SArray, cfg::GradientConfig) = gradient(f, x)
@inline gradient(f, x::SArray, cfg::GradientConfig, ::Val) = gradient(f, x)

@inline gradient!(result::Union{AbstractArray,DiffResult}, f, x::SArray) = vector_mode_gradient!(result, f, x)
@inline gradient!(result::Union{AbstractArray,DiffResult}, f, x::SArray, cfg::GradientConfig) = gradient!(result, f, x)
@inline gradient!(result::Union{AbstractArray,DiffResult}, f, x::SArray, cfg::GradientConfig, ::Val) = gradient!(result, f, x)

#####################
# result extraction #
#####################

@generated function extract_gradient(::Type{T}, y::Real, ::SArray{S,X,D,N}) where {T,S,X,D,N}
    result = Expr(:tuple, [:(partials(T, y, $i)) for i in 1:N]...)
    return quote
        $(Expr(:meta, :inline))
        return SArray{S}($result)
    end
end

function extract_gradient!(::Type{T}, result::DiffResult, y::Real) where {T}
    result = DiffResults.value!(result, y)
    grad = DiffResults.gradient(result)
    fill!(grad, zero(y))
    return result
end

function extract_gradient!(::Type{T}, result::DiffResult, dual::Dual) where {T}
    result = DiffResults.value!(result, value(T, dual))
    result = DiffResults.gradient!(result, partials(T, dual))
    return result
end

extract_gradient!(::Type{T}, result::AbstractArray, y::Real) where {T} = fill!(result, zero(y))
extract_gradient!(::Type{T}, result::AbstractArray, dual::Dual) where {T}= copyto!(result, partials(T, dual))

function extract_gradient_chunk!(::Type{T}, result, dual, index, chunksize) where {T}
    offset = index - 1
    for i in 1:chunksize
        result[i + offset] = partials(T, dual, i)
    end
    return result
end

function extract_gradient_chunk!(::Type{T}, result::DiffResult, dual, index, chunksize) where {T}
    extract_gradient_chunk!(T, DiffResults.gradient(result), dual, index, chunksize)
    return result
end

###############
# vector mode #
###############

function vector_mode_gradient(f::F, x, cfg::GradientConfig{T}) where {T, F}
    ydual = vector_mode_dual_eval(f, x, cfg)
    result = similar(x, valtype(ydual))
    return extract_gradient!(T, result, ydual)
end

function vector_mode_gradient!(result, f::F, x, cfg::GradientConfig{T}) where {T, F}
    ydual = vector_mode_dual_eval(f, x, cfg)
    result = extract_gradient!(T, result, ydual)
    return result
end

@inline function vector_mode_gradient(f::F, x::SArray{S,V}) where {F,S,V}
    T = typeof(Tag(f,V))
    return extract_gradient(T, static_dual_eval(T, f, x), x)
end

@inline function vector_mode_gradient!(result, f::F, x::SArray{S,V}) where {F,S,V}
    T = typeof(Tag(f,V))
    return extract_gradient!(T, result, static_dual_eval(T, f, x))
end

##############
# chunk mode #
##############

function chunk_mode_gradient_expr(result_definition::Expr)
    return quote
        @assert length(x) >= N "chunk size cannot be greater than length(x) ($(N) > $(length(x)))"

        # precalculate loop bounds
        xlen = length(x)
        remainder = xlen % N
        lastchunksize = ifelse(remainder == 0, N, remainder)
        lastchunkindex = xlen - lastchunksize + 1
        middlechunks = 2:div(xlen - lastchunksize, N)

        # seed work vectors
        xdual = cfg.duals
        seeds = cfg.seeds
        seed!(xdual, x)

        # do first chunk manually to calculate output type
        seed!(xdual, x, 1, seeds)
        ydual = f(xdual)
        $(result_definition)
        extract_gradient_chunk!(T, result, ydual, 1, N)
        seed!(xdual, x, 1)

        # do middle chunks
        for c in middlechunks
            i = ((c - 1) * N + 1)
            seed!(xdual, x, i, seeds)
            ydual = f(xdual)
            extract_gradient_chunk!(T, result, ydual, i, N)
            seed!(xdual, x, i)
        end

        # do final chunk
        seed!(xdual, x, lastchunkindex, seeds, lastchunksize)
        ydual = f(xdual)
        extract_gradient_chunk!(T, result, ydual, lastchunkindex, lastchunksize)

        # get the value, this is a no-op unless result is a DiffResult
        extract_value!(T, result, ydual)

        return result
    end
end

@eval function chunk_mode_gradient(f::F, x, cfg::GradientConfig{T,V,N}) where {F,T,V,N}
    $(chunk_mode_gradient_expr(:(result = similar(x, valtype(ydual)))))
end

@eval function chunk_mode_gradient!(result, f::F, x, cfg::GradientConfig{T,V,N}) where {F,T,V,N}
    $(chunk_mode_gradient_expr(:()))
end


text:
export BufferActor, buffer, getvalues

"""
    BufferActor{D}() where D

Buffer actor provides a storage actor. It copies last incoming successful `next` events in a `values` array, throws an ErrorException on `error!` event
and does nothing on completion event. Note: Actor does not check the size of incoming data.

# Examples
```jldoctest
using Rocket

source = of([ 1, 2, 3 ])
actor  = buffer(Int, 3)

subscribe!(source, actor)
show(getvalues(actor))

# output
[1, 2, 3]
```

See also: [`Actor`](@ref), [`buffer`](@ref)
"""
struct BufferActor{T, R} <: Actor{R}
    values :: R
end

function BufferActor(::Type{T}, size) where T
    storage = Array{T}(undef, size)
    return BufferActor{T, typeof(storage)}(storage)
end

getvalues(actor::BufferActor) = actor.values

on_next!(actor::BufferActor{T, R}, data::R) where {T, R} = copyto!(actor.values, data)
on_error!(actor::BufferActor, err)                       = error(err)
on_complete!(actor::BufferActor)                         = begin end

"""
    buffer(::Type{T}, size::Int) where T

# Arguments
- `::Type{T}`: Type of data in buffer
- `size::Int`: size of buffer

Creation operator for the `BufferActor` actor.

# Examples

```jldoctest
using Rocket

actor = buffer(Int, 3)
actor isa BufferActor{Int}

# output
true
```

See also: [`BufferActor`](@ref), [`AbstractActor`](@ref)
"""
buffer(::Type{T}, size...)     where T = BufferActor(T, size)
buffer(::Type{T}, size::Tuple) where T = BufferActor(T, size)

# Julia iterable interface

Base.IteratorSize(::Type{ <: BufferActor{T, R} })   where { T, R } = Base.IteratorSize(R)
Base.IteratorEltype(::Type{ <: BufferActor{T, R} }) where { T, R } = Base.IteratorEltype(R)

Base.IndexStyle(::Type{ <: BufferActor{T, R} }) where { T, R } = Base.IndexStyle(R)

Base.eltype(::Type{ <: BufferActor{T} }) where T = T

Base.iterate(actor::BufferActor)        = iterate(actor.values)
Base.iterate(actor::BufferActor, state) = iterate(actor.values, state)

Base.collect(actor::BufferActor)     = collect(actor.values)
Base.size(actor::BufferActor)        = size(actor.values)
Base.length(actor::BufferActor)      = length(actor.values)
Base.getindex(actor::BufferActor, I)    = Base.getindex(actor.values, I)
Base.getindex(actor::BufferActor, I...) = Base.getindex(actor.values, I...)

Base.firstindex(actor::BufferActor) = firstindex(actor.values)
Base.lastindex(actor::BufferActor)  = lastindex(actor.values)


text:
version = v"1.4.1"
api_version = v"1.2.13"

# Collection of sources required to build this package
#
# The level zero repository contains both the API headers and a loader,
# which are versioned independently.
sources = [
    GitSource("https://github.com/oneapi-src/level-zero.git",
              "551dd5810a3cea7a7e26ac4441da31878e804b53"),
]

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
# Boys.jl contains different implementations and approximations to the Boys function
using Interpolations

function boys_array_asymp(mmax,x)
    boys_array = zeros(Float64,mmax)
    denom = 0.5/x # Only used for large x, so don't check for small x
    boys_array[1] = sqrt(0.5*pi*denom)
    for m in 2:mmax
        boys_array[m] = boys_array[m-1]*(2m-1)*denom
    end
    return boys_array
end

function boys_array_gamma(mmax,x,SMALL=1e-18)
    x = max(x,SMALL) # needs underflow protection because of inverse
    boys_array = zeros(Float64,mmax)
    oox = 1/x
    denom = sqrt(oox)
    boys_array[1] = 0.5*denom*gammainc(0.5,x)
    for m in 2:mmax
        denom *= oox
        # Could speed this up more by expressing gamma(m) in terms of gamma(m±1)
        boys_array[m] = 0.5*denom*gammainc(m-0.5,x)
    end
    return boys_array
end

boys_array_fmref(mmax,x) = [fm_ref(m-1,x) for m in 1:mmax]
boys_array_Fgamma(mmax,x) = [Fgamma(m-1,x) for m in 1:mmax]

"Boys Fgamma function, using the lower incomplete gamma function."
@inline function Fgamma(m,T,SMALL=1e-18,Tcrit=20.0)
    # Note, most programs use a much larger value for Tcrit (117)
    mhalf = m+0.5
    T = max(T,SMALL) # Evidently needs underflow protection
    if T>Tcrit
        retval = sqrt(pi/2)*factorial2(2m-1)/(2T)^mhalf
    else
        #retval = 0.5*T^-mhalf*gammainc_nr(mhalf,T)
        retval = fmspline(m,T)
    end
    return retval
end

function Fm(m,T,SMALL=1e-18)
    mhalf = m+0.5
    T = max(T,SMALL)
    return 0.5*T^-mhalf*gammainc_nr(mhalf,T)
end

"gammainc returns the lower incomplete gamma function"
@inline gammainc(a,x) = gamma(a)*gamma_inc(a,x)[1]

# This is the ref function from libint
function fm_ref(m,T,eps = 1e-10)
    denom = (m + 0.5)
    term = exp(-T) /2denom
    old_term = 0.0
    sum = term
    while term > sum*eps || old_term < term
        denom += 1
        old_term = term
        term = old_term * T / denom
        sum += term
    end
    return sum
end

# These were taken from Numerical Recipes
@inline function gammainc_nr(a::Float64,x::Float64)
    if abs(a-1) < 1e-3
        println("Warning: gammainc_series is known to have problems for a ~ 1")
    end
    if x < (a+1.0)
        #Use the series representation
        gam,gln = gser(a,x)
    else
        #Use continued fractions
        gamc,gln = gcf(a,x)
        gam = 1-gamc
    end
    return exp(gln)*gam
end

@inline function gser(a,x,ITMAX=100,EPS=3e-9)
    # Series representation of Gamma. NumRec sect 6.1.
    gln=loggamma(a)
    if x == 0
        return 0,gln
    end
    ap = a
    delt = s = 1/a
    for i in 1:ITMAX
        ap += 1
        delt *= (x/ap)
        s += delt
        if abs(delt) < abs(s)*EPS
            break
        end
    end
    return s*exp(-x+a*log(x)-gln),gln
end

@inline function gcf(a::Float64,x::Float64,ITMAX::Int64=200,EPS::Float64=3e-9,FPMIN::Float64=1e-30)
    #Continued fraction representation of Gamma. NumRec sect 6.1"
    gln=loggamma(a)
    b=x+1.0-a
    c=1.0/FPMIN
    d=1.0/b
    h=d
    for i in 1:ITMAX
        an=-i*(i-a)
        b=b+2.
        d=an*d+b
        if abs(d) < FPMIN
            d=FPMIN
        end
        c=b+an/c
        if abs(c) < FPMIN
            c=FPMIN
        end
        d=1.0/d
        delt=d*c
        h=h*delt
        if abs(delt-1.) < EPS
            break
        end
    end
    gammcf = exp(-x+a*log(x)-gln)*h
    return gammcf,gln
end

function make_interpolator(mmax=10,Tmax=20.0)
    Tgrid = 0:0.005:Tmax
    mgrid = 0:mmax
    fmvalues = [Fm(m,T) for m in mgrid, T in Tgrid]
    itp = interpolate(fmvalues, BSpline(Cubic(Line(OnGrid()))))
    sitp = scale(itp,mgrid,Tgrid)
    return sitp
end
const fmspline = make_interpolator()


text:
abstract type BasisFunction end
abstract type Linear  <: BasisFunction end
abstract type Quadratic <: BasisFunction end

struct NodeLocation
    l::UInt8
    j::UInt16
    dj::UInt16
end
NodeLocation(x::Float64) = NodeLocation(position(x)...)

mutable struct AdaptiveGrid
    nodeinfo
    active::BitArray{1}
    overlap::Vector{Vector{Bool}}
end

mutable struct NGrid{D, B}
    L::Vector{Int}
    bounds::Array{Float64, 2}
    grid::Array{Float64, 2}
    covers::Array{UInt16, 2}
    covers_dM::Array{UInt16, 2}
    covers_loc::Vector{UInt32}

    adapt::AdaptiveGrid

    IDs::Vector{Vector{Int}}
    Bs::Vector{Vector{Float64}}
end

# Returns the tensor grid of a vector of univariate grids
function Base.kron(X::Vector{Vector{Float64}})
    D = length(X)
    lxs = Int[length(x) for x in X]
    G = zeros(prod(lxs), D)
    s = 1
    for d = 1:D
        snext = s * lxs[d]
        for j = 1:prod(lxs)
            G[j, d] = X[d][div(rem(j - 1, snext), s) + 1]
        end
        s = snext
    end
    return G
end

Base.kron(x::Vector{Float64}) = x

"""
    TensorGrid(L)

Returns a tensor grid of univariate grids of levels given by L = [l1,l2,...]
"""
function TensorGrid(L::Vector{Int})
    G = ndgrid(Vector{Float64}[cc_g(i) for i in L]...)
    G = hcat([vec(g) for g in G]...)
end

"""
    SmolyakSize(L[, mL])

Returns the number of nodes in levels mL of a Smolyak grid of maximum
of maximum level L = [l1,l2...]
"""
function SmolyakSize(L::Vector{Int}, mL::UnitRange{Int} = 0:maximum(L))
    D = length(L)
    m = Int[dM(l) for l = 1:maximum(L) + D]
    S = 0
    for l = mL
        for covering in comb(D, D + l)
            if all(covering .≤ L + 1)
                s = m[covering[1]]
                for i = 2:length(covering)
                    s *= m[covering[i]]
                end
                S += s
            end
        end
    end
    return S
end

"""
    SmolyakGrid(L[, mL])

 juli
"""
function SmolyakGrid(L::Vector{Int}, mL::UnitRange{Int} = 0:maximum(L))
    D = length(L)
    dg = Vector{Float64}[cc_dg(l) for l = 1:maximum(L) + D]
    G = Array{Array{Float64}}(0)
    index = Array{Array{Int}}(0)
    for l = mL
        for covering in comb(D, D + l)
            if all(covering .≤ L + 1)
                push!(G, kron(dg[covering]))
                push!(index, repmat(covering', size(G[end], 1)))
            end
        end
    end
    # G = vcat(G...)::Array{Float64,2}
    # index = vcat(index...)::Array{Int,2}
    # return G,index
    return (vcat(G...), vcat(index...))::Tuple{Array{Float64, 2}, Array{Int, 2}}
end

"""
    level(x)

Computes the minimum level of a point
"""
function level(x::Float64)
    l = 0
    if x == 0.5
        l =  1
    elseif x == 0.0 || x == 1.0
        l = 2
    else
        for l = 3:12
            mod(x, 1 / 2^(l - 1)) == 0.0 && break
        end
    end
    return l
end
level(G::NGrid) = level(G.grid)
level(X::Array{Float64}) = vec(sum(map(level, X), 2) - size(X, 2))

"""
    NGrid(L,bounds,B)

# Construction
Construct a Smolyak grid.
The vector of integers L determines the maximum level of the grid in each
dimension. The bounds can be optionally changed from the default of [0,1]^D and
basis function B can be either Linear or Quadratic.

# Interpolation
Grid objects are callable taking two arguments. The first is a vector containing
function values at the grid nodes. The second array contains rows of points at
which the interpolant is to be evaluated.
"""
function NGrid(L::Vector{Int}, bounds::Array{Float64} = [0, 1] .* ones(1, length(L));B::Type{BT} = Linear) where {BT <: BasisFunction}
    grid, ind = SmolyakGrid(L)
	   covers = UInt16.(unique(ind, 1))
    covers_loc = zeros(UInt32, size(covers, 1))
    hind = vec(mapslices(hash, ind, 2))
    hcovers = vec(mapslices(hash, covers, 2))
    for i = 1:size(covers, 1)
        covers_loc[i] = findfirst(hind, hcovers[i])
    end

    G = NGrid{length(L), B}(L,
                        bounds,
                        grid,
                        covers,
                        dM.(covers),
                        covers_loc,
                        AdaptiveGrid([], zeros(Bool, size(grid, 1)), []),
                        Vector{Int}[Int[] for i = 1:size(grid, 1)],
                        Vector{Float64}[Float64[] for i = 1:size(grid, 1)]
                        )
    G.adapt.active = (level(G) .== maximum(level(G)))
    buildW(G, hind, hcovers)
    return G
end

Base.show(io::IO, G::NGrid) = println(io, typeof(G), ": $(size(G.grid,1))pts")
Base.length(G::NGrid) = size(G.grid, 1)
Base.size(G::NGrid) = size(G.grid)
Base.values(G::NGrid) = nUtoX(G.grid, G.bounds)

Base.getindex(G::NGrid, args...) = getindex(G.grid, args...)


text:
###############################################################################
#
# SoilPlantAirContinuum system for simple mono species system
#
###############################################################################
"""
    mutable struct SPACSimple{FT}

Struct that simplest mono species SoilPlantAirContinuum system, with 1 root,
    stem, and leaf.

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct SPACSimple{FT<:AbstractFloat}
    # Hydraulic system
    "Hydraulic systems"
    hs::TreeSimple{FT} = TreeSimple{FT}()
    "Critical flow rate"
    ec::FT = FT(50);

    # Photosynthesis parameters
    "Photosynthesis system"
    ps::Leaf{FT} = Leaf{FT}()
    "Ratio between Vcmax25 and Jmax25"
    vtoj::FT = FT(1.67)

    # Surrounding AirLayer
    "Environmental conditions"
    envir::AirLayer{FT} = AirLayer{FT}()

    # local container for returned results
    "Container for gas exchange for a layer"
    container1L::SPACContainer1L{FT} = SPACContainer1L{FT}()
    "Container for gas exchange of sunlit and shaded layers"
    container2L::SPACContainer2L{FT} = SPACContainer2L{FT}()
    "Container for default hydraulic conductance"
    containerKS::Array{FT,1} = FT[hs.root.k_max, hs.stem.k_max, hs.leaf.k_sla]
    "Container for optimizer"
    containerOP::FT = FT(0)
    "Container for optimal sunlit layer flow rate"
    opt_f_sl   ::FT = FT(0)
    "Container for optimal shaded layer flow rate"
    opt_f_sh   ::FT = FT(0)
    "Container for optimal leaf area per basal area"
    opt_laba   ::FT = FT(2000)
    "Container for optimal Vcmax25"
    opt_vmax   ::FT = FT(5)

    # leaf related
    "Leaf area index"
    lai  ::FT = FT(3)
    "Leaf area per basal area"
    laba ::FT = 1500
    "Maximal stomatal conductance limit at 25 °C"
    g_max::FT = 0.8
    "Ground area per basal area"
    gaba ::FT = 500
    "Leaf width"
    width::FT = 0.05

    # soil related
    "Maximal soil water content"
    mswc  ::FT = (hs.root.sh).Θs
    "Current soil water content"
    swc   ::FT = (hs.root.sh).Θs
    "Soil matrical water potential"
    p_soil::FT = 0.0
    "Soil depth, 2X mean root depth"
    h_soil::FT = 2

    # leaf invest related
    "Leaf construction cost per leaf area"
    c_cons::FT = 2
    "Leaf nutrient cost per Vcmax25 per leaf area"
    c_vmax::FT = 0.04

    # geography related
    "Latitude `[°]`"
    latitude ::FT = 30
    "Longitude `[°]`"
    longitude::FT = 116
    "Elevation `[m]`"
    elevation::FT = 0
end

max_stars_repo_name: climate-machine/Land


text:
"""
    `nanoribbonN(p, (nx, ny); mono = false)`
builds a normal infinite graphene nanoribbon with Bravais vector: `br =(nx, ny)` and parameters:
`p::Params()`. `mono` specifies whether we have SLG (`mono = true`) or
BLG in Bernal stacking (`mono = false`).
    see: `Params()`, `modelN()`
"""
function nanoribbonN(p, (nx, ny); mono = false)
    (; Ln, Ls, U, a0) = p
    (; model0, field!) = modelN(p)
    lat = mono ? latSLG(p) : latBLG_unbounded(p)
    perpaxis = Quantica.normalize(bravais(lat) * SA[1 1; -1 1] * SA[-ny, nx])

    regionS(r) = abs(dot(perpaxis, r)) > ifelse(iszero(Ls), (round(0.5*Ln/a0)-0.9)*a0, Ln/2)
    smooth(r, β = 5) = iszero(Ls) ? 1.0 : tanh(β*(abs(dot(perpaxis, r)) - Ln/2) / Ls)

    # U! = @onsite!((o, r) -> o + U * sign(r[3]) * smooth(r) * σ0; region = regionS)

    h = lat |> hamiltonian(model0; orbitals = Val(2)) |>
               unitcell((1, -1), (1, 1)) |>
               unitcell((nx, ny), region = r -> abs(dot(perpaxis, r)) < Ln/2 + Ls,
                modifiers = (field!,))
    return h
end

"""
    `nanoribbonS(p, (nx, ny); mono = false, dxS = 10.)`
builds an infinite graphene nanoribbon with Bravais vector: `br =(nx, ny)`, parameters:
`p::Params()` with induced superconductivity along the edges with a characteristic
penetration given by: `dxS` in units of a0.
`mono` keyword specifies whether we have SLG (`mono = true`) or BLG in Bernal stacking
(`mono = false`).
    see: `Params()`, `modelS()`
"""
function nanoribbonS(p, (nx, ny); mono = false, dxS = 0.5)
    (; Ln, Ls, Δ, a0, τ, d, Ws) = p
    (; model0, field!) = modelS(p)
    lat = mono ? latSLG(p) : latBLG_unbounded(p)
    perpaxis = Quantica.normalize(bravais(lat) * SA[1 1; -1 1] * SA[-ny, nx])

    regionS(r) = abs(dot(perpaxis, r)) > Ln/2-Ls
    regionS(r, dr) = regionS(r)

    smooth(r) = iszero(d) ? ifelse(abs(dot(perpaxis, r)) < Ln/2 - dxS * a0, 0.0, 1.0) :
                            1-(1+tanh((Ln/2-abs(dot(perpaxis, r)))/d))/2

    SCo! = @onsite!((o, r) -> o + smooth(r) * Δ * σyτy)
    SCh! = @hopping!((t, r, dr) -> t * ifelse(iszero(d),  1.0, (1.0 - smooth(r)));
        sublats = (:At, :Bt ,:Ab, :Bb, :A, :B) .=> (:At, :Bt ,:Ab, :Bb, :A, :B))
    SCτ! = @hopping!(t -> t * τ; region = (r, dr) -> xor(regionS(r - dr/2),
        regionS(r + dr/2)))

    h = lat |> hamiltonian(model0; orbitals = Val(4)) |>
               unitcell((1, -1), (1, 1)) |>
               unitcell((nx, ny), region = r -> abs(dot(perpaxis, r)) < Ln/2 + Ls,
               modifiers = (field!, SCo!, SCh!, SCτ!))
    return h
end

nanoribbonNA(p = Params(); kw...) = nanoribbonN(p, (0, 1); kw...)
nanoribbonNZ(p = Params(); kw...) = nanoribbonN(p, (1, 0); kw...)
nanoribbonSA(p = Params(); kw...) = nanoribbonS(p, (0, 1); kw...)
nanoribbonSZ(p = Params(); kw...) = nanoribbonS(p, (1, 0); kw...)

##################################
# Deprecated: rotated nanoribbons
##################################

# function nanoribbonN(p, (nx, ny))
#     (; Ln, Ls) = p
#     (; model0, field!) = modelN(p)
#     lat_top, lat_bot = latBLG(p)
#     perpaxis = Quantica.normalize(bravais(lat) * SA[1 1; -1 1] * SA[-ny, nx])

#     h = lat |> hamiltonian(model0; orbitals = Val(2)) |>
#                unitcell((1, -1), (1, 1)) |>
#                unitcell((nx, ny), region = r -> abs(dot(perpaxis, r)) < Ln/2 + Ls,
#                 modifiers = (field!,))
#     return h
# end

# function nanoribbonS(p, θ)
#     (; Ln, Ls, Δ, a0, τ, d, Ws) = p
#     (; model0, field!, modelinter) = modelS(p)
#     lat_top, lat_bot = latBLG(p, θ)

#     # regionS(r) = abs(dot(perpaxis, r)) > ifelse(iszero(Ls),
            # (round(0.5*Ln/a0)-WS)*a0, Ln/2) #0.9
#     # regionS(r, dr) = regionS(r)
#     regionS(r) = abs(r[1]) >= Ln/2 - Ls || abs(r[2]) >= W/2 - Ws
#     regionS(r, dr) = regionS(r)

#     smooth(r) = iszero(d) ? ifelse(abs(dot(perpaxis, r)) < Ln/2 - dxS * a0, 0.0, 1.0) :
#                             1-(1+tanh((Ln/2-abs(dot(perpaxis, r)))/d))/2

#     SCo! = @onsite!((o, r) -> o + smooth(r) * Δ * σyτy)
#     SCh! = @hopping!((t, r, dr) -> t * ifelse(iszero(d),  1.0, (1.0 - smooth(r)));
#         sublats = (:At, :Bt ,:Ab, :Bb, :A, :B) .=> (:At, :Bt ,:Ab, :Bb, :A, :B))
#     SCτ! = @hopping!(t -> t * τ; region = (r, dr) -> xor(regionS(r - dr/2),
#         regionS(r + dr/2)))

#     h_top = lat_top |> hamiltonian(model0; orbitals = Val(4))
#     h_bot = lat_bot |> hamiltonian(model0; orbitals = Val(4))
#     h = Quantica.combine(h_top, h_bot, coupling = modelinter) |>
#         unitcell((1, -1), (1, 1), modifiers = (field!, SCo!, SCh!, SCτ!))

#     # h = lat |> hamiltonian(model0; orbitals = Val(4)) |>
#     #            unitcell((1, -1), (1, 1)) |>
#     #            unitcell((nx, ny), region = r -> abs(dot(perpaxis, r)) < Ln/2 + Ls,
#     #            modifiers = (field!, SCo!, SCh!, SCτ!))
#     return h
# end

# nanoribbonNA(p = Params(); kw...) = nanoribbonN(p, 0; kw...)
# nanoribbonNZ(p = Params(); kw...) = nanoribbonN(p, π/2; kw...)
# nanoribbonSA(p = Params(); kw...) = nanoribbonS(p, 0; kw...)
# nanoribbonSZ(p = Params(); kw...) = nanoribbonS(p, π/2; kw...)


text:
module MLPreprocessing

using Statistics
using StatsBase
using LearnBase
using DataFrames

using LearnBase: ObsDimension, default_obsdim

export

    ObsDim,
    expand_poly,

    center!,
    standardize!,
    fixedrange!,

    StandardScaler,
    FixedRangeScaler,
    fit,
    fit_transform,
    fit_transform!,
    transform,
    transform!

include("scaleselection.jl")
include("basis_expansion.jl")
include("center.jl")
include("standardize.jl")
include("fixedrange.jl")

end # module


text:
include("./Types.jl")
import Polyhedra, CDDLib
export primal
"""
    primal(Γ, P1, P2, D2, Cone, Dualcone ; ϵ=0.05)

implements Algorithm1 in
Löhne, A., Rudloff, B. & Ulus, F. Primal and dual approximation algorithms for convex vector optimization problems. J Glob Optim 60, 713–736 (2014). [Download](https://doi.org/10.1007/s10898-013-0136-0)
"""
function primal(Γ, P1, P2, D2, Cone, Dualcone; ϵ = 0.05)
    # any linearly independent q vectors in cone C in R^q
    # c1, c2, ..., cq-1, c where c in intC
    C = Cone
    # q geberators of vectors of the dual C+
    (_, q) = size(Cone)
    T = C
    Z = Dualcone
    # lines 1 and 3
    Xbar = [first(P1(Z[:, i])) for i = 1:q]
    # line 2
    Hs = [Polyhedra.HalfSpace(-Z[:, i], Z[:, i]' * Xbar[i]) for i = 1:q]
    P_hr = Polyhedra.hrep(Hs)
    P_poly = Polyhedra.polyhedron(P_hr)
    # line 3
    k = 1
    Tinv = inv(T)
    Tbar = [Float64.(T' * Z[:, i]) for i = 1:q]
    # line 4
    # D = D*(Tbar) = [Tbar₁, Tbar₂, ⋯, Tbar_q,   ]
    D = hcat([Tbar[i][1:q-1][1] for i in 1:q], [Z[:, i]' * Γ(Xbar[i]) for i = 1:q])

    # bags to collect the vertices
    Vknown = []
    # flag to break the loop
    go = true
    while go
        # line 7
        P_vr = Polyhedra.vrep(P_poly)
        # extract number of vecrtices from P_vr
        # that is |Pⱽ|
        num_v = Polyhedra.npoints(P_vr)
        # extract vecrtices from P_vr as row
        v = ((Polyhedra.MixedMatVRep(P_vr)).V)'
        vv = [v[:, i] for i = 1:num_v]
        ## by default stop the loop
        go = false

        # for each vertex found
        for v0 in vv
            # check if this vertex has been visited before?
            if ~(Vertex2D(v0...) ∈ Vertex2D.(Vknown))
                # if the vertex has not been visited
                # v0_as_vertic = Vertex2D(v0...)
                ## store this vertex as visited
                # Vknown[k] = v0_as_vertic
                ## collect it in the bag
                push!(Vknown, v0)
                # increase the index k
                k = k + 1
                # line 10
                (xv0, zv0) = P2(v0)
                (wv0, uv0) = D2(v0)
                # update Xbar and Tbar line 11
                push!(Xbar, xv0)
                tv0 = T'wv0
                push!(Tbar, tv0)
                Dstar_q = (wv0'*v0)[1] + zv0
                D = hcat(D, [tv0[1]; Dstar_q])
                # since zv0>ϵ we compute the huperplane in line 13
                if (zv0 > ϵ)
                    # H0Mat = [tv0[1:end-1][1] 1] * Tinv
                    H0Mat = (vcat(tv0[1:end-1], 1))' * Tinv
                    # if H = {y∈ Rq | aᵀ*y≤ b}, we construct it as
                    # Polyhedra.HalfSpace(a, b)
                    # if Hⱼ = {y∈ Rq | (zʲ)ᵗy≥ (zʲ)ᵗΓ(xⱼ)}
                    ## This halfspace is constructed as follows
                    # The HalfSpace of Polyhedra package demands that
                    #    Hⱼ = {y∈ Rq | ϕ(y,D*(Tᵀzʲ))≥0} = {y∈ Rq | ϕ(y,D*(Tᵀzʲ))≥0}
                    #    Hⱼ = {y∈ Rq | [(Tᵀzʲ)₁ (Tᵀzʲ)₂ ⋯ (Tᵀzʲ)q 1]*Tinv*y -D*(Tᵀzʲ)≥0}
                    #    Hⱼ = {y∈ Rq | [(Tᵀzʲ)₁ (Tᵀzʲ)₂ ⋯ (Tᵀzʲ)q 1]*Tinv*y≥D*(Tᵀzʲ)}
                    #    Hⱼ = {y∈ Rq | H0Mat*y≥Dstar_q}
                    # we write as
                    #    Hⱼ = {y∈ Rⁿ | - (zⱼ)ᵗy≤ - (zⱼ)ᵗΓ(xⱼ)}
                    # ∴  Hⱼ = {y∈ Rq | -H0Mat*y≤ -Dstar_q}
                    # i.e aᵀ = -H0Mat ⇒ a = -H0Mat'
                    #     b = -Dstar_q
                    H0 = Polyhedra.HalfSpace(-H0Mat', -Dstar_q)
                    P_hr = P_poly.hrep ∩ H0
                    ## we compute P_{k+1} in line 18 and break the loop
                    P_poly = Polyhedra.polyhedron(P_hr)
                    go = true
                end
            end
        end

    end
    return Xbar, Tbar, Vknown, P_poly
end


text:
using CompEcon
using Dierckx
using Statistics

# Create functions to test against
f1(x) = sin.(x)
fp1(x) = cos.(x)
f2(x) = exp.(x)
fp2(x) = exp.(x)
f3(x) = 1.5 .+ log.(x)
fp3(x) = 1.0 ./ x
fs = [f1, f2, f3]
fps = [fp1, fp2, fp3]
nf = length(fs)

x = collect(range(.5; stop=2, length=35))
finex = collect(range(1e-2; stop=2*pi - 1e-2, length=150))
ys = map(f->f(x), fs)
finey = map(f->f(finex), fs)

function compare_D_CE_levels(f::Function, x::Array{Float64, 1}, order::Int)

    # Create Dierckx Spline (and data)
    n = length(x)
    y = f(x)
    dspl = Spline1D(x, y; s=0., k=order)

    # Create CompEcon Spline (and data)
    cebasis = Basis(SplineParams(x, 0, order))
    xx = nodes(cebasis)[1]
    yy = f(xx)
    c, bs = funfitxy(cebasis, xx, yy)

    # Evaluate Splines on finer grid
    xfine = collect(range(x[1]; stop=x[end], length=2*n + 1))
    yfine = f(xfine)

    deval = Dierckx.evaluate(dspl, xfine)
    derr = abs.(yfine - deval)

    ceeval = CompEcon.funeval(c, cebasis, xfine)
    ceerr = abs.(yfine - ceeval)

    # # Evaluate Derivative
    # dderiv1 = Dierckx.derivative(dspl1, finex)
    # cederiv1 = CompEcon.evaluate(ceinterp1, finex, order=1)

    return derr, ceerr
end

function compare_D_CE_deriv(f::Function, fp::Function, x::Array{Float64, 1}, order::Int)

    # Create Dierckx Spline (and data)
    n = length(x)
    y = f(x)
    dspl = Spline1D(x, y; s=0., k=order)

    # Create CompEcon Spline (and data)
    cebasis = Basis(SplineParams(x, 0, order))
    xx = nodes(cebasis)[1]
    yy = f(xx)
    c, bs = funfitxy(cebasis, xx, yy)

    # Evaluate Splines on finer grid
    xfine = collect(range(x[1]; stop=x[end], length=2*n + 1))
    ypfine = fp(xfine)

    devalderiv = Dierckx.derivative(dspl, xfine)
    derr = abs.(ypfine - devalderiv)

    ceevalderiv = funeval(c, cebasis, xfine, ones(Int64, 1,1))
    ceerr = abs.(ypfine .- ceevalderiv)

    return derr, ceerr
end

for i=1:nf

    curr_f, curr_fp = fs[i], fps[i]

    derr1_levels, ceerr1_levels = compare_D_CE_levels(curr_f, x, 1)
    derr2_levels, ceerr2_levels = compare_D_CE_levels(curr_f, x, 2)
    derr3_levels, ceerr3_levels = compare_D_CE_levels(curr_f, x, 3)

    # derr1_deriv, ceerr1_deriv = compare_D_CE_deriv(curr_f, curr_fp, x, 1)
    derr2_deriv, ceerr2_deriv = compare_D_CE_deriv(curr_f, curr_fp, x, 2)
    derr3_deriv, ceerr3_deriv = compare_D_CE_deriv(curr_f, curr_fp, x, 3)

    println("\nLinear Approximation to Function\n")
    println("Dierckx Level Errors \n\t Max: $(maximum(derr1_levels)) \n\t Min: $(minimum(derr1_levels)) \n\t Mean: $(mean(derr1_levels))")
    # println("Dierckx Deriv Errors \n\t Max: $(maximum(derr1_deriv)) \n\t Min: $(minimum(derr1_deriv)) \n\t Mean: $(mean(derr1_deriv))")
    println("CompEcon Level Errors \n\t Max: $(maximum(ceerr1_levels)) \n\t Min: $(minimum(ceerr1_levels)) \n\t Mean: $(mean(ceerr1_levels))")
    # println("CompEcon Deriv Errors \n\t Max: $(maximum(ceerr1_deriv)) \n\t Min: $(minimum(ceerr1_deriv)) \n\t Mean: $(mean(ceerr1_deriv))")

    println("\n\nQuadratic Approximation to Function\n")
    println("Dierckx Level Errors \n\t Max: $(maximum(derr2_levels)) \n\t Min: $(minimum(derr2_levels)) \n\t Mean: $(mean(derr2_levels))")
    println("Dierckx Deriv Errors \n\t Max: $(maximum(derr2_deriv)) \n\t Min: $(minimum(derr2_deriv)) \n\t Mean: $(mean(derr2_deriv))")
    println("CompEcon Level Errors \n\t Max: $(maximum(ceerr2_levels)) \n\t Min: $(minimum(ceerr2_levels)) \n\t Mean: $(mean(ceerr2_levels))")
    println("CompEcon Deriv Errors \n\t Max: $(maximum(ceerr2_deriv)) \n\t Min: $(minimum(ceerr2_deriv)) \n\t Mean: $(mean(ceerr2_deriv))")

    println("\n\nCubic Approximation to Function\n")
    println("Dierckx Level Errors \n\t Max: $(maximum(derr3_levels)) \n\t Min: $(minimum(derr3_levels)) \n\t Mean: $(mean(derr3_levels))")
    println("Dierckx Deriv Errors \n\t Max: $(maximum(derr3_deriv)) \n\t Min: $(minimum(derr3_deriv)) \n\t Mean: $(mean(derr3_deriv))")
    println("CompEcon Level Errors \n\t Max: $(maximum(ceerr3_levels)) \n\t Min: $(minimum(ceerr3_levels)) \n\t Mean: $(mean(ceerr3_levels))")
    println("CompEcon Deriv Errors \n\t Max: $(maximum(ceerr3_deriv)) \n\t Min: $(minimum(ceerr3_deriv)) \n\t Mean: $(mean(ceerr3_deriv))")

end


text:
mutable struct Swarm{T<:AbstractFloat}

    # Vector of particles
    particles::Vector{Particle{T}}

    # Preallocated vector of UInt16 for neighborhood selection
    nVec::Vector{Int}

    # Global best objective function value
    b::T

    # Location of global best objective function value
    d::Vector{T}

    n::Int      # Neighborhood size
    w::T        # Inertia
    c::Int      # Adaptive inertia counter

    y₁::T   # Self adjustment weight
    y₂::T   # Social adjustment weight

end

function Swarm{T}(::UndefInitializer) where {T}
        return Swarm{T}(Vector{Particle{T}}(undef, 0),
            Vector{Int}(undef, 0), T(0.0), Vector{T}(undef, 0),
            Int(0), T(0.0), Int(0), T(0.0), T(0.0))
end

function Swarm{T}(nDims::Integer, nParticles::Integer) where {T}
    if nDims < 0
        throw(ArgumentError("nDims must be greater than 0."))
    end
    if nParticles < 0
        throw(ArgumentError("nParticles must be greater than 0."))
    end
    return Swarm{T}(Vector{Particle{T}}([Particle{T}(nDims) for i in 1:nParticles]),
        Vector{Int}(1:nParticles), T(0.0), Vector{T}(undef, nDims),
        Int(0), T(0.0), Int(0), T(0.0), T(0.0))
end

# ===== Methods

Base.length(s::Swarm) = length(s.particles)

function Base.getindex(s::Swarm, i::Int)
    1 <= i <= length(s.particles) || throw(BoundsError(s, i))
    return s.particles[i]
end

function Base.setindex!(s::Swarm, v, i::Int)
    1 <= i <= length(s.particles) || throw(BoundsError(s, i))
    s.particles[i] = v
end

function feval!(s::Swarm, f, opts::Options; init = false)
    # Evaluate objective functions
    @inbounds begin
        if opts.useParallel
            @qthreads for i in 1:length(s)
                s[i].fx = f(s[i].x)
            end
        else
            @simd for i in 1:length(s)
                s[i].fx = f(s[i].x)
            end
        end
    end

    # Check objective function values if desired
    if opts.funValCheck
        fValCheck(s)
    end

    # Update each particles best objective function value and its location
    if !init
        @inbounds @simd for i in 1:length(s)
            if s[i].fx < s[i].fp
                s[i].p .= s[i].x
                s[i].fp = s[i].fx
            end
        end
    else
        @inbounds @simd for i in 1:length(s)
            s[i].fp = s[i].fx
        end
    end
end

function setGlobalBest!(s::Swarm)
    updated = false
    @inbounds begin
        @simd for i in 1:length(s)
            if s[i].fp < s.b
                s.b = s[i].fp
                s.d .= s[i].p
                updated = true
            end
        end
    end
    return updated
end

function fValCheck(s::Swarm)
    @inbounds begin
        @simd for i in 1:length(s)
            if isinf(s[i].fx) || isnan(s[i].fx)
                throw(ErrorException("Objective function returned Inf or NaN!"))
            end
        end
    end
end

function updateVelocities!(s::Swarm)
    n = length(s.d)
    m = length(s)
    @inbounds begin
        @simd for i in 1:m

            # Shuffle vector containing integers 1:n
            # first m != i will be neighborhood
            shuffle!(s.nVec)

            # Determine fbest(S)
            fbest = Inf
            best = 0
            incr = 0
            for j in 1:s.n
                s.nVec[j] == i ? incr = 1 : ()
                k = j + incr
                if s[k].fp < fbest
                    fbest = s[k].fp
                    best = k
                end
            end

            # Update i's velocity
            @simd for j in 1:n
                s[i].v[j] = s.w*s[i].v[j] + s.y₁*rand()*(s[i].p[j] - s[i].x[j]) +
                                s.y₂*rand()*(s[best].p[j] - s[i].x[j])
            end
        end
    end
end

function step!(s::Swarm)
    @inbounds @simd for i in 1:length(s)
        s[i].x .= s[i].x .+ s[i].v
    end
end

function enforceBounds!(s::Swarm, LB, UB)
    @inbounds begin
        @simd for i in 1:length(s)
            @simd for j in 1:length(s.d)
                if s[i].x[j] > UB[j]
                    s[i].x[j] = UB[j]
                    s[i].v[j] = 0.0
                elseif s[i].x[j] < LB[j]
                    s[i].x[j] = LB[j]
                    s[i].v[j] = 0.0
                end
            end
        end
    end
end

function printStatus(s::Swarm, time::AbstractFloat, iter::Int, stallCount::Int)
    fspec1 = FormatExpr("Time Elapsed: {1:f} sec, Iteration Number: {2:d}, Function Evaluations: {3:d}")
    fspec2 = FormatExpr("Stall Iterations: {1:d}, Global Best: {2:e}")
    printfmtln(fspec1, time, iter, (iter + 1)*length(s))
    printfmtln(fspec2, stallCount, s.b)
    println(" ")
end

# Initializes position and velocities of particles by sampling from a
# uniform distribution
function uniformInitialization!(swarm::Swarm, prob::Problem, opts::Options)

    # Get N: Number of diamensions and M: Swarm Size
    M = length(swarm)
    N = length(swarm[1])

    # Get Boundary Constraints
    LB  = prob.LB
    UB  = prob.UB

    # Check if initial bounds on positions have been set
    useInitBnds = false
    if length(opts.iLB) == N && length(opts.iUB) == N
        useInitBnds = true
        iLB = opts.iLB
        iUB = opts.iUB
    end

    # Initialize particle positions and velocities
    @inbounds begin
        @simd for d in 1:N
            # Get local bounds for d-axis
            lLB = useInitBnds ? (LB[d] < iLB[d] ? iLB[d] : LB[d]) : LB[d]
            lUB = useInitBnds ? (UB[d] > iUB[d] ? iUB[d] : UB[d]) : UB[d]
            @simd for p in 1:M
                # Position information
                swarm[p].x[d] = lLB + (lUB - lLB)*rand()
                swarm[p].p[d] = swarm[p].x[d]

                # Velocity
                r = useInitBnds ? min(lUB-lLB,UB[d]-LB[d]) : lUB - lLB
                swarm[p].v[d] = -r + 2*r*rand()
            end
        end
    end

    return nothing
end

# Initializes position and velocities of particles using logistic map
function logisticsMapInitialization!(swarm::Swarm, prob::Problem, opts::Options)

    # Get N: Number of diamensions and M: Swarm Size
    M = length(swarm)
    N = length(swarm[1])

    # Get Boundary Constraints
    LB  = prob.LB
    UB  = prob.UB

    # Check if initial bounds on positions have been set
    useInitBnds = false
    if length(opts.iLB) == N && length(opts.iUB) == N
        useInitBnds = true
        iLB = opts.iLB
        iUB = opts.iUB
    end

    # Logistics map initialization
    fixedPointTol = 1e-14
    maxPert = 1e-12
    lMapIters = 3000
    @inbounds begin
        @simd for j in 1:M
            @simd for k in 1:N
                swarm[j].x[k] = 0.4567 + 2*(rand() - 0.5)*maxPert
            end
        end
        @simd for i in 1:lMapIters
            @simd for j in 1:M
                @simd for k in 1:N
                    val = swarm[j].x[k]
                    if val < fixedPointTol ||
                       abs(val - 0.25) < fixedPointTol ||
                       abs(val - 0.50) < fixedPointTol ||
                       abs(val - 0.75) < fixedPointTol ||
                       abs(val - 1.00) < fixedPointTol

                       swarm[j].x[k] += maxPert*rand()
                    end
                    swarm[j].x[k] = lMap(swarm[j].x[k])
                    if isinf(swarm[j].x[k])
                        throw(ErrorException("Inf or NaN"))
                    end
                end
            end
        end
    end

    # Scale particle positions and initialize velocities
    @inbounds begin
        @simd for d in 1:N
            # Get local bounds for d-axis
            lLB = useInitBnds ? (LB[d] < iLB[d] ? iLB[d] : LB[d]) : LB[d]
            lUB = useInitBnds ? (UB[d] > iUB[d] ? iUB[d] : UB[d]) : UB[d]
            @simd for p in 1:M
                # Position information
                swarm[p].x[d] = lLB + (lUB - lLB)*swarm[p].x[d]
                swarm[p].p[d] = swarm[p].x[d]

                # Velocity
                r = useInitBnds ? min(lUB-lLB,UB[d]-LB[d]) : lUB - lLB
                swarm[p].v[d] = -r + 2*r*rand()
            end
        end
    end

    return nothing
end