text:
module ArrayInterfaceCore

using LinearAlgebra
using LinearAlgebra: AbstractTriangular
using SparseArrays
using SuiteSparse

@static if isdefined(Base, Symbol("@assume_effects"))
    using Base: @assume_effects
else
    macro assume_effects(_, ex)
        Base.@pure ex
    end
end

@assume_effects :total __parameterless_type(T) = Base.typename(T).wrapper
parameterless_type(x) = parameterless_type(typeof(x))
parameterless_type(x::Type) = __parameterless_type(x)

const VecAdjTrans{T,V<:AbstractVector{T}} = Union{Transpose{T,V},Adjoint{T,V}}
const MatAdjTrans{T,M<:AbstractMatrix{T}} = Union{Transpose{T,M},Adjoint{T,M}}
const UpTri{T,M} = Union{UpperTriangular{T,M},UnitUpperTriangular{T,M}}
const LoTri{T,M} = Union{LowerTriangular{T,M},UnitLowerTriangular{T,M}}

"""
    parent_type(::Type{T}) -> Type

Returns the parent array that type `T` wraps.
"""
parent_type(x) = parent_type(typeof(x))
parent_type(::Type{Symmetric{T,S}}) where {T,S} = S
parent_type(::Type{<:AbstractTriangular{T,S}}) where {T,S} = S
parent_type(@nospecialize T::Type{<:PermutedDimsArray}) = fieldtype(T, :parent)
parent_type(@nospecialize T::Type{<:Adjoint}) = fieldtype(T, :parent)
parent_type(@nospecialize T::Type{<:Transpose}) = fieldtype(T, :parent)
parent_type(@nospecialize T::Type{<:SubArray}) = fieldtype(T, :parent)
parent_type(@nospecialize T::Type{<:Base.ReinterpretArray}) = fieldtype(T, :parent)
parent_type(@nospecialize T::Type{<:Base.ReshapedArray}) = fieldtype(T, :parent)
parent_type(@nospecialize T::Type{<:Union{Base.Slice,Base.IdentityUnitRange}}) = fieldtype(T, :indices)
parent_type(::Type{Diagonal{T,V}}) where {T,V} = V
parent_type(T::Type) = T

"""
    buffer(x)

Return the buffer data that `x` points to. Unlike `parent(x::AbstractArray)`, `buffer(x)`
may not return another array type.
"""
buffer(x) = parent(x)
buffer(x::SparseMatrixCSC) = getfield(x, :nzval)
buffer(x::SparseVector) = getfield(x, :nzval)
buffer(@nospecialize x::Union{Base.Slice,Base.IdentityUnitRange}) = getfield(x, :indices)

"""
    is_forwarding_wrapper(::Type{T}) -> Bool

Returns `true` if the type `T` wraps another data type and does not alter any of its
standard interface. For example, if `T` were an array then its size, indices, and elements
would all be equivalent to its wrapped data.
"""
is_forwarding_wrapper(T::Type) = false
is_forwarding_wrapper(@nospecialize T::Type{<:Base.Slice}) = true
is_forwarding_wrapper(@nospecialize x) = is_forwarding_wrapper(typeof(x))

"""
    GetIndex(buffer) = GetIndex{true}(buffer)
    GetIndex{check}(buffer) -> g

Wraps an indexable buffer in a function type that is indexed when called, so that `g(inds..)`
is equivalent to `buffer[inds...]`. If `check` is `false`, then all indexing arguments are
considered in-bounds. The default value for `check` is `true`, requiring bounds checking for
each index.

!!! Warning
    Passing `false` as `check` may result in incorrect results/crashes/corruption for
    out-of-bounds indices, similar to inappropriate use of `@inbounds`. The user is
    responsible for ensuring this is correctly used.

# Examples

```julia
julia> ArrayInterfaceCore.GetIndex(1:10)[3]
3

julia> ArrayInterfaceCore.GetIndex{false}(1:10)[11]  # shouldn't be in-bounds
11

```

"""
struct GetIndex{CB,B} <: Function
    buffer::B

    GetIndex{true,B}(b) where {B} = new{true,B}(b)
    GetIndex{false,B}(b) where {B} = new{false,B}(b)
    GetIndex{check}(b::B) where {check,B} = GetIndex{check,B}(b)
    GetIndex(b) = GetIndex{true}(b)
end

buffer(g::GetIndex) = getfield(g, :buffer)

Base.@propagate_inbounds @inline (g::GetIndex{true})(inds...) = buffer(g)[inds...]
@inline (g::GetIndex{false})(inds...) = @inbounds(buffer(g)[inds...])

"""
    can_change_size(::Type{T}) -> Bool

Returns `true` if the Base.size of `T` can change, in which case operations
such as `pop!` and `popfirst!` are available for collections of type `T`.
"""
can_change_size(x) = can_change_size(typeof(x))
function can_change_size(::Type{T}) where {T}
    is_forwarding_wrapper(T) ? can_change_size(parent_type(T)) : false
end
can_change_size(::Type{<:Vector}) = true
can_change_size(::Type{<:AbstractDict}) = true
can_change_size(::Type{<:Base.ImmutableDict}) = false

function ismutable end

"""
    ismutable(::Type{T}) -> Bool

Query whether instances of type `T` are mutable or not, see
https://github.com/JuliaDiffEq/RecursiveArrayTools.jl/issues/19.
"""
ismutable(x) = ismutable(typeof(x))
function ismutable(::Type{T}) where {T<:AbstractArray}
    if parent_type(T) <: T
        return true
    else
        return ismutable(parent_type(T))
    end
end
ismutable(::Type{<:AbstractRange}) = false
ismutable(::Type{<:AbstractDict}) = true
ismutable(::Type{<:Base.ImmutableDict}) = false
ismutable(::Type{BigFloat}) = false
ismutable(::Type{BigInt}) = false
function ismutable(::Type{T}) where {T}
    if parent_type(T) <: T
        @static if VERSION ≥ v"1.7.0-DEV.1208"
            return Base.ismutabletype(T)
        else
            return T.mutable
        end
    else
        return ismutable(parent_type(T))
    end
end

# Piracy
function Base.setindex(x::AbstractArray, v, i...)
    _x = Base.copymutable(x)
    _x[i...] = v
    return _x
end

function Base.setindex(x::AbstractVector, v, i::Int)
    n = length(x)
    x .* (i .!== 1:n) .+ v .* (i .== 1:n)
end

function Base.setindex(x::AbstractMatrix, v, i::Int, j::Int)
    n, m = Base.size(x)
    x .* (i .!== 1:n) .* (j .!== i:m)' .+ v .* (i .== 1:n) .* (j .== i:m)'
end

"""
    can_setindex(::Type{T}) -> Bool

Query whether a type can use `setindex!`.
"""
can_setindex(x) = can_setindex(typeof(x))
can_setindex(T::Type) = is_forwarding_wrapper(T) ? can_setindex(parent_type(T)) : true
can_setindex(@nospecialize T::Type{<:AbstractRange}) = false
can_setindex(::Type{<:AbstractDict}) = true
can_setindex(::Type{<:Base.ImmutableDict}) = false
can_setindex(@nospecialize T::Type{<:Tuple}) = false
can_setindex(@nospecialize T::Type{<:NamedTuple}) = false
can_setindex(::Type{<:Base.Iterators.Pairs{<:Any,<:Any,P}}) where {P} = can_setindex(P)

"""
    aos_to_soa(x)

Converts an array of structs formulation to a struct of array.
"""
aos_to_soa(x) = x

"""
    isstructured(::Type{T}) -> Bool

Query whether a type is a representation of a structured matrix.
"""
isstructured(x) = isstructured(typeof(x))
isstructured(::Type) = false
isstructured(::Type{<:Symmetric}) = true
isstructured(::Type{<:Hermitian}) = true
isstructured(::Type{<:UpperTriangular}) = true
isstructured(::Type{<:LowerTriangular}) = true
isstructured(::Type{<:Tridiagonal}) = true
isstructured(::Type{<:SymTridiagonal}) = true
isstructured(::Type{<:Bidiagonal}) = true
isstructured(::Type{<:Diagonal}) = true

"""
    has_sparsestruct(x::AbstractArray) -> Bool

Determine whether `findstructralnz` accepts the parameter `x`.
"""
has_sparsestruct(x) = has_sparsestruct(typeof(x))
has_sparsestruct(::Type) = false
has_sparsestruct(::Type{<:AbstractArray}) = false
has_sparsestruct(::Type{<:SparseMatrixCSC}) = true
has_sparsestruct(::Type{<:Diagonal}) = true
has_sparsestruct(::Type{<:Bidiagonal}) = true
has_sparsestruct(::Type{<:Tridiagonal}) = true
has_sparsestruct(::Type{<:SymTridiagonal}) = true

"""
    issingular(A::AbstractMatrix) -> Bool

Determine whether a given abstract matrix is singular.
"""
issingular(A::AbstractMatrix) = issingular(Matrix(A))
issingular(A::AbstractSparseMatrix) = !issuccess(lu(A, check=false))
issingular(A::Matrix) = !issuccess(lu(A, check=false))
issingular(A::UniformScaling) = A.λ == 0
issingular(A::Diagonal) = any(iszero, A.diag)
issingular(A::Bidiagonal) = any(iszero, A.dv)
issingular(A::SymTridiagonal) = diaganyzero(ldlt(A).data)
issingular(A::Tridiagonal) = !issuccess(lu(A, check=false))
issingular(A::Union{Hermitian,Symmetric}) = diaganyzero(bunchkaufman(A, check=false).LD)
issingular(A::Union{LowerTriangular,UpperTriangular}) = diaganyzero(A.data)
issingular(A::Union{UnitLowerTriangular,UnitUpperTriangular}) = false
issingular(A::Union{Adjoint,Transpose}) = issingular(parent(A))
diaganyzero(A) = any(iszero, view(A, diagind(A)))

"""
    findstructralnz(x::AbstractArray)

Return: (I,J) #indexable objects
Find sparsity pattern of special matrices, the same as the first two elements of findnz(::SparseMatrixCSC).
"""
function findstructralnz(x::Diagonal)
    n = Base.size(x, 1)
    (1:n, 1:n)
end

function findstructralnz(x::Bidiagonal)
    n = Base.size(x, 1)
    isup = x.uplo == 'U' ? true : false
    rowind = BidiagonalIndex(n + n - 1, isup)
    colind = BidiagonalIndex(n + n - 1, !isup)
    (rowind, colind)
end

function findstructralnz(x::Union{Tridiagonal,SymTridiagonal})
    n = Base.size(x, 1)
    rowind = TridiagonalIndex(n + n - 1 + n - 1, n, true)
    colind = TridiagonalIndex(n + n - 1 + n - 1, n, false)
    (rowind, colind)
end

function findstructralnz(x::SparseMatrixCSC)
    rowind, colind, _ = findnz(x)
    (rowind, colind)
end

abstract type ColoringAlgorithm end

"""
    fast_matrix_colors(A)

Query whether a matrix has a fast algorithm for getting the structural
colors of the matrix.
"""
fast_matrix_colors(A) = false
fast_matrix_colors(A::AbstractArray) = fast_matrix_colors(typeof(A))
fast_matrix_colors(A::Type{<:Union{Diagonal,Bidiagonal,Tridiagonal,SymTridiagonal}}) = true

"""
    matrix_colors(A::Union{Array,UpperTriangular,LowerTriangular})

The color vector for dense matrix and triangular matrix is simply
`[1,2,3,..., Base.size(A,2)]`.
"""
function matrix_colors(A::Union{Array,UpperTriangular,LowerTriangular})
    eachindex(1:Base.size(A, 2)) # Vector Base.size matches number of rows
end
matrix_colors(A::Diagonal) = fill(1, Base.size(A, 2))
matrix_colors(A::Bidiagonal) = _cycle(1:2, Base.size(A, 2))
matrix_colors(A::Union{Tridiagonal,SymTridiagonal}) = _cycle(1:3, Base.size(A, 2))
_cycle(repetend, len) = repeat(repetend, div(len, length(repetend)) + 1)[1:len]

"""
  lu_instance(A) -> lu_factorization_instance

Returns an instance of the LU factorization object with the correct type
cheaply.
"""
function lu_instance(A::Matrix{T}) where {T}
    noUnitT = typeof(zero(T))
    luT = LinearAlgebra.lutype(noUnitT)
    ipiv = Vector{LinearAlgebra.BlasInt}(undef, 0)
    info = zero(LinearAlgebra.BlasInt)
    return LU{luT}(similar(A, 0, 0), ipiv, info)
end
function lu_instance(jac_prototype::SparseMatrixCSC)
    SuiteSparse.UMFPACK.UmfpackLU(
        Ptr{Cvoid}(),
        Ptr{Cvoid}(),
        1,
        1,
        jac_prototype.colptr[1:1],
        jac_prototype.rowval[1:1],
        jac_prototype.nzval[1:1],
        0,
    )
end

"""
  lu_instance(a::Number) -> a

Returns the number.
"""
lu_instance(a::Number) = a

"""
    lu_instance(a::Any) -> lu(a, check=false)

Returns the number.
"""
lu_instance(a::Any) = lu(a, check=false)

"""
    safevec(v)

It is a form of `vec` which is safe for all values in vector spaces, i.e., if it
is already a vector, like an AbstractVector or Number, it will return said
AbstractVector or Number.
"""
safevec(v) = vec(v)
safevec(v::Number) = v
safevec(v::AbstractVector) = v

"""
    zeromatrix(u::AbstractVector)

Creates the zero'd matrix version of `u`. Note that this is unique because
`similar(u,length(u),length(u))` returns a mutable type, so it is not type-matching,
while `fill(zero(eltype(u)),length(u),length(u))` doesn't match the array type,
i.e., you'll get a CPU array from a GPU array. The generic fallback is
`u .* u' .* false`, which works on a surprising number of types, but can be broken
with weird (recursive) broadcast overloads. For higher-order tensors, this
returns the matrix linear operator type which acts on the `vec` of the array.
"""
function zeromatrix(u)
    x = safevec(u)
    x .* x' .* false
end

# Reduces compile time burdens
function zeromatrix(u::Array{T}) where {T}
    out = Matrix{T}(undef, length(u), length(u))
    fill!(out, false)
end

"""
    restructure(x,y)

Restructures the object `y` into a shape of `x`, keeping its values intact. For
simple objects like an `Array`, this simply amounts to a reshape. However, for
more complex objects such as an `ArrayPartition`, not all of the structural
information is adequately contained in the type for standard tools to work. In
these cases, `restructure` gives a way to convert for example an `Array` into
a matching `ArrayPartition`.
"""
function restructure(x, y)
    out = similar(x, eltype(y))
    vec(out) .= vec(y)
    out
end

function restructure(x::Array, y)
    reshape(convert(Array, y), Base.size(x)...)
end

abstract type AbstractDevice end
abstract type AbstractCPU <: AbstractDevice end
struct CPUPointer <: AbstractCPU end
struct CPUTuple <: AbstractCPU end
struct CheckParent end
struct CPUIndex <: AbstractCPU end
struct GPU <: AbstractDevice end

"""
    can_avx(f) -> Bool

Returns `true` if the function `f` is guaranteed to be compatible with
`LoopVectorization.@avx` for supported element and array types. While a return
value of `false` does not indicate the function isn't supported, this allows a
library to conservatively apply `@avx` only when it is known to be safe to do so.

```julia
function mymap!(f, y, args...)
    if can_avx(f)
        @avx @. y = f(args...)
    else
        @. y = f(args...)
    end
end
```
"""
can_avx(::Any) = false

"""
    fast_scalar_indexing(::Type{T}) -> Bool

Query whether an array type has fast scalar indexing.
"""
fast_scalar_indexing(x) = fast_scalar_indexing(typeof(x))
fast_scalar_indexing(::Type) = true
fast_scalar_indexing(::Type{<:LinearAlgebra.AbstractQ}) = false
fast_scalar_indexing(::Type{<:LinearAlgebra.LQPackedQ}) = false

"""
    allowed_getindex(x,i...)

A scalar `getindex` which is always allowed.
"""
allowed_getindex(x, i...) = x[i...]

"""
    allowed_setindex!(x,v,i...)

A scalar `setindex!` which is always allowed.
"""
allowed_setindex!(x, v, i...) = Base.setindex!(x, v, i...)

"""
    ArrayIndex{N}

Subtypes of `ArrayIndex` represent series of transformations for a provided index to some
buffer which is typically accomplished with square brackets (e.g., `buffer[index[inds...]]`).
The only behavior that is required of a subtype of `ArrayIndex` is the ability to transform
individual index elements (i.e. not collections). This does not guarantee bounds checking or
the ability to iterate (although additional functionallity may be provided for specific
types).
"""
abstract type ArrayIndex{N} end

const MatrixIndex = ArrayIndex{2}

const VectorIndex = ArrayIndex{1}

Base.ndims(::Type{<:ArrayIndex{N}}) where {N} = N

struct BidiagonalIndex <: MatrixIndex
    count::Int
    isup::Bool
end

struct TridiagonalIndex <: MatrixIndex
    count::Int# count==nsize+nsize-1+nsize-1
    nsize::Int
    isrow::Bool
end

Base.firstindex(i::Union{BidiagonalIndex,TridiagonalIndex}) = 1
Base.lastindex(i::Union{BidiagonalIndex,TridiagonalIndex}) = i.count
Base.length(i::Union{BidiagonalIndex,TridiagonalIndex}) = lastindex(i)

Base.@propagate_inbounds function Base.getindex(ind::BidiagonalIndex, i::Int)
    @boundscheck 1 <= i <= ind.count || throw(BoundsError(ind, i))
    if ind.isup
        ii = i + 1
    else
        ii = i + 1 + 1
    end
    convert(Int, floor(ii / 2))
end

Base.@propagate_inbounds function Base.getindex(ind::TridiagonalIndex, i::Int)
    @boundscheck 1 <= i <= ind.count || throw(BoundsError(ind, i))
    offsetu = ind.isrow ? 0 : 1
    offsetl = ind.isrow ? 1 : 0
    if 1 <= i <= ind.nsize
        return i
    elseif ind.nsize < i <= ind.nsize + ind.nsize - 1
        return i - ind.nsize + offsetu
    else
        return i - (ind.nsize + ind.nsize - 1) + offsetl
    end
end

_cartesian_index(i::Tuple{Vararg{Int}}) = CartesianIndex(i)
_cartesian_index(::Any) = nothing

"""
    known_first(::Type{T}) -> Union{Int,Nothing}

If `first` of an instance of type `T` is known at compile time, return it.
Otherwise, return `nothing`.

```julia
julia> ArrayInterface.known_first(typeof(1:4))
nothing

julia> ArrayInterface.known_first(typeof(Base.OneTo(4)))
1
```
"""
known_first(x) = known_first(typeof(x))
known_first(T::Type) = is_forwarding_wrapper(T) ? known_first(parent_type(T)) : nothing
known_first(::Type{<:Base.OneTo}) = 1
known_first(@nospecialize T::Type{<:LinearIndices}) = 1
known_first(@nospecialize T::Type{<:Base.IdentityUnitRange}) = known_first(parent_type(T))
function known_first(::Type{<:CartesianIndices{N,R}}) where {N,R}
    _cartesian_index(ntuple(i -> known_first(R.parameters[i]), Val(N)))
end

"""
    known_last(::Type{T}) -> Union{Int,Nothing}

If `last` of an instance of type `T` is known at compile time, return it.
Otherwise, return `nothing`.

```julia
julia> ArrayInterfaceCore.known_last(typeof(1:4))
nothing

julia> ArrayInterfaceCore.known_first(typeof(static(1):static(4)))
4

```
"""
known_last(x) = known_last(typeof(x))
known_last(T::Type) = is_forwarding_wrapper(T) ? known_last(parent_type(T)) : nothing
function known_last(::Type{<:CartesianIndices{N,R}}) where {N,R}
    _cartesian_index(ntuple(i -> known_last(R.parameters[i]), Val(N)))
end

"""
    known_step(::Type{T}) -> Union{Int,Nothing}

If `step` of an instance of type `T` is known at compile time, return it.
Otherwise, return `nothing`.

```julia
julia> ArrayInterface.known_step(typeof(1:2:8))
nothing

julia> ArrayInterface.known_step(typeof(1:4))
1

```
"""
known_step(x) = known_step(typeof(x))
known_step(T::Type) = is_forwarding_wrapper(T) ? known_step(parent_type(T)) : nothing
known_step(@nospecialize T::Type{<:AbstractUnitRange}) = 1

"""
    is_splat_index(::Type{T}) -> Bool
Returns `static(true)` if `T` is a type that splats across multiple dimensions.
"""
is_splat_index(T::Type) = false
is_splat_index(@nospecialize(x)) = is_splat_index(typeof(x))

"""
    ndims_index(::Type{I}) -> Int

Returns the number of dimension that an instance of `I` maps to when indexing. For example,
`CartesianIndex{3}` maps to 3 dimensions. If this method is not explicitly defined, then `1`
is returned.
"""
ndims_index(::Type{<:Base.AbstractCartesianIndex{N}}) where {N} = N
# preserve CartesianIndices{0} as they consume a dimension.
ndims_index(::Type{CartesianIndices{0,Tuple{}}}) = 1
ndims_index(@nospecialize T::Type{<:AbstractArray{Bool}}) = ndims(T)
ndims_index(@nospecialize T::Type{<:AbstractArray}) = ndims_index(eltype(T))
ndims_index(@nospecialize T::Type{<:Base.LogicalIndex}) = ndims(fieldtype(T, :mask))
ndims_index(T::Type) = 1
ndims_index(@nospecialize(i)) = ndims_index(typeof(i))

"""
    ndims_shape(::Type{I}) -> Union{Int,Tuple{Vararg{Int}}}

Returns the number of dimension that are represented in shape of the returned array when
indexing with an instance of `I`.
"""
ndims_shape(T::DataType) = ndims_index(T)
ndims_shape(::Type{Colon}) = 1
ndims_shape(T::Type{<:Base.AbstractCartesianIndex{N}}) where {N} = ntuple(zero, Val{N}())
ndims_shape(@nospecialize T::Type{<:CartesianIndices}) = ntuple(one, Val{ndims(T)}())
ndims_shape(@nospecialize T::Type{<:Number}) = 0
ndims_shape(@nospecialize T::Type{<:AbstractArray}) = ndims(T)
ndims_shape(x) = ndims_shape(typeof(x))

"""
    IndicesInfo(T::Type{<:Tuple}) -> IndicesInfo{NI,NS,IS}()

Provides basic trait information for each index type in in the tuple `T`. `NI`, `NS`, and
`IS` are tuples of [`ndims_index`](@ref), [`ndims_shape`](@ref), and
[`is_splat_index`](@ref) (respectively) for each field of `T`.
"""
struct IndicesInfo{NI,NS,IS} end
IndicesInfo(@nospecialize x::Tuple) = IndicesInfo(typeof(x))
@generated function IndicesInfo(::Type{T}) where {T<:Tuple}
    NI = Expr(:tuple)
    NS = Expr(:tuple)
    IS = Expr(:tuple)
    for i in 1:fieldcount(T)
        T_i = fieldtype(T, i)
        push!(NI.args, :(ndims_index($(T_i))))
        push!(NS.args, :(ndims_shape($(T_i))))
        push!(IS.args, :(is_splat_index($(T_i))))
    end
    Expr(:block, Expr(:meta, :inline), :(IndicesInfo{$(NI),$(NS),$(IS)}()))
end

"""
    instances_do_not_alias(::Type{T}) -> Bool

Is it safe to `ivdep` arrays containing elements of type `T`?
That is, would it be safe to write to an array full of `T` in parallel?
This is not true for `mutable struct`s in general, where editing one index
could edit other indices.
That is, it is not safe when different instances may alias the same memory.
"""
instances_do_not_alias(::Type{T}) where {T} = Base.isbitstype(T)

"""
    indices_do_not_alias(::Type{T<:AbstractArray}) -> Bool

Is it safe to `ivdep` arrays of type `T`?
That is, would it be safe to write to an array of type `T` in parallel?
Examples where this is not true are `BitArray`s or `view(rand(6), [1,2,3,1,2,3])`.
That is, it is not safe whenever different indices may alias the same memory.
"""
indices_do_not_alias(::Type) = false
indices_do_not_alias(::Type{A}) where {T, A<:Base.StridedArray{T}} = instances_do_not_alias(T)
indices_do_not_alias(::Type{Adjoint{T,A}}) where {T, A <: AbstractArray{T}} = indices_do_not_alias(A)
indices_do_not_alias(::Type{Transpose{T,A}}) where {T, A <: AbstractArray{T}} = indices_do_not_alias(A)
indices_do_not_alias(::Type{<:SubArray{<:Any,<:Any,A,I}}) where {
  A,I<:Tuple{Vararg{Union{Integer, UnitRange, Base.ReshapedUnitRange, Base.AbstractCartesianIndex}}}} = indices_do_not_alias(A)

end # module


text:
include("abstract_buffer.jl")
include("circular_array_buffer.jl")
include("circular_turn_buffer.jl")
include("episode_turn_buffer.jl")

ext: jl


text:
reload("SSH")
isdefined(:sock) && close(sock)
gc()
port = length(ARGS) >= 1 ? parse(Int, ARGS[1]) : 10000
sock = listen(port)

function open_fake_pty()
    const O_RDWR = Base.Filesystem.JL_O_RDWR
    const O_NOCTTY = Base.Filesystem.JL_O_NOCTTY

    fdm = ccall(:posix_openpt, Cint, (Cint,), O_RDWR|O_NOCTTY)
    fdm == -1 && error("Failed to open PTY master")
    rc = ccall(:grantpt, Cint, (Cint,), fdm)
    rc != 0 && error("grantpt failed")
    rc = ccall(:unlockpt, Cint, (Cint,), fdm)
    rc != 0 && error("unlockpt")

    fds = ccall(:open, Cint, (Ptr{UInt8}, Cint),
        ccall(:ptsname, Ptr{UInt8}, (Cint,), fdm), O_RDWR|O_NOCTTY)

    # slave
    slave   = RawFD(fds)
    master = Base.TTY(RawFD(fdm); readable = true)
    slave, master, RawFD(fdm)
end

while true
    client = accept(sock)
    @async begin
        session = connect(SSH.Session, client; client = false)
        SSH.negotiate_algorithms!(session)
        SSH.server_dh_kex!(session,
            joinpath(ENV["HOME"],".ssh","id_rsa.pub"),
            joinpath(ENV["HOME"],".ssh","id_rsa"))
        SSH.wait_for_userauth(session, publickey = function (username, algorithm, blob)
            return true
        end)
        SSH.perform_ssh_connection(session) do kind, channel
            if kind == "session"
                c = Condition()
                local encoded_termios
                SSH.on_channel_request(channel) do kind, packet
                    want_reply = read(packet, UInt8) != 0
                    @show kind
                    if kind == "pty-req"
                        TERM_var = SSH.read_string(packet)
                        termwidthchars = bswap(read(packet, UInt32))
                        termheightchars = bswap(read(packet, UInt32))
                        termwidthpixs = bswap(read(packet, UInt32))
                        termheightpixs = bswap(read(packet, UInt32))
                        encoded_modes = SSH.read_string(packet)
                        encoded_termios = IOBuffer(encoded_modes)
                        notify(c)
                    elseif kind == "signal"
                        SSH.disconnect(channel.session)
                    end
                    want_reply
                end
                open(channel)
                TIOCSCTTY_str = """
                    Libc.systemerror("ioctl",
                    0 != ccall(:ioctl, Cint, (Cint, Cint, Int64), 0,
                    (is_bsd() || is_apple()) ? 0x20007461 : is_linux() ? 0x540E :
                    error("Fill in TIOCSCTTY for this OS here"), 0))
                """
                cmd = """
                    $TIOCSCTTY_str
                    try; run(`lua $(joinpath(ENV["HOME"],"termtris/termtris.lua"))`); end
                """
                @async begin
                    wait(c)
                    slave, master, masterfd = open_fake_pty()
                    new_termios = Ref{SSH.termios}()
                    systemerror("tcgetattr",
                        -1 == ccall(:tcgetattr, Cint, (Cint, Ptr{Void}), slave, new_termios))
                    new_termios[] = SSH.decode_modes(encoded_termios, new_termios[])
                    systemerror("tcsetattr",
                        -1 == ccall(:tcsetattr, Cint, (Cint, Cint, Ptr{Void}), slave, 0, new_termios))
                    p = spawn(detach(`$(Base.julia_cmd()) -e $cmd`), slave, slave, slave)
                    @async while true
                        write(channel, readavailable(master))
                    end
                    @async while true
                        data = readavailable(channel)
                        @show data
                        write(master, data)
                    end
                    wait(p)
                    SSH.disconnect(session)
                end
            end
        end
    end
    println("Done")
end


text:
# Advent of Code 2020 - Day 15

function spoken(list, turn)
    spoken_init = parse.(Int,(split(list, ',')))
    spoken_list = fill(0, turn -1)
    for j=1:length(spoken_init)-1
        spoken_list[spoken_init[j]+1] += j
    end
    last_spoken = spoken_init[end]
    for i=length(spoken_init):turn-1
        if spoken_list[last_spoken+1] == 0
            spoken_list[last_spoken+1] = i
            last_spoken = 0
        else
            difference = i - spoken_list[last_spoken+1]
            spoken_list[last_spoken+1] = i
            last_spoken = difference
        end
    end
    return last_spoken
end

# The next Int32 version reduces allocations by half and time
function spoken32(list, turn)
    spoken_init = parse.(Int32,(split(list, ',')))
    spoken_list = fill(Int32(0), turn -1)
    for j=1:length(spoken_init)-1
        spoken_list[spoken_init[j]+1] += Int32(j)
    end
    last_spoken = spoken_init[end]
    for i=length(spoken_init):turn-1
        if spoken_list[last_spoken+1] == 0
            spoken_list[last_spoken+1] = Int32(i)
            last_spoken = 0
        else
            difference = i - spoken_list[last_spoken+1]
            spoken_list[last_spoken+1] = Int32(i)
            last_spoken = difference
        end
    end
    return last_spoken
end

@show spoken("0,3,6", 2020) == 436
@show spoken("1,3,2", 2020) == 1
@show spoken("2,1,3", 2020) == 10
@show spoken("1,2,3", 2020) == 27
@show spoken("2,3,1", 2020) == 78
@show spoken("3,2,1", 2020) == 438
@show spoken("3,1,2", 2020) == 1836
@show spoken("10,16,6,0,1,17", 2020) == 412

@show spoken("0,3,6", 30000000) == 175594
@show spoken("1,3,2", 30000000) == 2578
@show spoken("2,1,3", 30000000) == 3544142
@show spoken("1,2,3", 30000000) == 261214
@show spoken("2,3,1", 30000000) == 6895259
@show spoken("3,2,1", 30000000) == 18
@show spoken("3,1,2", 30000000) == 362
@show spoken("10,16,6,0,1,17", 30000000) == 243
@show spoken32("10,16,6,0,1,17", 30000000) == 243

#=
```julia
julia> @btime spoken("10,16,6,0,1,17", 30000000)
  1.568 s (6 allocations: 228.88 MiB)
243

julia> @btime spoken32("10,16,6,0,1,17", 30000000)
  1.328 s (6 allocations: 114.44 MiB)
243
```
 =#


text:
# Miscellaneous helper macros and functions
#

# Helper macro to iterate over the DG grid. Generates the needed loops
# and indices: `eh`, `ev`, `e`, `k,`, `j`, `i`, `ijk`.
macro visitQ(nhorzelem, nvertelem, Nqk, Nq, expr)
    return esc(
        quote
            for eh in 1:($nhorzelem)
                for ev in 1:($nvertelem)
                    e = ev + (eh - 1) * $nvertelem
                    for k in 1:($Nqk)
                        for j in 1:($Nq)
                            for i in 1:($Nq)
                                ijk = i + $Nq * ((j - 1) + $Nq * (k - 1))
                                $expr
                            end
                        end
                    end
                end
            end
        end,
    )
end

# Helper macro to iterate over a 3D array. Used for walking the
# interpolated grid.
macro visitI(nlong, nlat, nlevel, expr)
    return esc(quote
        for lo in 1:($nlong)
            for la in 1:($nlat)
                for le in 1:($nlevel)
                    $expr
                end
            end
        end
    end)
end

# Helpers to extract data from the various state arrays
function extract_state_conservative(dg, state_conservative, ijk, e)
    bl = dg.balance_law
    FT = eltype(state_conservative)
    num_state_conservative = number_state_conservative(bl, FT)
    local_state_conservative = MArray{Tuple{num_state_conservative}, FT}(undef)
    for s in 1:num_state_conservative
        local_state_conservative[s] = state_conservative[ijk, s, e]
    end
    return Vars{vars_state_conservative(bl, FT)}(local_state_conservative)
end
function extract_state_auxiliary(dg, state_auxiliary, ijk, e)
    bl = dg.balance_law
    FT = eltype(state_auxiliary)
    num_state_auxiliary = number_state_auxiliary(bl, FT)
    local_state_auxiliary = MArray{Tuple{num_state_auxiliary}, FT}(undef)
    for s in 1:num_state_auxiliary
        local_state_auxiliary[s] = state_auxiliary[ijk, s, e]
    end
    return Vars{vars_state_auxiliary(bl, FT)}(local_state_auxiliary)
end
function extract_state_gradient_flux(dg, state_gradient_flux, ijk, e)
    bl = dg.balance_law
    FT = eltype(state_gradient_flux)
    num_state_gradient_flux = number_state_gradient_flux(bl, FT)
    local_state_gradient_flux =
        MArray{Tuple{num_state_gradient_flux}, FT}(undef)
    for s in 1:num_state_gradient_flux
        local_state_gradient_flux[s] = state_gradient_flux[ijk, s, e]
    end
    return Vars{vars_state_gradient_flux(bl, FT)}(local_state_gradient_flux)
end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
# This file was generated, do not modify it. # hide
@from_network CompositeModel(std=std_model, box=box_model,
                             ridge=ridge_model) <= ŷ;


text:
export fwi_objective

function fwi_objective(model_full::ModelRho, source::joData, dObs::joData, srcnum::Int64; options=Options())
# Setup time-domain linear or nonlinear foward and adjoint modeling and interface to OPESCI/devito

	# Load full geometry for out-of-core geometry containers
	typeof(dObs.geometry) == GeometryOOC && (dObs.geometry = Geometry(dObs.geometry))
	typeof(source.geometry) == GeometryOOC && (source.geometry = Geometry(source.geometry))

	# for 3D modeling, limit model to area with sources/receivers
	if options.limit_m == true && model_full.n[3] > 1	# only supported for 3D
		model = deepcopy(model_full)
		model = limit_model_to_receiver_area(source.geometry,dObs.geometry,model,options.buffer_size)
	else
		model = model_full
	end

 	# Dimensions for arrays initialization
	if length(model.n) == 2
		nx, nz = model.n
		damp = zeros(Float32,nx, nz)
		dimsflip = (2,1)
		dimsfull = nz+model.nb*2, nx+model.nb*2	# dimensions for devito (z,x)
	else
		nx, ny, nz = model.n
		damp = zeros(Float32,nx, ny, nz)
		dimsflip = (3,2,1)
		dimsfull = nz+model.nb*2, ny+model.nb*2, nx+model.nb*2	# dimensions for devito (z,y,x)
	end

	# Set up model structure
	modelPy = ct.IGrid(model.o, model.d, model.n, sqrt.(1./model.m), rho=model.rho, nbpml=model.nb)
	dtComp = modelPy[:critical_dt]
	h = modelPy[:get_spacing]()

	# Initiate source PyObject
	src = ct.IShot()
	tmaxSrc = source.geometry.t[1]
	src[:set_time_axis](dtComp, tmaxSrc)
	numTracesSrc = length(source.geometry.zloc[1])

	# Initiate receiver PyObject
	rec = ct.IShot()
	tmaxRec = dObs.geometry.t[1]
	rec[:set_time_axis](dtComp, tmaxRec)
	numTracesRec = length(dObs.geometry.zloc[1])

	# Extrapolate input data to computational grid
	qIn = (src[:reinterpolateD](reshape(source.data[1], source.geometry.nt[1], numTracesSrc), source.geometry.dt[1], dtComp))'
	if typeof(dObs.data[1]) == SeisIO.SeisCon
		data = convert(Array{Float32,2},dObs.data[1][1].data)
		dObs = joData(dObs.geometry,data)
	end
	dObserved = (rec[:reinterpolateD](reshape(dObs.data[1], dObs.geometry.nt[1], numTracesRec), dObs.geometry.dt[1], dtComp))'

	ntComp = size(dObserved,2)
	ntSrc = Int(trunc(tmaxSrc/dtComp+1))
	ntRec = Int(trunc(tmaxRec/dtComp+1))
	isempty(options.save_rate) ? save_rate=dtComp : save_rate=options.save_rate

	# Set up source coordinates
	src_coords = setup_grid(source.geometry, model.n, model.o)
	src[:set_receiver_pos](src_coords)
	src[:set_shape](ntComp, numTracesSrc)

	# Set up receiver coordinates
	rec_coords = setup_grid(dObs.geometry, model.n, model.o)
	rec[:set_receiver_pos](rec_coords)
	rec[:set_shape](ntComp, numTracesRec)

	# Initiate acoustic modeling object
	Acoustic = ac.JAcoustic_cg(modelPy, rec, src, s_order=options.space_order, tsave=save_rate)

	# Initialize dampening array and extend velocity model
	damp = damp_boundary(add_pml(damp, model.nb), model)
	m = add_pml(model.m, model.nb)
	rho = add_pml(model.rho, model.nb)

	# Reshape from Julia to Python array ordering (x,y,z) -> (z,y,x)
	damp = permutedims(damp,dimsflip)
	m = permutedims(m,dimsflip)
	rho = permutedims(rho,dimsflip)

	# Forward modeling to generate synthetic data and background wavefields
	u = zeros(Float32,(prod(dimsfull), 3))
	nsave = trunc(Int,ntComp/(save_rate/dtComp) +1)
	up = zeros(Float32,(prod(dimsfull),nsave))
	dPredicted = zeros(Float32,rec[:shape])'
	op = Acoustic[:Forward](save=true, free_surface=options.free_surface)
	types, args = process_args(op, Dict(:damp=>damp, :m=>m, :rec=>dPredicted, :rec_coords=>rec_coords', :src=>qIn, :src_coords=>src_coords', :u=>u, :rho=>rho, :usave=>up))
	lib = op[:compile]
	status = ccallexpr(lib, Symbol(op[:name]), types, args)

	# zero-padd output if shot record has more samples than source and observed data
	if ntRec > ntComp
		dPredicted = [dPredicted zeros(size(dOut,1),ntRec - ntComp)]
	end

	# Data misfit
	argout1 = .5f0*norm(vec(dPredicted) - vec(dObserved),2)^2.f0

	# Backpropagation of data residual
	v = zeros(Float32,(prod(dimsfull), 3))
	residual = dPredicted - dObserved
	gradient = zeros(Float32, size(m))
	opg = Acoustic[:Gradient](free_surface=options.free_surface)
	types, args = process_args(opg, Dict(:damp=>damp, :m=>m, :rec=>residual, :rec_coords=>rec_coords', :u=>up, :grad=>gradient, :v=>v, :rho=>rho))
	lib = opg[:compile]
	status = ccallexpr(lib, Symbol(opg[:name]), types, args)
	argout2 = remove_pml(permutedims(reshape(gradient,dimsfull),dimsflip),model.nb)

	# Extend gradient back to original size
	if options.limit_m==true && model_full.n[3] > 1
		argout2 = extend_gradient(model_full,model,argout2)
	end
	u=[]; v=[]; gc();
	return [argout1; vec(argout2)]
end

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
# rotate3.jl
# time 3D image rotation using @thread and foreach

using Revise
using SPECTrecon
using BenchmarkTools
using Plots; default(markerstrokecolor=:auto)

nx, nz = 128, 100
T = Float32
image3 = randn(T, nx,nx,nz)
out1 = similar(image3)
out2 = similar(image3)

plans = plan_rotate(nx; T) # [Threads.nthreads()]

θ = π/6
imrotate!(out1, image3, θ, plans, 17) # foreach
imrotate!(out2, image3, θ, plans, :thread)
@assert out1 == out2

#@btime imrotate!($out1, $image3, $θ, $plans, 20) # 5.9 ms (888 allocations: 63.53 KiB)
#@btime imrotate!($out1, $image3, $θ, $plans) # 6.1 ms (802 allocations: 56.77 KiB)
#@btime imrotate!($out2, $image3, $θ, $plans, :thread) # 5.8 ms (41 allocations: 3.17 KiB)
nothing

if !@isdefined(elapse0)
    elapse0 = @belapsed imrotate!($out1, $image3, $θ, $plans, :thread)
end

if !@isdefined(elapse)
	ntask = 1:20
	elapse = zeros(length(ntask))
	for (i,nt) in enumerate(ntask)
		@show i, nt
    	elapse[i] = @belapsed imrotate!($out1, $image3, $θ, $plans, $nt)
	end
end

@show minimum(elapse), elapse0

ncore = Threads.nthreads()
i = findall(==(ncore), ntask)[1]
scatter(ntask, elapse, xlabel="ntask", ylabel="time", label="foreach",
ylim=(0.005,0.03),
xlim=(1,20),
ytick=[0.006,0.010,0.015,0.029],
xtick=2 .^ (0:5),
)
scatter!([ntask[i]], [elapse[i]], color=:red, label="@threads", title="rotate $nx × $nx × $nz")
#savefig("rotate.pdf")


text:
mutable struct generic_variable_t <: LCMCore.LCMType
  utime::Int64

  variablelabel::String
  variabletype::String

  datadescription::String
  datalength::Int32
  data::Vector{UInt8}
end

@lcmtypesetup(generic_variable_t,
  data => (datalength,),
)


text:
# This file was generated, do not modify it. # hide
r3 = x -> round(x, sigdigits=3)
r3(pi)


text:
using IteratorInterfaceExtensions, TableTraits

import IteratorInterfaceExtensions: isiterable, getiterator
import TableTraits: isiterabletable
import Base: length, eltype, iterate

# T stores the type info for a NamedTuple
# S stores the SamplePath type
struct SamplePathIterator{T,S}
    xw::S
end

# required methods for SamplePath
IteratorInterfaceExtensions.isiterable(xw::SamplePath) = true
TableTraits.isiterabletable(xw::SamplePath) = true

# TODO: how should we handle Configuration data?
function IteratorInterfaceExtensions.getiterator(xw::SamplePath{T}) where T <: Number
    t_type = eltype(xw.t)

    # format:
    # t X1 X2 X3 ...
    value_types = Type[]
    value_names = Symbol[]

    push!(value_types, t_type)
    push!(value_names, :t)

    # use the first record as a representative
    example = xw.u[1]

    # fill in the types and names for each component; should be the same
    # should be the same (for now)
    for k in 1:length(example)
        push!(value_types, eltype(example))
        push!(value_names, Symbol("X$(k)"))
    end

    # construct our iterator
    iter = SamplePathIterator{NamedTuple{(value_names...,), Tuple{value_types...}}, typeof(xw)}(xw)

    return iter
end

# iteration interface implementation for SamplePathIterator
Base.length(iter::SamplePathIterator) = length(iter.xw)
Base.eltype(iter::SamplePathIterator{T,S}) where {T,S} = T

@generated function Base.iterate(iter::SamplePathIterator{T,S}, state=1) where {T,S}
    record = []

    # add t value for the record
    push!(record, :(iter.xw.t[state]))

    # populate the record with values for each component
    n = length(fieldnames(T))
    for i in 1:n-1
        push!(record, :(iter.xw.u[state][$i]))
    end

    # function body
    # each iterate is a named tuple
    quote
        if state > length(iter)
            return nothing
        else
            return $(T)(($(record...),)), state+1
        end
    end
end

# workaround for broken integrations in IterableTables.jl
tablefy(xw::SamplePath) = IteratorInterfaceExtensions.getiterator(xw)

# T stores the type info for a NamedTuple
# S stores the SamplePath type
struct EnsembleIterator{T,S}
    x::S
    trials::Int
end

# required methods for SamplePath
IteratorInterfaceExtensions.isiterable(x::Ensemble) = true
TableTraits.isiterabletable(x::Ensemble) = true

# TODO: how should we handle Configuration data?
function IteratorInterfaceExtensions.getiterator(x::Ensemble{T}) where T <: Number
    t_type = eltype(x[1].t)

    # format:
    # trial t X1 X2 X3 ...
    value_types = Type[]
    value_names = Symbol[]

    # first column is iteration number
    push!(value_types, Int)
    push!(value_names, :trial)

    # second column is a timestamp
    push!(value_types, t_type)
    push!(value_names, :t)

    # use the first record as a representative
    example = x[1].u[1]

    # fill in the types and names for each component; should be the same
    # should be the same (for now)
    for k in 1:length(example)
        push!(value_types, eltype(example))
        push!(value_names, Symbol("X$(k)"))
    end

    # construct our iterator
    iter = EnsembleIterator{NamedTuple{(value_names...,), Tuple{value_types...}}, typeof(x)}(x, length(x))

    return iter
end

# iteration interface implementation for SamplePathIterator
Base.length(iter::EnsembleIterator) = sum(length, iter.x)
Base.eltype(iter::EnsembleIterator{T,S}) where {T,S} = T

@generated function Base.iterate(iter::EnsembleIterator{T,S}, state=(1,1)) where {T,S}
    record = []
    # add trial number
    push!(record, :(state[1]))

    # add t value for the record
    push!(record, :(iter.x[state[1]].t[state[2]]))

    # populate the record with values for each component
    n = length(fieldnames(T))

    for i in 1:n-2
        push!(record, :(iter.x[state[1]].u[state[2]][$i]))
    end

    # function body
    # each iterate is a named tuple
    quote
        if state[1] > iter.trials
            return nothing
        else
            # get number of states for current trajectory
            m = length(iter.x[state[1]])

            # if we've touched all the values for this trajectory,
            # move on to the next one
            if state[2] ≥ m
                next_state = (state[1]+1, 1)
            else
            # otherwise, advance to the next time point
                next_state = (state[1], state[2]+1)
            end

            return $(T)(($(record...),)), next_state
        end
    end
end

# workaround for broken integrations in IterableTables.jl
tablefy(x::Ensemble) = IteratorInterfaceExtensions.getiterator(x)


text:
"""
    MeshSelectionModule

Module for  selection of mesh entities.
"""
module MeshSelectionModule

__precompile__(true)

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import ..FESetModule: AbstractFESet, bfundpar, nodesperelem, manifdim, connasarray
import ..FENodeSetModule: FENodeSet, spacedim
import ..BoxModule: inflatebox!, initbox!, updatebox!, boxesoverlap
import ..FENodeToFEMapModule: FENodeToFEMap
import LinearAlgebra: norm, dot, cross
import Statistics: mean

"""
    connectednodes(fes::AbstractFESet)

Extract the node numbers of the nodes connected by given finite elements.

Extract the list of unique node numbers for the nodes that are connected
by the finite element set `fes`. Note that it is assumed that all the
FEs are of the same type (the same number of connected nodes by each
cell).
"""
function connectednodes(fes::AbstractFESet)
    return unique(connasarray(fes)[:]);
end

"""
    connectedelems(fes::AbstractFESet, node_list::FIntVec)

Extract the list of numbers for the fes  that are connected to given
nodes.
"""
function connectedelems(fes::AbstractFESet, node_list::FIntVec, nmax::FInt)
    f2fm = FENodeToFEMap(fes.conn, nmax)
    cg = zeros(FInt, count(fes)); # No elements are part of the group to begin with
    for j = node_list # for all nodes in the list
        for i = 1:length(f2fm.map[j])
            cg[f2fm.map[j][i]] = 1   # Mark element as being part of the group
        end
    end
    return findall(v -> (v != 0), cg)
end

"""
    selectnode(fens::FENodeSet; kwargs...)

Select nodes using some criterion.

# Arguments
- `v` = array of locations, one location per row
- `kwargs` = pairs of keyword argument/value

## Selection criteria

### box
```
nLx = vselect(fens.xyz, box = [0.0 Lx  0.0 0.0 0.0 0.0], inflate = Lx/1.0e5)
```

The keyword 'inflate' may be used to increase or decrease the extent of
the box (or the distance) to make sure some nodes which would be on the
boundary are either excluded or included.

### distance
```
list = selectnode(fens.xyz, distance=1.0+0.1/2^nref, from=[0. 0.],
        inflate=tolerance);
```

### plane
```
candidates = selectnode(fens, plane = [0.0 0.0 1.0 0.0], thickness = h/1000)
```
The keyword `plane` defines the plane by its normal (the first two or
three numbers) and its distance from the origin (the last number). Nodes
are selected they lie on the plane,  or near the plane within the
distance `thickness` from the plane. The normal is assumed to be of unit
length, if it isn't apply as such, it will be normalized internally.

### nearestto
Find the node nearest to the location given.
```
nh = selectnode(fens, nearestto = [R+Ro/2, 0.0, 0.0] )
```
"""
function selectnode(fens::FENodeSet; kwargs...)
    nodelist = vselect(fens.xyz; kwargs...)
    nodelist = dropdims(reshape(nodelist,1,length(nodelist)), dims=1);
    return nodelist
end

"""
    selectelem(fens::FENodeSet, fes::T; kwargs...) where {T<:AbstractFESet}

Select finite elements.

# Arguments
- `fens` = finite element node set
- `fes` = finite element set
- `kwargs` = keyword arguments to specify the selection criteria

## Selection criteria

### facing
Select all "boundary" elements that "face" a certain direction.
```
exteriorbfl = selectelem(fens, bdryfes, facing=true, direction=[1.0, 1.0, 0.0]);
```
or
```
exteriorbfl = selectelem(fens, bdryfes, facing=true, direction=dout, dotmin = 0.99);
```
where
```
function dout(xyz)
    return xyz/norm(xyz)
end
```
and `xyz` is the location of the centroid  of  a boundary element.
Here the finite element is considered "facing" in the given direction if the dot
product of its normal and the direction vector is greater than `dotmin`.
The default value for `dotmin` is 0.01 (this corresponds to  almost 90 degrees
between the normal to the finite element  and the given direction).

This selection method makes sense only for elements that are  surface-like (i. e.
for boundary mmeshes).

### label
Select elements based on their label.
```
rl1 = selectelem(fens, fes, label=1)
```

### box, distance
Select elements based on some criteria that their nodes satisfy.  See the
function `selectnode()`.

Example:
Select all  elements whose nodes are closer than `R+inflate` from the point `from`.
```
linner = selectelem(fens, bfes, distance = R, from = [0.0 0.0 0.0],
  inflate = tolerance)
```

Example:

```
exteriorbfl = selectelem(fens, bdryfes,
   box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);
```

### withnodes

Select elements whose nodes are in a given list of node numbers.

Example:
```julia
l = selectelem(fens, fes, withnodes = [13, 14])
```

### flood
Select all FEs connected together, starting from a given node.
Connections through a vertex (node) are sufficient.

Example:
Select all FEs connected together (Starting from node 13):
```julia
l = selectelem(fens, fes, flood = true, startnode = 13)
```

### Optional keyword arguments
Should we consider the element only if all its nodes are in?
- `allin` = Boolean: if true, then all nodes of an element must satisfy
the criterion; otherwise  one is enough.

# Output
`felist` = list of finite elements from the set that satisfy the criteria
"""
function selectelem(fens::FENodeSet, fes::T; kwargs...) where {T<:AbstractFESet}

    # smoothpatch
    #
    # Select all FEs that are part of a smooth surface.
    # For instance, starting from the finite element number 13, select all finite elements
    # whose normals defer from the normal of the Neighbor element by less than
    # 0.05 in the sense that dot(n1,n2)>sqrt(1-0.05^2).
    #     fe_select(fens,fes,struct ('flood', true, 'startfe', 13,  'normaldelta',0.05))
    # # Select all FEs "facing" in the direction [x(1),x(2),0] (away from the z-axis):
    #     fe_select(fens,fes,struct ('facing', true, 'direction', @(x)(+[x(1:2),0])))
    # Here x is the centroid of the nodes of each selected fe.
    # # Select all FEs "facing" in the direction [x(1),x(2),0] (away from the z-axis):
    #     fe_select(fens,fes,struct ('facing',1,'direction',@(x)(+[x(1:2),0]),'tolerance',0.01))
    # Here the fe is considered facing in the given direction if the dot
    # product of its normal and the direction vector is greater than tolerance.

    # Extract arguments
    allin = nothing; flood = nothing; facing = nothing; label = nothing; withnodes = nothing
    inflate = 0.0
    # nearestto = nothing; smoothpatch = nothing;
    startnode = 0; dotmin = 0.01
    overlappingbox = nothing
    for apair in pairs(kwargs)
        sy, val = apair
        if sy == :flood
            flood = val
        elseif sy == :facing
            facing = val
        elseif sy == :label
            label = val
        elseif sy == :overlappingbox
            overlappingbox = val
        elseif sy == :nearestto
            nearestto = val
        elseif sy == :inflate
            inflate = val
        elseif sy == :withnodes
            withnodes = val
        elseif sy == :allin
            allin = val
        end
    end

    if flood !== nothing
        for apair in pairs(kwargs)
            sy, val = apair
            if sy == :startnode
                startnode = val
            end
        end
    end

    if facing != nothing
        facing = true;
        direction = nothing
        dotmin = 0.01;
        for apair in pairs(kwargs)
            sy, val = apair
            if sy == :direction
                direction = val
            elseif (sy == :dotmin) || (sy == :tolerance)# allow for obsolete keyword  to work
                dotmin = val
            end
        end
    end

    # if isfield(options,'smoothpatch')
    #     smoothpatch = true;
    #     if isfield(options,'normaldelta')
    #         normaldelta = options.normaldelta;
    #     else
    #         normaldelta= 0.05;
    #     end
    #     if isfield(options,'startfe')
    #         startfe = options.startfe;
    #     else
    #         error('Need the number of the starting finite element');
    #     end
    # end

    # if isfield(options, 'nearestto')
    #     nearestto = true;
    #     locations = options.nearestto;
    # end

    # The  elements of this array are flipped from zero  when the element satisfies
    # the search condition.. This list is  eventually purged of the zero elements and
    # returned.
    felist = zeros(FInt,length(fes.conn));

    #     Select based on fe label
    if label!= nothing
        @assert length(fes.label) == length(fes.conn)
        for i=1:length(fes.conn)
            if label==fes.label[i]
                felist[i] =i;   # matched this element
            end
        end
        return  felist[findall(x->x!=0, felist)]; # return the nonzero element numbers
    end

    # Select by flooding
    if flood != nothing && (flood)
        @assert startnode > 0
        fen2fe = FENodeToFEMap(connasarray(fes), count(fens))
        felist = zeros(FInt, count(fes));
        pfelist = zeros(FInt, count(fes));
        felist[fen2fe.map[startnode]] .= 1;
        while true
            copyto!(pfelist, felist);
            markedl = findall(x -> x != 0, felist)
            for j = markedl
                for k = fes.conn[j]
                    felist[fen2fe.map[k]] .= 1;
                end
            end
            if sum(pfelist-felist) == 0 # If there are no more changes in this pass, we are done
                break;
            end
        end
        return findall(x -> x != 0, felist); # return the nonzero element numbers;
    end

    # Helper Function: are all nodes of an element in the list?
    function nodesin(c, withnodes)
        for i in c
            if !(i in withnodes)
                return false
            end
        end
        return true
    end

    # Select all elements whose nodes are in a given list
    if withnodes != nothing
        for i=1:length(fes.conn)
            if nodesin(fes.conn[i], withnodes)
                felist[i] =i;   # matched this element
            end
        end
        return  felist[findall(x->x!=0, felist)]; # return the nonzero element numbers
    end

    # Helper function: calculate the normal to a boundary finite element
    function normal(Tangents)
        sdim, mdim = size(Tangents);
        @assert (mdim==1) || (mdim==2)
        if     mdim==1 # 1-D fe
            N = [Tangents[2,1],-Tangents[1,1]];
            return N/norm(N);
        elseif mdim==2 # 2-D fe
            N= cross(Tangents[:,1],Tangents[:,2])
            return N/norm(N);
        end
    end

    # Select by in which direction the normal of the fes face
    if (facing != nothing) && (facing)
        xs =fens.xyz;
        sd =spacedim(fens);
        md =manifdim(fes);
        @assert (md == sd-1) "'Facing': only for Manifold dim. == Space dim.-1"
        param_coords =zeros(FFlt,1,md);
        Need_Evaluation = (typeof(direction) <: Function);
        if ( !Need_Evaluation)
            d = reshape(direction,1,sd)/norm(direction);
        end
        Nder = bfundpar(fes, vec(param_coords));
        for i=1:length(fes.conn)
            xyz =xs[[j for j in fes.conn[i]],:];
            Tangents =xyz'*Nder;
            N = normal(Tangents);
            if (Need_Evaluation)
                d = direction(mean(xyz, dims = 1));
                d = reshape(d,1,sd)/norm(d);
            end
            if (dot(vec(N),vec(d)) > dotmin)
                felist[i]=i;
            end
        end
        return  felist[findall(x->x!=0, felist)]; # return the nonzero element numbers
    end

    #   # Select by the change in normal
    #   if (smoothpatch)
    #     mincos=sqrt(1-normaldelta^2);
    #     xs=fens.xyz;
    #     sdim =size(xs,2);
    #     mdim=fes.dim;
    #     if (mdim~=sdim-1)
    #       error ('"Smoothpatch" selection of fes make sense only for Manifold dimension ==Space dimension-1')
    #     end
    #     fen2fe_map=fenode_to_fe_map (struct ('fes',fes));
    #     gmap=fen2fe_map.map;
    #     conn=fes.conn;
    #     param_coords =zeros(1,mdim);# This is a hack: this may not be the proper location for all element types
    #     Nder = bfundpar (fes, param_coords);
    #     felist= zeros(1, size(conn,1));
    #     felist(startfe)=1;
    #     while true
    #       pfelist=felist;
    #       markedl=find(felist~=0);
    #       for j=markedl
    #         xyz =xs(conn(j,:),:);
    #         Tangents =xyz'*Nder;
    #         Nj = normal (Tangents,sdim, mdim);
    #         for k=conn(j,:)
    #           for ke=gmap{k}
    #             xyz =xs(conn(ke,:),:);
    #             Tangents =xyz'*Nder;
    #             Nk = normal (Tangents,sdim, mdim);
    #             if (dot(Nj,Nk)>mincos)
    #               felist(ke)=1;
    #             end
    #           end
    #         end
    #       end
    #       if sum(pfelist-felist)==0, break; end
    #     end
    #     felist =find(felist~=0);
    #     return;
    #   end

    # Select all FEs whose bounding box overlaps given box
    if (overlappingbox !== nothing)
        obox = deepcopy(overlappingbox) # do not mess with the array on input
        inflatebox!(obox, inflate)
        bbox = zeros(2 * size(fens.xyz, 2))
        for i in 1:length(fes.conn)
            bbox = initbox!(bbox, vec(fens.xyz[fes.conn[i][1], :]))
            for k in 2:length(fes.conn[i])
                bbox = updatebox!(bbox, fens.xyz[fes.conn[i][k], :])
            end
            if boxesoverlap(obox, bbox)
                felist[i]=i; # this element overlaps the box
            end
        end
        return  felist[findall(x->x!=0, felist)]; # return the nonzero element numbers
    end

    #   # get the fe nearest to the supplied point
    #   if (nearestto)

    #     felist = [];

    #     # get the locations of all the nodes
    #     xs = fens.xyz;
    #     #     Disqualify all the nodes that are not connected to the finite
    #     #     elements on input
    #     cn = connected_nodes(fes);
    #     xs(setdiff(1:size(xs,1),cn),:) =Inf;

    #     # get the connectivity of all the FEs
    #     conns = fes.conn;

    #     for i = 1:size(locations,1)

    #       # Get the smallest distances between the node locations and the
    #       # desired location using ipdm (by John D'Errico)
    #       distance = ipdm(xs, locations(i,:));

    #       # Get array index of smallest distance to location from
    #       # all nodes
    #       [junk,IX] = min(distance);

    #       # Find the fes connected to the nearest node
    #       [crows, ccols] = find(conns == IX(1));

    #       if numel(crows) == 1
    #         # if there's only one, this is superb, we can add it to the
    #         # list and continue to the next location
    #         felist(end+1) = crows;
    #       else

    #         # otherwise we must determine which cell is closest based on
    #         # the other connected nodes
    #         nearconns = conns(crows, :);

    #         nearconndist = zeros(size(nearconns));

    #         # Get the distances between all the nodes in the nearest
    #         # FEs and the nearest node to the location
    #         for j = 1:size(nearconns, 1)
    #           for k = 1:size(nearconns, 2)
    #             nearconndist(j,k) = ipdm(xs(nearconns(j,k), :), locations(i,:));
    #           end
    #         end

    #         # set the distance from the nearest node to in the cell to
    #         # the location to infinity
    #         nearconndist(nearconns == IX(1)) = inf;

    #         # order the FEs by the smallest distance of any connected
    #         # node from the nearest node
    #         [junk,IX] = sort(min(nearconndist,[],2));

    #         # return the index in conns to the first of these ordered
    #         # FEs
    #         if (~isempty(crows))
    #           felist(end+1) = crows(IX(1));
    #         end

    #       end

    #     end

    #     return;

    #   end

    #  Default:   Select based on location of nodes
    #   Should we consider the element only if all its nodes are in?
    allinvalue = (allin == nothing) || ((allin != nothing) && (allin))
    # Select elements whose nodes are in the selected node list
    nodelist = selectnode(fens; kwargs...);
    nper = nodesperelem(fes)
    for i = 1:length(fes.conn)
        found = 0
        for nd in fes.conn[i]
            if nd in nodelist
                found = found + 1
            end
        end
        if allinvalue
            if found == nper
                felist[i] =i;
            end
        else
            if found >= 1
                felist[i] =i;
            end
        end
    end
    return felist[findall(x->x!=0, felist)]; # return the nonzero element numbers

end

function _compute_vlist!(vlist::Vector{FInt}, abox::T, sdim::FInt, v::FFltMat) where {T}
    # Helper functions
    @inline inrange(rangelo,rangehi,x) = (rangelo <= x <= rangehi)
    nn = 0
    for j in 1:size(v,1)
        matches = true
        for i in 1:sdim
            if !inrange(abox[2*i-1], abox[2*i], v[j, i])
                matches = false; break
            end
        end
        if matches
            nn = nn + 1; vlist[nn] = j;
        end
    end
    return vlist, nn
end

"""
    vselect(v::FFltMat; kwargs...)

Select locations (vertices) from the array based on some criterion.

See the function `selectnode()` for examples of the criteria that can be
used to search vertices.
"""
function vselect(v::FFltMat; kwargs...)
    # Extract arguments
    box = nothing; distance = nothing; from = nothing; plane  =  nothing;
    thickness = nothing; nearestto = nothing; inflate = 0.0;
    for apair in pairs(kwargs)
        sy, val = apair
        if sy == :box
            box = val
        elseif sy == :distance
            distance = val
        elseif sy == :from
            from = val
        elseif sy == :plane
            plane = val
        elseif sy == :thickness
            thickness = val
        elseif sy == :nearestto
            nearestto = val
        elseif sy == :inflate
            inflate = val
        end
    end

    # Did we get an inflate value
    inflatevalue = 0.0;
    if inflate != nothing
        inflatevalue = FFlt(inflate);
    end

    # Initialize the output list
    vlist = zeros(FInt, size(v,1)); nn = 0;

    # Process the different options
    if box != nothing
        sdim = size(v,2)
        dim = Int(round(length(box)/2.))::FInt;
        @assert dim == sdim "Dimension of box not matched to dimension of array of vertices"
        abox = vec(box)::FFltVec
        inflatebox!(abox, inflatevalue)
        vlist, nn = _compute_vlist!(vlist, abox, sdim, v)
    elseif distance != nothing
        fromvalue =fill!(deepcopy(v[1,:]), 0.0);
        if from!=nothing
            fromvalue = from;
        end
        fromvalue = reshape(fromvalue, (size(v[1,:])))
        d=distance+inflatevalue;
        for i=1:size(v,1)
            if norm(fromvalue-v[i,:])<d
                nn =nn +1; vlist[nn] =i;
            end
        end
    elseif plane != nothing
        normal = plane[1:end-1];
        normal = vec(normal/norm(normal));
        thicknessvalue = 0.0;
        if thickness != nothing
            thicknessvalue = thickness;
        end
        t = thicknessvalue+inflatevalue;
        distance = plane[end];
        for i = 1:size(v,1)
            ad = dot(v[i,:], normal);
            if abs(distance-ad)<t
                nn = nn +1; vlist[nn] =i;
            end
        end
    elseif nearestto != nothing
        location = vec(nearestto);
        distance =  zeros(size(v,1));
        for i=1:size(v,1)
            distance[i] = norm(location-vec(v[i,:]))
        end
        Mv,j = findmin(distance)
        vlist[1] = j;
        nn=1;
    end
    if (nn==0)
        vlist = FInt[];# nothing matched
    else
        vlist = vlist[1:nn];
    end
    return vlist
end

"""
    findunconnnodes(fens::FENodeSet, fes::AbstractFESet)

Find nodes that are not connected to any finite element.

connected = array is returned which is for the node k either true (node k is
     connected), or false (node k is not connected).

Let us say there are nodes not connected to any finite element that you
would like to remove from the mesh: here is how that would be
accomplished.
"""
function findunconnnodes(fens::FENodeSet, fes::AbstractFESet)
    connected = trues(count(fens));
    fen2fem = FENodeToFEMap(connasarray(fes), count(fens))
    for i=1:length(fen2fem.map),
        connected[i] = (!isempty(fen2fem.map[i]));
    end
    return connected
end

end

# NOTE: This operation is probably best done with the  node-to-element  map.
# """
#     connectedelems(fes::AbstractFESet, node_list::FIntVec)
#
# Extract the numbers of the finite elements connected to given nodes.
#
# Extract the list of numbers for the fes  that are connected to given nodes.
#
# Warning: this tends to be an expensive operation.
# """
# function connectedelems(fes::AbstractFESet, node_list::FIntVec)
#   cg=zeros(FInt,size(fes.conn,1));
#   for j=1:size(fes.conn,1)
#     cg[j]= length( intersect(fes.conn[j,:], node_list) );
#   end
#   cg =find_n(cg);
# end


text:
# A test for abort in the generic callback.
# Solve:
    # minimize x
    # s.t. x is integral: 0 <= x <= 1

@testset "GenCB_01" begin
    env = CPLEX.Env()
    CPLEX.set_param!(env, "CPX_PARAM_THREADS", 1)
    model = CPLEX.Model(env, "callback_test")
    CPLEX.add_vars!(model, [1.0], [0.0], [1.0])
    CPLEX.c_api_chgctype(model, Int32[1], Cchar['I'])
    function my_callback(cb_context::CPLEX.CallbackContext, context_id::Clong)
        CPLEX.cbabort(cb_context)
    end
    CPLEX.cbsetfunc(model, CPLEX.CPX_CALLBACKCONTEXT_CANDIDATE, my_callback)
    CPLEX.optimize!(model)
    @test CPLEX.get_status(model) == :CPXMIP_ABORT_INFEAS
end


text:
using Mocking

describe_active_receipt_rule_set_1_patch = @patch function describe_active_receipt_rule_set(a...; b...)
    return Dict{String, Any}("Rules" => Any[Dict{String, Any}("Actions" => Any[Dict{String, Any}("S3Action" => Dict{String, Any}("ObjectKeyPrefix" => "email", "BucketName" => "MyBucket"))], "TlsPolicy" => "Optional", "Enabled" => true, "Name" => "MyRule", "ScanEnabled" => true)], "Metadata" => Dict{String, Any}("Name" => "default-rule-set", "CreatedTimestamp" => "2016-07-15T16:25:59.607Z"))
end

describe_receipt_rule_1_patch = @patch function describe_receipt_rule(a...; b...)
    return Dict{String, Any}("Rule" => Dict{String, Any}("Actions" => Any[Dict{String, Any}("S3Action" => Dict{String, Any}("ObjectKeyPrefix" => "email", "BucketName" => "MyBucket"))], "TlsPolicy" => "Optional", "Enabled" => true, "Name" => "MyRule", "ScanEnabled" => true))
end

describe_receipt_rule_set_1_patch = @patch function describe_receipt_rule_set(a...; b...)
    return Dict{String, Any}("Rules" => Any[Dict{String, Any}("Actions" => Any[Dict{String, Any}("S3Action" => Dict{String, Any}("ObjectKeyPrefix" => "email", "BucketName" => "MyBucket"))], "TlsPolicy" => "Optional", "Enabled" => true, "Name" => "MyRule", "ScanEnabled" => true)], "Metadata" => Dict{String, Any}("Name" => "MyRuleSet", "CreatedTimestamp" => "2016-07-15T16:25:59.607Z"))
end

get_account_sending_enabled_1_patch = @patch function get_account_sending_enabled(a...; b...)
    return Dict{String, Any}("Enabled" => true)
end

get_identity_dkim_attributes_1_patch = @patch function get_identity_dkim_attributes(a...; b...)
    return Dict{String, Any}("DkimAttributes" => Dict{String, Any}("user@example.com" => Dict{String, Any}("DkimVerificationStatus" => "NotStarted", "DkimEnabled" => false), "example.com" => Dict{String, Any}("DkimVerificationStatus" => "Success", "DkimTokens" => Any["EXAMPLEjcs5xoyqytjsotsijas7236gr", "EXAMPLEjr76cvoc6mysspnioorxsn6ep", "EXAMPLEkbmkqkhlm2lyz77ppkulerm4k"], "DkimEnabled" => true)))
end

get_identity_mail_from_domain_attributes_1_patch = @patch function get_identity_mail_from_domain_attributes(a...; b...)
    return Dict{String, Any}("MailFromDomainAttributes" => Dict{String, Any}("example.com" => Dict{String, Any}("BehaviorOnMXFailure" => "UseDefaultValue", "MailFromDomain" => "bounces.example.com", "MailFromDomainStatus" => "Success")))
end

get_identity_notification_attributes_1_patch = @patch function get_identity_notification_attributes(a...; b...)
    return Dict{String, Any}("NotificationAttributes" => Dict{String, Any}("example.com" => Dict{String, Any}("BounceTopic" => "arn:aws:sns:us-east-1:EXAMPLE65304:ExampleTopic", "ForwardingEnabled" => true, "HeadersInBounceNotificationsEnabled" => false, "HeadersInComplaintNotificationsEnabled" => false, "HeadersInDeliveryNotificationsEnabled" => false, "ComplaintTopic" => "arn:aws:sns:us-east-1:EXAMPLE65304:ExampleTopic", "DeliveryTopic" => "arn:aws:sns:us-east-1:EXAMPLE65304:ExampleTopic")))
end

get_identity_policies_1_patch = @patch function get_identity_policies(a...; b...)
    return Dict{String, Any}("Policies" => Dict{String, Any}("MyPolicy" => "{\"Version\":\"2008-10-17\",\"Statement\":[{\"Sid\":\"stmt1469123904194\",\"Effect\":\"Allow\",\"Principal\":{\"AWS\":\"arn:aws:iam::123456789012:root\"},\"Action\":[\"ses:SendEmail\",\"ses:SendRawEmail\"],\"Resource\":\"arn:aws:ses:us-east-1:EXAMPLE65304:identity/example.com\"}]}"))
end

get_identity_verification_attributes_1_patch = @patch function get_identity_verification_attributes(a...; b...)
    return Dict{String, Any}("VerificationAttributes" => Dict{String, Any}("example.com" => Dict{String, Any}("VerificationToken" => "EXAMPLE3VYb9EDI2nTOQRi/Tf6MI/6bD6THIGiP1MVY=", "VerificationStatus" => "Success")))
end

get_send_quota_1_patch = @patch function get_send_quota(a...; b...)
    return Dict{String, Any}("SentLast24Hours" => 1, "Max24HourSend" => 200, "MaxSendRate" => 1)
end

get_send_statistics_1_patch = @patch function get_send_statistics(a...; b...)
    return Dict{String, Any}("SendDataPoints" => Any[Dict{String, Any}("DeliveryAttempts" => 5, "Complaints" => 0, "Rejects" => 0, "Bounces" => 0, "Timestamp" => "2016-07-13T22:43:00Z"), Dict{String, Any}("DeliveryAttempts" => 3, "Complaints" => 0, "Rejects" => 0, "Bounces" => 0, "Timestamp" => "2016-07-13T23:13:00Z"), Dict{String, Any}("DeliveryAttempts" => 1, "Complaints" => 0, "Rejects" => 0, "Bounces" => 0, "Timestamp" => "2016-07-13T21:13:00Z")])
end

list_identities_1_patch = @patch function list_identities(a...; b...)
    return Dict{String, Any}("NextToken" => "", "Identities" => Any["user@example.com"])
end

list_identity_policies_1_patch = @patch function list_identity_policies(a...; b...)
    return Dict{String, Any}("PolicyNames" => Any["MyPolicy"])
end

list_receipt_filters_1_patch = @patch function list_receipt_filters(a...; b...)
    return Dict{String, Any}("Filters" => Any[Dict{String, Any}("Name" => "MyFilter", "IpFilter" => Dict{String, Any}("Cidr" => "1.2.3.4/24", "Policy" => "Block"))])
end

list_receipt_rule_sets_1_patch = @patch function list_receipt_rule_sets(a...; b...)
    return Dict{String, Any}("NextToken" => "", "RuleSets" => Any[Dict{String, Any}("Name" => "MyRuleSet", "CreatedTimestamp" => "2016-07-15T16:25:59.607Z")])
end

list_verified_email_addresses_1_patch = @patch function list_verified_email_addresses(a...; b...)
    return Dict{String, Any}("VerifiedEmailAddresses" => Any["user1@example.com", "user2@example.com"])
end

send_email_1_patch = @patch function send_email(a...; b...)
    return Dict{String, Any}("MessageId" => "EXAMPLE78603177f-7a5433e7-8edb-42ae-af10-f0181f34d6ee-000000")
end

send_raw_email_1_patch = @patch function send_raw_email(a...; b...)
    return Dict{String, Any}("MessageId" => "EXAMPLEf3f73d99b-c63fb06f-d263-41f8-a0fb-d0dc67d56c07-000000")
end

verify_domain_dkim_1_patch = @patch function verify_domain_dkim(a...; b...)
    return Dict{String, Any}("DkimTokens" => Any["EXAMPLEq76owjnks3lnluwg65scbemvw", "EXAMPLEi3dnsj67hstzaj673klariwx2", "EXAMPLEwfbtcukvimehexktmdtaz6naj"])
end

verify_domain_identity_1_patch = @patch function verify_domain_identity(a...; b...)
    return Dict{String, Any}("VerificationToken" => "eoEmxw+YaYhb3h3iVJHuXMJXqeu1q1/wwmvjuEXAMPLE")
end


text:
#########
# Width
getmid(img) = fld1(size(img,1),2)
getmid(exi::ExpImgs) = fld1(size(exi)[1],2)

top2ind(top, img_or_exi) = getmid(img_or_exi) - top + 1
bottom2ind(bottom, img_or_exi) = bottom + getmid(img_or_exi)
top_inds(img_or_exi) = getmid(img_or_exi):-1:1
bottom_inds(img_or_exi) = getmid(img_or_exi)+1:size(img,1)

################
# width via color tresholding

function findtrough(img, title; verbose=false)
    nbins = 100
    hi = fit(Histogram, img[:], closed=:left, nbins=nbins)
    peaks, ips, proms, doms = findpeaks(-hi.weights[1:(end-end÷4)])
    i = findmax(proms[ips.<(nbins÷2)])[2] # look only if first half of histogram
    if proms[i]<length(img)/500
        warn("No color boundary found in $title")
        return 0
    end
    verbose && P.plot(1:length(hi.weights), -hi.weights)
    #verbose && P.bar(hi.edges[1][1:end-1], hi.weights)
    verbose && @show ips[i]
    ips[i]
end

"""

"""
function keep_only_big_regions(img, loc, siz=round(Int,(sqrt(length(img))/20)^2),
                               k=10; verbose=false)
    segments = felzenszwalb(img, k, siz)
    verbose && @show segments

    # map(i->segment_mean(segments,i), labels_map(segments))
    tmp = labels_map(segments)
    # set region at `loc` to true, rest to false.
    tmp.==tmp[loc...]
end

function channel_width_thresholding(img_, loc, median_filter_region, title; verbose=false)
    trough = findtrough(img_, title, verbose=verbose)
    trough==0 && return zeros(Int,size(img_,2)),zeros(Int,size(img_,2)),zeros(Int,size(img_,2))
    img_ = img_.>trough
    if median_filter_region!=(0,0)
        img_ = median_filter(img_, median_filter_region)
    end
    img_3 = keep_only_big_regions(img_, loc,
                                  round(Int,(sqrt(length(img_))/50)^2),
                                  verbose=verbose);
    verbose && imshow(img_);
    verbose && imshow(img_3);

    mid = getmid(img_)
    # top = squeeze(sum(img_3[1:mid,:],1),1)
    # bottom = squeeze(sum(img_3[mid+1:end,:],1),1)
    top = [findlast(img_3[mid:-1:1,i]) for i=1:size(img_,2)]
    bottom = [findlast(img_3[mid+1:end,i]) for i=1:size(img_,2)]

    # remove outliers
    quantile_filter!(top, 0.05, 5, 1, 5, 6)
    quantile_filter!(bottom, 0.05, 5, 1, 5, 6)

    return top, bottom, top+bottom
end

###############
# Edge detection
"""
    channel_width_edgedetection(img_, minhalfwidth, gauss_w = 10, quant=0.8, gap=11; verbose=false)

Determine width using edge detection.
"""
channel_width_edgedetection(img_, exi; verbose=false) =
    channel_width_edgedetection(img_, exi.minhalfwidth, exi.gauss_w, exi.quant, exi.gap, verbose=verbose)
function channel_width_edgedetection(img_, minhalfwidth, gauss_w,
                                     quant, gap, morphgrad=false; verbose=false)
    @assert !iseven(gap) "gap must be a odd number"
    # works on img_
    verbose && imshow(img_);

    mid = getmid(img_)
    ncols = size(img_,2)

    # Gauss filter:
    img_f = imfilter(img_, Kernel.gaussian(gauss_w));

    # eges:
    if morphgrad
        mag = morphogradient(img_f)
    else
        grad_x, grad_y, mag, orient = imedge(img_f);
        # Only look at gradient in x-direction
        # and only when going from more green to less green
        # in outward direction.
        grad_x[mid:end,:] = -grad_x[mid:end,:]
        grad_x = -grad_x
        grad_x[grad_x.<0] = 0
        mag = grad_x
    end
    verbose && imshow(mag);
    #verbose && imshow(grad_x);

    # # Find significant peaks on each column:
    # out = falses(mag);
    # out_ = zeros(mag);
    # for i=1:ncols
    #     peaks,locs,proms,doms = findpeaks(mag[:,i]);
    #     big = find(x->x>quantile(peaks, quant) && x>quantile(proms, quant), peaks)
    #     out[locs[big],i] = true
    #     out_[locs,i] = proms
    # end
    # verbose && imshow(out);

    # Find significant peaks on each column, version 2:
    out = falses(img_);
    vverbose = false;
    for i=1:ncols
        # right (lower) border
        # peaks,locs,proms,doms = findpeaks(mag[:,i]);
        peaks_t,locs_t,proms_t,doms_t = findpeaks(mag[mid:-1:1,i]);
        peaks_b,locs_b,proms_b,doms_b = findpeaks(mag[mid+1:end,i]);
        locs_tt = mid +1 - locs_t
        locs_bb = locs_b-1+mid

        big_t = find_spikes_quantile(proms_t, >,
                                     0.9, mean(proms_t), 1,
                                     length(proms_t)÷10, 0,
                                     verbose=vverbose)[1]
        big_b = find_spikes_quantile(proms_b, >,
                                     0.9, mean(proms_b), 1,
                                     length(proms_b)÷10, 0,
                                     verbose=vverbose)[1]

        out[locs_tt[big_t],i] = true
        out[locs_bb[big_b],i] = true
    end
    verbose && imshow(out);

    # # find most prominent peak on each half:
    # out_p = falses(mag);
    # for i=1:ncols
    #     # right (lower) border
    #     peaks,locs,proms,doms = findpeaks(mag[:,i]);
    #     big = find(x->x>quantile(peaks, quant) && x>quantile(proms, quant), peaks)
    #     out[locs[big],i] = true
    #     out_[locs,i] = proms
    # end
    # verbose && imshow(out);

    # # Just find the two biggest peaks in both halves
    # top = zeros(Int, ncols);
    # bottom = zeros(Int, ncols);
    # for i=1:ncols
    #     top[i] = findmax(mag[1:mid-minhalfwidth,i])[2]
    #     bottom[i] = findmax(mag[mid+minhalfwidth+1:end,i])[2] + minhalfwidth
    # end
    # return top, bottom, top+bottom

    ## remove small bits
    # labs0 = label_components(out, trues(3,3));
    # inds0 = component_indices(labs0);
    # len0 = component_lengths(labs0);
    # # drop single & double pixels
    # for (i,l) in enumerate(len0)
    #     if l<3
    #         labs0[inds0[i]] = 0
    #     end
    # end
    # out = labs0.>0;
    # verbose && imshow(out);

    ## join bigger ones and remove left-over smaller ones
    labs = label_components(out, trues(gap,gap));
    inds = component_indices(labs);
    len = component_lengths(labs);
    # drop all components with are below certain length
    tlen = 50
    for (i,l) in enumerate(len)
        if l<tlen
            labs[inds[i]] = 0
        end
    end
    out = labs.>0;
    verbose && imshow(labs.>0);

    ## TODO: actually join the big pieces

    # now determine distance from middle+minhalfwidth to first peak
    top = zeros(Int, size(mag,2))
    bottom = zeros(Int, size(mag,2))
    for i=1:size(mag,2)
        tmp = findfirst(@view out[mid-minhalfwidth:-1:1,i])
        top[i] = tmp==0 ? -1 : tmp + minhalfwidth
        tmp = findfirst(@view out[mid+minhalfwidth+1:end, i])
        bottom[i] = tmp==0 ? -1 : tmp + minhalfwidth
    end
    verbose && PyPlot.figure()
    verbose && PyPlot.plot(top, ":")
    verbose && PyPlot.plot(bottom, ":")

    # remove outliers
    quantile_filter!(top,
                     0.05, 5, 1, 5, 6)
    quantile_filter!(bottom,
                     0.05, 5, 1, 5, 6)

    # top[top.==-1] = NaN
    # bottom[bottom.==-1] = NaN

    verbose && PyPlot.plot(top)
    verbose && PyPlot.plot(bottom)
    return top, bottom, top+bottom
end

"""
    filter1d_outliers!(vec::Vector{<:Integer}, halfwindow=1)

Filter outliers using median.  Probably better ways to do this.
"""
function filter1d_outliers!(vec::Vector{<:Integer}, halfwindow=1)
    len = length(vec)
    for i=eachindex(vec)
        vec[i] = round.(eltype(vec), median(@view vec[max(1,i-halfwindow):min(len,i+halfwindow)]))
    end
    vec
end

# function plottop()
#     for i=1:size(mag,2)
#         PyPlot.plot(top[i]*0+i, top[i],".")
#         PyPlot.plot(bottom[i]*0+i, bottom[i],".")
#     end
# end

###############
#

"""
    channel_width(exi::ExpImgs; verbose=false)

Returns the best of both outlines for all images in a exi.  Returns

- tops, bottoms: pixel distance to top and bottom outline in terms of
  pixels of the cropped image (i.e. not corrected with center_dist)
- t_dist, b_dist: distance in meters to top and bottom from a line
  through the original image center and parallel to (exi.p1, exi.p2).

"""

function channel_width(exi::ExpImgs; verbose=false, vverbose=false, saveit=true,
                       overwrite=false, store_imgs=false)
    @unpack dir, ns, p1, p2, thin_num, algo,
            minhalfwidth_orig, gauss_w, quant, gap, median_filter_region = exi

    println("\n Figuring out R-channel width for pictures in dir: $dir\n")
    minhalfwidth = minhalfwidth_orig÷thin_num

    tops = Int[]
    bottoms = Int[]
    last_top = [minhalfwidth]
    last_bottom = [minhalfwidth]
    ncol = size(exi)[2]
    center_dists = Int[]
    thumbs = []
    imgs = []
    time = DateTime[]
    @showprogress for (n,fl) in zip(ns, image_files(exi))
        img, center_dist, capture_time = prep_img(fl, exi; verbose=vverbose)
        t, b = _channel_width(img, last_top, last_bottom, exi, "$n:  $fl",
                             algo,
                             verbose=verbose,
                             vverbose=vverbose)
        last_top = t
        last_bottom = b
        append!(tops, t)
        append!(bottoms,b)
        push!(center_dists, center_dist)
        push!(thumbs, thin(img, 16))
        store_imgs && push!(imgs, img)
        push!(time, capture_time)
    end
    tops = reshape(tops, size(exi)[2], length(tops)÷size(exi)[2])
    bottoms = reshape(bottoms, size(exi)[2], length(tops)÷size(exi)[2])
    ts_dist = pixel2meter.(tops.+center_dists', exi)
    bs_dist = pixel2meter.(bottoms.-center_dists', exi)
    error()
    # add support for
    res = ExpImgsResults(exi,
                         time,
                         t,
                         trun,
                         irun,
                         center_dists,
                         tops,
                         bottoms,
                         ts_dist,
                         bs_dist,
                         rchannel_stats(ts_dist,bs_dist)...,
                         thumbs,
                         imgs,
                         Dict()
                         )
    saveit && save_result(res, overwrite=overwrite, store_imgs=store_imgs)
    return res
end

function _channel_width(img, last_top, last_bottom, exi::ExpImgs, title, algo;
                       verbose=false, vverbose=false)
    @unpack minhalfwidth_orig, gauss_w, quant, gap, median_filter_region, thin_num = exi
    minhalfwidth = minhalfwidth_orig÷thin_num

    if length(last_bottom)==1
        last_bottom = zeros(Int, size(img,2)) .+ last_bottom
    end
    if length(last_top)==1
        last_top = zeros(Int, size(img,2)) .+ last_top
    end

    # minhalfwidth_cur = max(min(minimum(last_top), minimum(last_bottom)) - minhalfwidth÷5, minhalfwidth)
    minhalfwidth_cur = minhalfwidth
    if algo in [:both, :thresh]
        top_t, bottom_t, tot_t = channel_width_thresholding(img,
                                                            (size(img,1)÷2, size(img,2)÷2),
                                                            median_filter_region, title;
                                                            verbose=vverbose)
        if algo==:thresh
            new_top, new_bottom = top_t, bottom_t
        end
    end
    if algo in [:both, :edge]
        top_e, bottom_e, tot_e = channel_width_edgedetection(img, minhalfwidth_cur,
                                                             gauss_w, quant, gap;
                                                             verbose=vverbose)
        if algo==:edge
            new_top, new_bottom = top_e, bottom_e
        end
    end
    if !(algo in [:both, :edge, :thresh])
        error("algo=$algo not recognized")
    end

    if algo==:both
        # now produce a best vector somehow:

        # Go through to find points of agreement which are farther than v_old.
        # Of those pick the closest to v_old, this will serve as starting point
        tv, bv = -1, -1
        ta, ba = 10^10, 10^10
        ti, bi = 0, 0
        for i=1:length(top_e)
            tmp_,vv = how_close(top_t[i], top_e[i], last_top[i])
            if tmp_<ta
                ta = tmp_
                ti = i
                tv = vv
            end
            tmp_,vv = how_close(bottom_t[i], bottom_e[i], last_bottom[i])
            if tmp_<ba
                ba = tmp_
                bi = i
                bv = vv
            end
        end
        if tv<1
            # no starting point found!
            error("no starting point found for tops in: $title")
        end
        if bv<1
            error("no starting point found for bottoms in: $title")
        end

        new_top = zeros(last_top)
        new_top[ti] = tv
        if ti<length(new_top)
            new_top[ti+1] = pick_next_point(top_t[ti+1], top_e[ti+1], new_top[ti], last_top[ti])
            for i=ti+2:length(new_top)
                new_top[i] = pick_next_point(top_t[i], top_e[i], new_top[i-1], new_top[i-2], last_top[i])
            end
        end
        if ti>1
            new_top[ti-1] = pick_next_point(top_t[ti-1], top_e[ti-1], new_top[ti], last_top[ti])
            for i=ti-2:-1:1
                new_top[i] = pick_next_point(top_t[i], top_e[i], new_top[i+1], new_top[i+2], last_top[i])
            end
        end

        new_bottom = zeros(last_bottom)
        new_bottom[bi] = bv
        if bi<length(new_bottom)
            # start
            new_bottom[bi+1] = pick_next_point(bottom_t[bi+1], bottom_e[bi+1], new_bottom[bi], last_bottom[bi])
            for i=bi+2:length(new_bottom)
                new_bottom[i] = pick_next_point(bottom_t[i], bottom_e[i], new_bottom[i-1], new_bottom[i-2], last_bottom[i])
            end
        end
        if bi>1
            # start
            new_bottom[bi-1] = pick_next_point(bottom_t[bi-1], bottom_e[bi-1], new_bottom[bi], last_bottom[bi])
            for i=bi-2:-1:1
                new_bottom[i] = pick_next_point(bottom_t[i], bottom_e[i], new_bottom[i+1], new_bottom[i+2], last_bottom[i])
            end
        end
    end

    window = 10
    niter = 30
    quantile_filter!(new_top, 0.1, 2, 1, window, window+1, niter)
    quantile_filter!(new_bottom, 0.1, 2, 1, window, window+1, niter)

    window = 11
    new_top = round.(Int,mapwindow(mean, new_top, window))
    new_bottom = round.(Int,mapwindow(mean, new_bottom, window))

    if verbose
        if algo==:both
            plot_all_n_new(img, new_top, new_bottom, top_t, bottom_t, top_e, bottom_e;
                           col="r", label="", title=title,
                           legend=false)
        elseif algo==:edge
            plot_all_n_new(img, new_top, new_bottom, top_e, bottom_e, top_e, bottom_e;
                           col="r", label="", title=title,
                           legend=false)
        elseif algo==:thresh
            plot_all_n_new(img, new_top, new_bottom, top_t, bottom_t, top_t, bottom_t;
                           col="r", label="", title=title,
                           legend=false)
        else
            error()
        end
        P.draw()
    end
    return new_top, new_bottom
end

"""
    how_close(v1, v2, v_old, weightdiv=5)

Squared distance of three points, the distance to the last has only
1/weightdiv the weight.  Returns also the v1 or v1 which is closer to
v_old.
"""
function how_close(v1, v2, v_old, weightdiv=5)
    if v_old<0
        return 10^6, -1
    else
        fit = ((v1-v_old)^2 + (v2-v_old)^2)÷weightdiv + (v2-v1)^2
        tmp = (v1-v_old)^2>(v2-v_old)^2 ? v2 : v1
        return fit, tmp
        # if v1==v2
        #     # weight closeness to v_old
        #     return sqrt((v1-v_old)^2)/5000, v1
        # else
        #     tmp = (v1-v_old)^2>(v2-v_old)^2 ? v2 : v1
        #     return (v1-v2)^2.0, tmp
        # end
    end
end

"""
    pick_next_point(v1, v2, vp, v_old)
    pick_next_point(v1, v2, vp, vpp, v_old)

Try to pick the better of v1 and v2.
"""
function pick_next_point(v1, v2, vp, v_old)
    fits = zeros(2)
    for (i,v) in enumerate((v1,v2))
        # oldness = v>v_old ? (v-v_old)^2 : 2*(v-v_old)^2
        oldness = (v-v_old)^2
        fits[i] = (v-vp)^2 + 2*oldness÷3
    end
    ii = findmin(fits)[2]
    if ii==1
        return v1
    elseif ii==2
        return v2
    else # mean is best -> try harder to find which v is better
        error()
    end
    ## old:
    # if v1>=v_old && v2>=v_old
    #     # both good.  Pick the one closer to
    #     # previous point:
    #     return (v1-vp)^2>(v2-vp)^2 ? v2 : v1
    # elseif v1<v_old && v2<v_old
    #     # oh no!
    #     return v_old
    # elseif v1<v_old
    #     return v2
    # elseif v2<v_old
    #     return v1
    # else
    #     error()
    # end
end

function pick_next_point(v1, v2, vp, vpp, v_old)
    fits = zeros(2)
    dv = vpp-vp
    for (i,v) in enumerate((v1,v2))
        # oldness = v>v_old ? (v-v_old)^2 : 2*(v-v_old)^2
        oldness = (v-v_old)^2
        dvv = v-vp
        fits[i] = (dv-dvv)^2 + oldness÷5
    end
    ii = findmin(fits)[2]
    if ii==1
        return v1
    elseif ii==2
        return v2
    else
        error()
    end
end

function rchannel_stats(tdist, bdist)
    d = tdist+bdist
    dia_mean = mean(d)
    dia_std = std(d)
    dia_quant = quantile(d, [0.1, 0.5, 0.9])
    scalopping = (dia_quant[3]-dia_quant[1])/dia_quant[2]
    scalopping2 = dia_std./dia_mean
    tb_cor = cor(tdist, bdist)
    return dia_mean, dia_quant, scalopping, tb_cor
end
function rchannel_stats(tdists::Matrix, bdists::Matrix)::Tuple{
    Vector{Float64}, Matrix{Float64}, Vector{Float64}, Vector{Float64}}
    tmp = hcat([[rchannel_stats(tdists[:,i],bdists[:,i])...] for i=1:size(bdists,2)]...)
    return tmp[1,:], hcat(tmp[2,:]...), tmp[3,:], tmp[4,:]
end


text:
function snakecase(cur_string::AbstractString)
  _delimcase(cur_string, "_")
end


text:
module HOOMD

# Signal Revise.jl that this module should be tracked as a package
__revise_mode__ = :eval

# Dependecies
using CUDA
using DLPack
using PyCall
using StaticArrays

DLExt = pyimport("hoomd.dlext")

# Types
struct ContextWrapper
    context::PyObject
    sysview::PyObject
    synchronize::PyObject

    function ContextWrapper(context::PyObject)
        sysview = @pycall DLExt.SystemView(context."system_definition")::PyObject
        return new(context, sysview, sysview.synchronize)
    end
end

@pydef mutable struct Sampler <: DLExt.HalfStepHook
    function __init__(self, sample_and_bias)
        py"super"(Sampler, self).__init__()
        self.sample_and_bias = sample_and_bias
        return nothing
    end
    function update(self, timestep)
        self.sample_and_bias()
        return nothing
    end
end

# Methods
is_on_gpu(context) = pycall(context.on_gpu, Bool)

location(::Type{Array}) = DLExt.AccessLocation."OnHost"
location(::Type{CuArray}) = DLExt.AccessLocation."OnDevice"

vector(::Type{A}, ::Type{T}, v) where {A, T} = unsafe_wrap(A, DLVector{T}(v))

function matrix(::Type{A}, ::Type{T}, M) where {A, T}
    array = unsafe_wrap(A, DLMatrix{T}(M))
    # HOOMD's particle data is stored in row-major format
    return reshape(array, reverse(size(array)))
end

function take_snapshot(wrapped_context; location = default_location())
    if is_on_gpu(wrapped_context.context) && location == DLExt.AccessLocation."OnDevice"
        return take_snapshot(CuArray, wrapped_context)
    else
        return take_snapshot(Array, wrapped_context)
    end
end

function take_snapshot(A::Union{Type{Array}, Type{CuArray}}, wrapped_context)
    context = wrapped_context.context
    sysview = wrapped_context.sysview
    mode = DLExt.AccessMode."ReadWrite"
    loc = location(A)
    T = Float64

    positions = matrix(A, T, @pycall DLExt.positions_types(sysview, loc, mode)::PyObject)
    momenta = matrix(A, T, @pycall DLExt.velocities_masses(sysview, loc, mode)::PyObject)
    forces = matrix(A, T, @pycall DLExt.net_forces(sysview, loc, mode)::PyObject)
    ids = vector(A, UInt32, @pycall DLExt.rtags(sysview, loc, mode)::PyObject)

    pdata = pycall(sysview."particle_data", PyObject)
    box = pycall(pdata."getGlobalBox", PyObject)

    xy = pycall(box."getTiltFactorXY", T)
    xz = pycall(box."getTiltFactorXZ", T)
    yz = pycall(box."getTiltFactorYZ", T)
    lo = pycall(box."getLo", PyObject)
    L = pycall(box."getL", PyObject)
    Lx, Ly, Lz = L."x", L."y", L."z"

    H = @SMatrix T[
        Lx   Ly * xy   Lz * xz ;
        0    Ly        Lz * yz ;
        0    0         Lz      ;
    ]
    origin = @SVector T[lo."x", lo."y", lo."z"]

    dt = convert(T, context."integrator"."dt")

    return (  # NamedTuple
        positions  = view(positions, 1:3, :),
        velocities = view(momenta, 1:3, :),
        masses     = view(momenta, 4, :),
        forces     = view(forces, 1:3, :),
        ids        = ids,
        box        = (H = H, origin = origin),
        dt         = dt,
    )
end

function bind(context, sampling_method, callback; kwargs...)
    wrapped_context = ContextWrapper(context)
    snapshot = take_snapshot(wrapped_context)
    update = sampling_method(snapshot)
    sample_and_bias = () -> update(; sync_backend = wrapped_context.synchronize)

    sampler = Sampler(sample_and_bias)
    context."integrator"."cpp_integrator"."setHalfStepHook"(sampler)
    return sampler
end

# Module Initialization
function __init__()
    if hasproperty(DLExt.AccessLocation, :OnDevice)
        @eval default_location() = DLExt.AccessLocation."OnDevice"
    else
        @eval default_location() = DLExt.AccessLocation."OnHost"
    end
end

end  # module HOOMD