text:
#
# SAOImageDS9.jl --
#
# Implement communication with SAOImage/DS9 (http://ds9.si.edu) via the XPA
# protocol.
#
#------------------------------------------------------------------------------
#
# This file is part of SAOImageDS9.jl released under the MIT "expat" license.
# Copyright (C) 2016-2020, Éric Thiébaut (https://github.com/emmt).
#

module SAOImageDS9

export DS9

const DS9 = SAOImageDS9

using XPA
using XPA: TupleOf, connection, join_arguments

using TwoDimensional

using Base: ENV
using Base.Iterators: Pairs

# FITS pixel types.
const PixelTypes = Union{UInt8,Int16,Int32,Int64,Float32,Float64}
const PIXELTYPES = (UInt8, Int16, Int32, Int64, Float32, Float64)

"""
    SAOImageDS9.accesspoint()

yields the XPA access point which identifies the SAOImage/DS9 server.  This
access point can be set by calling the [`SAOImageDS9.connect`](@ref) method.
An empty string is returned if no access point has been chosen.  To
automatically connect to SAOImage/DS9 if not yet done, you can do:

    if SAOImageDS9.accesspoint() == ""; SAOImageDS9.connect(); end

See also [`SAOImageDS9.connect`](@ref) and [`SAOImageDS9.accesspoint`](@ref).

""" accesspoint
const _ACCESSPOINT = Ref("")
accesspoint() = _ACCESSPOINT[]

# Same as `connection()` but attempt to automatically connect an access point
# has not yet been chosen.
function _apt() # @btime -> 3.706 ns (0 allocations: 0 bytes)
    apt = accesspoint()
    if apt == ""
        try
            apt = connect()
        catch err
            _warn("Failed to automatically connect to SAOImage/DS9.\n",
                  "Launch ds9 then do `ds9connect()`")
        end
    end
    return apt
end

"""
    SAOImageDS9.connect(ident="DS9:*") -> apt

set the access point for further SAOImage/DS9 commands.  Argument `ident`
identifies the XPA access point, it can be a template string like `"DS9:*"`
which is the default value or a regular expression.  The returned value is the
name of the access point.

To retrieve the name of the current SAOImage/DS9 access point, call the
[`SAOImageDS9.accesspoint`](@ref) method.

"""
function connect(ident::Union{Regex,AbstractString} = "DS9:*"; kwds...)
    apt = XPA.find(ident; kwds...)
    apt === nothing && error("no matching SAOImage/DS9 server found")
    rep = XPA.get(apt, "version"; nmax=1)
    if length(rep) != 1 || ! XPA.verify(rep)
        error("XPA server at address \"" * apt *
              "\" does not seem to be a SAOImage/DS9 server")
    end
    addr = XPA.address(apt)
    _ACCESSPOINT[] = addr
    return addr
end

_warn(args...) = printstyled(stderr, "WARNING: ", args..., "\n";
                             color=:yellow)

"""
    SAOImageDS9.get([T, [dims,]] args...)

sends a "get" request to the SAOImage/DS9 server.  The request is made of
arguments `args...` converted into strings and merged with separating spaces.
An exception is thrown in case of error.

The returned value depends on the optional arguments `T` and `dims`:

* If neither `T` nor `dims` are specified, an instance of `XPA.Reply` is
  returned with at most one answer (see documentation for `XPA.get` for more
  details).

* If only `T` is specified, it can be:

  * `String` to return the answer as a single string;

  * `Vector{String}}` or `Tuple{Vararg{String}}` to return the answer split in
    words as a vector or as a tuple of strings;

  * `T` where `T<:Real` to return a value of type `T` obtained by parsing the
    textual answer.

  * `Tuple{Vararg{T}}` where `T<:Real` to return a value of type `T` obtained
    by parsing the textual answer;

  * `Vector{T}` where `T` is not `String` to return the binary contents
    of the answer as a vector of type `T`;

* If both `T` and `dims` are specified, `T` can be an array type like
  `Array{S}` or `Array{S,N}` and `dims` a list of `N` dimensions to retrieve
  the binary contents of the answer as an array of type `Array{S,N}`.

As a special case:

    SAOImageDS9.get(Array; endian=:native) -> arr

yields the contents of current SAOImage/DS9 frame as an array (or as `nothing`
if the frame is empty). Keyword `endian` can be used to specify the byte order
of the received values (see [`SAOImageDS9.byte_order`](@ref)).

To retrieve the version of the SAOImage/DS9 program:

    SAOImageDS9.get(VersionNumber)

See also [`SAOImageDS9.connect`](@ref), [`SAOImageDS9.set`](@ref) and
`XPA.get`.

"""
get(args...) = XPA.get(_apt(), join_arguments(args); nmax=1, throwerrors=true)

# Yields result as a vector of numerical values extracted from the binary
# contents of the reply.
get(::Type{Vector{T}}, args...) where {T} =
    XPA.get(Vector{T}, _apt(), join_arguments(args); nmax=1, throwerrors=true)

# Idem with given number of elements.
get(::Type{Vector{T}}, dim::Integer, args...) where {T} =
    XPA.get(Vector{T}, (dim,), _apt(), join_arguments(args);
            nmax=1, throwerrors=true)

# Yields result as an array of numerical values with given dimensions
# and extracted from the binary contents of the reply.
get(::Type{Array{T,N}}, dims::NTuple{N,Integer}, args...) where {T,N} =
    XPA.get(Array{T,N}, dims, _apt(), join_arguments(args);
            nmax=1, throwerrors=true)

# Idem but Array number of dimensions not specified.
get(::Type{Array{T}}, dims::NTuple{N,Integer}, args...) where {T,N} =
    get(Array{T,N}, dims, args...)

# Yields result as a single string.
get(::Type{String}, args...) =
    XPA.get(String, _apt(), join_arguments(args); nmax=1, throwerrors=true)

# Yields result as a vector of strings split out of the textual contents of the
# reply.
function get(::Type{Vector{String}}, args...;
             delim = isspace,
             keepempty::Bool=false)
    return split(chomp(get(String, args...)), delim; keepempty=keepempty)
end

# Yields result as a tuple of strings split out of the textual contents of the
# reply.
get(::Type{TupleOf{String}}, args...; kdws...) =
    Tuple(get(Vector{String}, args...; kdws...))

# Yields result as a numerical value parsed from the textual contents of the
# reply.
function get(::Type{T}, args...) :: T where {T<:Real}
    _parse(T, get(String, args...))
end

# Yields result as a tuple of numerical values parsed from the textual contents
# of the reply.
get(::Type{TupleOf{T}}, args...) where {T<:Real} =
    _parse(TupleOf{T}, get(String, args...))

_parse(::Type{T}, str::AbstractString) where {T<:Real} = parse(T, str)

function _parse(::Type{Bool}, str::AbstractString)
    key = strip(str)
    haskey(_BOOLEANS, key) ||
        throw(ArgumentError("invalid boolean textual value"))
    return _BOOLEANS[key]
end

const _BOOLEANS = Dict{String,Bool}("true"  => true, "yes" => true,
                                    "false" => false, "no" => false)

function _parse(::Type{TupleOf{T}},
                list::TupleOf{AbstractString}) where {T<:Real}
    return map(s -> _parse(T, s), list)
end

function _parse(::Type{TupleOf{T}},
                str::AbstractString) where {T<:Real}
    return Tuple(map(s -> _parse(T, s), split(str; keepempty=false)))
end

function get(::Type{Array}; endian::Union{Symbol,AbstractString}=:native)
    T = bitpix_to_type(get(Int, "fits bitpix"))
    if T === Nothing; return nothing; end
    dims = get(TupleOf{Int}, "fits size")
    return get(Array{T}, dims, "array", byte_order(endian))
end

function get(::Type{VersionNumber})
    str = get(String, "version")
    m = match(r"^ds9 +([.0-9]+[a-z]*)\s*$", str)
    m === nothing && error("unknown version number in \"$str\"")
    return VersionNumber(m.captures[1])
end

"""
    SAOImageDS9.set(args...; data=nothing)

sends command and/or data to the SAOImage/DS9 server.  The command is made of
arguments `args...` converted into strings and merged with a separating spaces.
Keyword `data` can be used to specify the data to send.  An exception is thrown
in case of error.

As a special case:

    SAOImageDS9.set(arr; mask=false, new=false, endian=:native)

set the contents of the current SAOImage/DS9 frame to be array `arr`.  Keyword
`new` can be set true to create a new frame for displyaing the array.  Keyword
`endian` can be used to specify the byte order of the values in `arr` (see
[`SAOImageDS9.byte_order`](@ref).

See also [`SAOImageDS9.connect`](@ref), [`SAOImageDS9.get`](@ref) and
`XPA.set`.

"""
function set(args...; data=nothing)
    XPA.set(_apt(), join_arguments(args); nmax=1, throwerrors=true, data=data)
    return nothing
end

function set(arr::DenseArray{T,N};
             endian::Symbol=:native,
             mask::Bool=false,
             new::Bool=false) where {T<:PixelTypes,N}
    args = String[]
    push!(args, "array")
    if new; push!(args, "new"); end
    if mask; push!(args, "mask"); end
    set(args..., _arraydescriptor(arr; endian=endian); data=arr)
end

# Convert other pixel types.
for (T, S) in ((Int8,   Int16),
               (UInt16, Float32),
               (UInt32, Float32),
               (UInt64, Float32))
    @eval set(arr::AbstractArray{$T,N}; kwds...) where {N} =
        set(convert(Array{$S,N}, arr); kwds...)
end

# Convert non-dense array types.
for T in PIXELTYPES
    @eval set(arr::AbstractArray{$T,N}; kwds...) where {N} =
        set(convert(Array{$T,N}, arr); kwds...)
end

function _arraydescriptor(arr::DenseArray{T,2};
                          endian::Symbol=:native) where {T}
    bp = bitpix_of(T)
    bp != 0 || error("unsupported data type")
    return string("[xdim=",size(arr,1),",ydim=",size(arr,2),
                  ",bitpix=",bp,",endian=",endian,"]")
end

function _arraydescriptor(arr::DenseArray{T,3};
                          endian::Symbol=:native) where {T}
    bp = bitpix_of(T)
    bp != 0 || error("unsupported data type")
    return string("[xdim=",size(arr,1),",ydim=",size(arr,2),",zdim=",size(arr,3),
                  ",bitpix=",bp,",endian=",endian,"]")
end

_arraydescriptor(arr::DenseArray; kdws...) =
    error("only 2D and 3D arrays are supported")

"""
    SAOImageDS9.bitpix_of(x) -> bp

yields FITS bits-per-pixel (BITPIX) value for `x` which can be an array or a
type.  A value of 0 is returned if `x` is not of a supported type.

See also [`SAOImageDS9.bitpix_to_type`](@ref).

"""
bitpix_of(::DenseArray{T}) where {T} = bitpix_of(T)
for T in PIXELTYPES
    bp = (T <: Integer ? 8 : -8)*sizeof(T)
    @eval bitpix_of(::Type{$T}) = $bp
    @eval bitpix_of(::$T) = $bp
end
bitpix_of(::Any) = 0

"""
    SAOImageDS9.bitpix_to_type(bp) -> T

yields Julia type corresponding to FITS bits-per-pixel (BITPIX) value `bp`.
The value `Nothing` is returned if `bp` is unknown.

See also [`SAOImageDS9.bitpix_of`](@ref).

"""
bitpix_to_type(bitpix::Int) =
    (bitpix ==   8 ? UInt8   :
     bitpix ==  16 ? Int16   :
     bitpix ==  32 ? Int32   :
     bitpix ==  64 ? Int64   :
     bitpix == -32 ? Float32 :
     bitpix == -64 ? Float64 : Nothing)
bitpix_to_type(bitpix::Integer) = bitpix_to_type(Int(bitpix))
bitpix_to_type(::Any) = Nothing

"""
    SAOImageDS9.byte_order(endian)

yields the byte order for retrieving the elements of a SAOImage/DS9 array.
Argument can be one of the strings (or the equivalent symbol): `"big"` for most
significant byte first, `"little"` for least significant byte first or
`"native"` to yield the byte order of the machine.

See also [`SAOImageDS9.get`](@ref), [`SAOImageDS9.set`](@ref).

"""
byte_order(endian::Symbol) =
    (endian == :native ? (ENDIAN_BOM == 0x01020304 ? "big" :
                          ENDIAN_BOM == 0x04030201 ? "little" :
                          error("unknown byte order")) :
     endian == :big ? "big" :
     endian == :little ? "little" :
     error("invalid byte order"))

byte_order(endian::AbstractString) =
    (endian == "native" ? (ENDIAN_BOM == 0x01020304 ? "big" :
                           ENDIAN_BOM == 0x04030201 ? "little" :
                           error("unknown byte order")) :
     endian == "big" || endian == "little" ? endian :
     error("invalid byte order"))

#------------------------------------------------------------------------------
# DRAWING

"""
    SAOImageDS9.draw(args...; kwds...)

draws something in SAOImage/DS9 application.  The operation depends on the type
of the arguments.

---
    SAOImageDS9.draw(img; kwds...)

displays image `img` (a 2-dimensional Julia array) in SAOImage/DS9.
The following keywords are possible:

- Keyword `frame` can be used to specify the frame number.

- Keyword `cmap` can be used to specify the name of the colormap.  For
  instance, `cmap="gist_stern"`.

- Keyword `zoom` can be used to specify the zoom factor.

- Keywords `min` and/or `max` can be used to specify the scale limits.

---
    SAOImageDS9.draw(pnt; kwds...)

draws `pnt` as point(s) in SAOImage/DS9, `pnt` is a `Point`, an array or a
tuple of `Point`.

---
    SAOImageDS9.draw(box; kwds...)

draws `box` as rectangle(s) in SAOImage/DS9, `box` is a `BoundingBox`, an array
or a tuple of `BoundingBox`.

"""
draw(args...; kwds...) = draw(args; kwds...)
draw(::Tuple{}; kwds...) = nothing
draw(::T; kwds...) where {T} = error("unexpected type of argument(s): $T")

function draw(img::AbstractMatrix;
              min::Union{Real,Nothing} = nothing,
              max::Union{Real,Nothing} = nothing,
              cmap = nothing,
              frame = nothing,
              zoom = nothing)
    # FIXME: pack all commands into a single one.
    frame === nothing || ds9set("frame", frame)
    zoom === nothing || ds9set("zoom to", zoom)
    set(img)
    if min !== nothing || max !== nothing
        set("scale limits", limits(img, min, max)...)
    end
    cmap === nothing || set("cmap", cmap)
    return nothing
end

# For multiple points/circles/... we just send a command to SAOImage/DS9 for
# each item to draw.  Packing multiple commands (separated by semi-columns)
# does not really speed-up things and is more complicated because there is a
# limit to the total length of an XPA command (given by XPA.SZ_LINE I guess).

draw(A::Point; kwds...) = _draw(_region(Val(:point), kwds), A)
function draw(A::Union{Tuple{Vararg{Point}},
                       AbstractArray{<:Point}}; kwds...)
    cmd = _region(Val(:point), kwds)
    for i in eachindex(A)
        _draw(cmd, A[i])
    end
end

draw(A::BoundingBox; kwds...) =  _draw(_region(Val(:polygon), kwds), A)
function draw(A::Union{Tuple{Vararg{BoundingBox}},
                       AbstractArray{<:BoundingBox}}; kwds...)
    cmd = _region(Val(:polygon), kwds)
    for i in eachindex(A)
        _draw(cmd, A[i])
    end
end

function _draw(cmd::NTuple{2,AbstractString}, A::Point)
    set(cmd[1], A.x, A.y, cmd[2])
    nothing
end

function _draw(cmd::NTuple{2,AbstractString}, A::BoundingBox)
    x0, x1, y0, y1 = Tuple(A)
    set(cmd[1], x0, y0, x1, y0, x1, y1, x0, y1, x0, y0, cmd[2])
    nothing
end

_set_options(kwds::Pairs) = (isempty(kwds) ? "" :
                             throw(ArgumentError("unexpected argument")))

function _set_options(kwds::Pairs{Symbol})
    io = IOBuffer()
    print(io, " #")
    for (key, val) in kwds
        _set_option(io, key, val)
    end
    print(io, " }")
    return String(take!(io))
end

_set_option(io::IOBuffer, name::Union{Symbol,AbstractString}, value::Nothing) =
    nothing
_set_option(io::IOBuffer, name::Union{Symbol,AbstractString}, value) =
    print(io, " ", name, "=", value)
_set_option(io::IOBuffer, name::Union{Symbol,AbstractString}, value::Bool) =
    print(io, " ", name, (value ? "=1" : "=0"))

for id in (:circle, :ellipse, :box, :polygon, :point, :line,
           :vector, :text, :ruler, :compass, :projection, :annulus,
           :panda, :epanda, :bpanda)
    V = Val{id}
    cmd = "regions command { $id"
    @eval _region(::$V, kwds::Pairs) = ($cmd, _set_options(kwds))
end

"""
    SAOImageDS9.message([apt=SAOImageDS9.accesspoint(),] msg; cancel=false)

displays a message dialog with text `msg` in SAOImage/DS9 application referred
by `apt` and returns a boolean.  If keyword `cancel` is true, a *Cancel* button
is added to the dialog and `false` maybe returned if the dialog is not closed
by clicking the *OK* button; otherwise `true` is returned.

"""
message(msg::AbstractString; kwds...) =  message(_apt(), msg; kwds...)
function message(apt, msg::AbstractString; cancel::Bool = false)
    btn = (cancel ? "okcancel" : "ok")
    cmd = "analysis message $btn {$msg}"
    tryparse(Int, XPA.get(String, apt, cmd)) == 1
end

"""
    SAOImageDS9.select([apt=SAOImageDS9.accesspoint(),];
                       text="", key=false, cancel=false) -> (k,x,y,v)

returns the position selected by the user in SAOImage/DS9 application referred
by `apt`.  If keyword `text` is set, a dialog message is first displayed,
possibly with a *Cancel* button if keyword `cancel` is true.  If keyword `key`
is true, the position is selected by pressing a key; otherwise, the position is
selected by clicking the first mouse button.  The result is either `nothing`
(for instance if the *Cancel* button of the dialog is clicked) or a 4-tuple
`(k,x,y,v)` with `k` the pressed key (an empty string if `key` is false),
`(x,y)` are the coordinates of the selected position and `v` is the
corresponding value in the data.

"""
function select(apt = _apt();
                text::AbstractString = "",
                cancel::Bool = false,
                key::Bool = false)
    if length(text) > 0
        message(apt, text; cancel=cancel) || return nothing
    end
    cmd = (key ? "iexam key {\$x \$y \$value}" : "iexam {\$x \$y \$value}")
    rep = split(XPA.get(String, apt, cmd))
    off = (key ? 1 : 0)
    length(rep) == 3+off || return nothing
    k = (key ? string(rep[off+0]) : "")
    x = tryparse(Float64, rep[off+1])
    y = tryparse(Float64, rep[off+2])
    v = tryparse(Float64, rep[off+3])
    (x === nothing || y === nothing || v === nothing) && return nothing
    return k, x, y, v
end

#------------------------------------------------------------------------------
# LIMITS

"""
    limits(A, cmin=nothing, cmax=nothing) -> (lo, hi)

yields the clipping limits of values in array `A`.  The result is a 2-tuple of
double precision floats `(lo,hi)`.  If `cmin` is `nothing`, `lo` is the minimal
finite value found in `A` and converted to `Cdouble`; otherwise `lo =
Cdouble(cmin)`.  If `cmax` is `nothing`, `hi` is the maximal finite value found
in `A` and converted to `Cdouble`; otherwise `hi = Cdouble(cmax)`.

"""
limits(A::AbstractArray{<:Real}, ::Nothing, ::Nothing) =
    to_limits(finite_extrema(A))

limits(A::AbstractArray{<:Real}, cmin::Real, ::Nothing) =
    to_limits(cmin, finite_maximum(A))

limits(A::AbstractArray{<:Real}, ::Nothing, cmax::Real) =
    to_limits(finite_minimum(A), cmax)

limits(A::AbstractArray{<:Real}, cmin::Real, cmax::Real) =
    to_limits(cmin, cmax)

to_limits(c::Tuple{Real,Real}) = to_limits(c...)
to_limits(cmin::Cdouble, cmax::Cdouble) = (cmin, cmax)
to_limits(cmin::Real, cmax::Real) = to_limits(convert(Cdouble, cmin),
                                              convert(Cdouble, cmax))

"""
    finite_extrema(A) -> (vmin, vmax)

yields the minimum and maximum finite values in array `A`.  The result is such
that `vmin ≤ vmax` (both values being finite) unless there are no finite values
in `A` in which case `vmin > vmax`.

"""
finite_extrema(A::AbstractArray{<:Real}) = valid_extrema(isfinite, A)

"""
    finite_minimum(A) -> vmin

yields the minimum finite value in array `A`.  The result is never a NaN but
may be `typemax(eltype(A))` if there are no finite values in `A`.

"""
finite_minimum(A::AbstractArray{<:Real}) = valid_minimum(isfinite, A)

"""
    finite_maximum(A) -> vmax

yields the maximum finite value in array `A`.  The result is never a NaN but
may be `typemin(eltype(A))` if there are no finite values in `A`.

"""
finite_maximum(A::AbstractArray{<:Real}) = valid_maximum(isfinite, A)

"""
    valid_extrema(pred, A) -> (vmin, vmax)

yields the minimum and maximum valid values in array `A`.  Valid values are
those for which predicate `pred` yields `true`.  The result is such that `vmin
≤ vmax` (both values being valid) unless there are no valid values in `A` in
which case `vmin > vmax`.  The predicate function is assumed to take care of
NaN's.

"""
function valid_extrema(pred, A::AbstractArray{T}) where {T}
    vmin = typemax(T)
    vmax = typemin(T)
    @inbounds @simd for i in eachindex(A)
        val = A[i]
        vmin = ifelse(pred(val)&(val < vmin), val, vmin)
        vmax = ifelse(pred(val)&(val > vmax), val, vmax)
    end
    return (vmin, vmax)
end
valid_extrema(::typeof(isfinite), A::AbstractRange{<:Real}) = extrema(A)

"""
    valid_minimum(A) -> vmin

yields the minimum valid value in array `A`.  Valid values are those for which
predicate `pred` yields `true`.  The result is a valid value but may be
`typemax(eltype(A))` if there are no valid values in `A`.  The predicate
function is assumed to take care of NaN's.

"""
function valid_minimum(pred, A::AbstractArray{T}) where {T}
    vmin = typemax(T)
    @inbounds @simd for i in eachindex(A)
        val = A[i]
        vmin = ifelse(pred(val)&(val < vmin), val, vmin)
    end
    return vmin
end
valid_minimum(::typeof(isfinite), A::AbstractRange{<:Real}) = minimum(A)

"""
    valid_maximum(A) -> vmax

yields the maximum valid value in array `A`.  Valid values are those for which
predicate `pred` yields `true`.  The result is a valid value but may be
`typemin(eltype(A))` if there are no valid values in `A`.  The predicate
function is assumed to take care of NaN's.

"""
function valid_maximum(pred, A::AbstractArray{T}) where {T}
    vmax = typemin(T)
    @inbounds @simd for i in eachindex(A)
        val = A[i]
        vmax = ifelse(pred(val)&(val > vmax), val, vmax)
    end
    return vmax
end
valid_maximum(::typeof(isfinite), A::AbstractRange{<:Real}) = maximum(A)

end # module


text:
# Fast parallel reduction for Kepler hardware
# - uses shuffle and shared memory to reduce efficiently
# - support for large arrays
#
# Based on devblogs.nvidia.com/parallelforall/faster-parallel-reductions-kepler/

using CUDAdrv, CUDAnative

#
# Main implementation
#

# Reduce a value across a warp
@inline function reduce_warp{F<:Function,T}(op::F, val::T)::T
    offset = CUDAnative.warpsize() ÷ UInt32(2)
    # TODO: this can be unrolled if warpsize is known...
    while offset > 0
        val = op(val, shfl_down(val, offset))
        offset ÷= UInt32(2)
    end
    return val
end

# Reduce a value across a block, using shared memory for communication
@inline function reduce_block{F<:Function,T}(op::F, val::T)::T
    # shared mem for 32 partial sums
    shared = @cuStaticSharedMem(T, 32)

    # TODO: use fldmod1 (JuliaGPU/CUDAnative.jl#28)
    wid  = div(threadIdx().x-UInt32(1), CUDAnative.warpsize()) + UInt32(1)
    lane = rem(threadIdx().x-UInt32(1), CUDAnative.warpsize()) + UInt32(1)

    # each warp performs partial reduction
    val = reduce_warp(op, val)

    # write reduced value to shared memory
    if lane == 1
        @inbounds shared[wid] = val
    end

    # wait for all partial reductions
    sync_threads()

    # read from shared memory only if that warp existed
    @inbounds val = (threadIdx().x <= fld(blockDim().x, CUDAnative.warpsize())) ? shared[lane] : zero(T)

    # final reduce within first warp
    if wid == 1
        val = reduce_warp(op, val)
    end

    return val
end

# Reduce an array across a complete grid
function reduce_grid{F<:Function,T}(op::F, input::CuDeviceArray{T,1},
                                    output::CuDeviceArray{T,1}, len::Integer)

    # TODO: neutral element depends on the operator (see Base's 2 and 3 argument `reduce`)
    val = zero(T)

    # reduce multiple elements per thread (grid-stride loop)
    # TODO: step range (see JuliaGPU/CUDAnative.jl#12)
    i = (blockIdx().x-UInt32(1)) * blockDim().x + threadIdx().x
    step = blockDim().x * gridDim().x
    while i <= len
        @inbounds val = op(val, input[i])
        i += step
    end

    val = reduce_block(op, val)

    if threadIdx().x == UInt32(1)
        @inbounds output[blockIdx().x] = val
    end

    return
end

"""
Reduce a large array.

Kepler-specific implementation, ie. you need sm_30 or higher to run this code.
"""
function gpu_reduce{F<:Function,T}(op::F, input::CuArray{T,1}, output::CuArray{T,1})
    len = length(input)

    # TODO: these values are hardware-dependent, with recent GPUs supporting more threads
    threads = 512
    blocks = min((len + threads - 1) ÷ threads, 1024)

    # the output array must have a size equal to or larger than the number of thread blocks
    # in the grid because each block writes to a unique location within the array.
    if length(output) < blocks
        throw(ArgumentError("output array too small, should be at least $blocks elements"))
    end

    @cuda (blocks,threads) reduce_grid(op, input, output, Int32(len))
    @cuda (1,1024) reduce_grid(op, output, output, Int32(blocks))

    return
end

# FURTHER IMPROVEMENTS:
# - use atomic memory operations
# - dynamic block/grid size based on device capabilities
# - vectorized memory access
#   devblogs.nvidia.com/parallelforall/cuda-pro-tip-increase-performance-with-vectorized-memory-access/

ext: jl


text:
"""
$(SIGNATURES)

Create customized distinguishable colormap for interior regions
"""
region_cmap(n)=distinguishable_colors(max(5,n),
                                      [RGB(0.85,0.6,0.6), RGB(0.6,0.85,0.6),RGB(0.6,0.6,0.85)],
                                      lchoices = range(70, stop=80, length=5),
                                      cchoices = range(25, stop=65, length=15),
                                      hchoices = range(20, stop=360, length=15)
                                      )

"""
$(SIGNATURES)

Create customized distinguishable colormap for boundary regions
"""
bregion_cmap(n)=distinguishable_colors(max(5,n),
                                      [RGB(1.0,0.0,0.0), RGB(0.0,1.0,0.0), RGB(0.0,0.0,1.0)],
                                      lchoices = range(50, stop=75, length=10),
                                      cchoices = range(75, stop=100, length=10),
                                      hchoices = range(20, stop=360, length=30)
                                      )

function Colors.RGB(c::String)
    c64=Colors.color_names[c]
    RGB(c64[1]/255,c64[2]/255, c64[3]/255)
end

Colors.RGB(c::Symbol)=Colors.RGB(String(c))
Colors.RGB(c::Tuple)=Colors.RGB(c...)
rgbtuple(c)=rgbtuple(Colors.RGB(c))
rgbtuple(c::RGB)=(red(c),green(c),blue(c))

"""
$(SIGNATURES)

Extract visible tetrahedra - those intersecting with the planes
`x=xyzcut[1]` or `y=xyzcut[2]`  or `z=xyzcut[3]`.

Return corresponding points and facets for each region for drawing as mesh (Makie,MeshCat)
or trisurf (pyplot)
"""
function extract_visible_cells3D(grid::ExtendableGrid,xyzcut; primepoints=zeros(0,0),Tp=SVector{3,Float32},Tf=SVector{3,Int32})
    coord=grid[Coordinates]
    cellnodes=grid[CellNodes]
    cellregions=grid[CellRegions]
    nregions=grid[NumCellRegions]
    extract_visible_cells3D(coord,cellnodes,cellregions,nregions,xyzcut;
                            primepoints=primepoints,
                            Tp=Tp,Tf=Tf)
end

function extract_visible_cells3D(coord,cellnodes,cellregions,nregions,xyzcut;
                                 primepoints=zeros(0,0),Tp=SVector{3,Float32},Tf=SVector{3,Int32})

    function take(coord,simplex,xyzcut)
        all_lt=@MVector ones(Bool,3)
        all_gt=@MVector ones(Bool,3)
        for idim=1:3
            for inode=1:4
                c=coord[idim,simplex[inode]]-xyzcut[idim]
                all_lt[idim]=all_lt[idim] && (c<0.0)
                all_gt[idim]=all_gt[idim] && (c>0.0)
            end
        end
        tke=false
        tke=tke  ||   (!all_lt[1])  &&  (!all_gt[1]) && (!all_gt[2]) && (!all_gt[3])
        tke=tke  ||   (!all_lt[2])  &&  (!all_gt[2]) && (!all_gt[1]) && (!all_gt[3])
        tke=tke  ||   (!all_lt[3])  &&  (!all_gt[3]) && (!all_gt[1]) && (!all_gt[2])
    end

    faces=[Vector{Tf}(undef,0) for iregion=1:nregions]
    points=[Vector{Tp}(undef,0) for iregion=1:nregions]

    for iregion=1:nregions
        for iprime=1:size(primepoints,2)
            @views push!(points[iregion],Tp(primepoints[:,iprime]))
        end
    end
    tet=zeros(Int32,4)

    for itet=1:size(cellnodes,2)
        iregion=cellregions[itet]
        for i=1:4
            tet[i]=cellnodes[i,itet]
        end
        if take(coord,tet,xyzcut)
            npts=size(points[iregion],1)
            @views begin
                push!(points[iregion],coord[:,cellnodes[1,itet]])
                push!(points[iregion],coord[:,cellnodes[2,itet]])
                push!(points[iregion],coord[:,cellnodes[3,itet]])
                push!(points[iregion],coord[:,cellnodes[4,itet]])
                push!(faces[iregion],(npts+1,npts+2,npts+3))
                push!(faces[iregion],(npts+1,npts+2,npts+4))
                push!(faces[iregion],(npts+2,npts+3,npts+4))
                push!(faces[iregion],(npts+3,npts+1,npts+4))
            end
        end
    end
    points,faces
end

function extract_visible_bfaces3D(grid::ExtendableGrid,xyzcut; primepoints=zeros(0,0), Tp=SVector{3,Float32},Tf=SVector{3,Int32})
    coord=grid[Coordinates]
    bfacenodes=grid[BFaceNodes]
    bfaceregions=grid[BFaceRegions]
    nbregions=grid[NumBFaceRegions]

    extract_visible_bfaces3D(coord,bfacenodes,bfaceregions, nbregions, xyzcut;
                             primepoints=primepoints,Tp=Tp,Tf=Tf)
end

function extract_visible_bfaces3D(coord,bfacenodes,bfaceregions, nbregions, xyzcut;
                                  primepoints=zeros(0,0), Tp=SVector{3,Float32},Tf=SVector{3,Int32})

    nbfaces=size(bfacenodes,2)
    cutcoord=zeros(3)

    function take(coord,simplex,xyzcut)
        for idim=1:3
            all_gt=true
            for inode=1:3
                c=coord[idim,simplex[inode]]-xyzcut[idim]
                all_gt= all_gt && c>0
            end
            if all_gt
                return false
            end
        end
        return true
    end

    Tc=SVector{3,eltype(coord)}
    xcoord=reinterpret(Tc,reshape(coord,(length(coord),)))

    faces=[Vector{Tf}(undef,0) for iregion=1:nbregions]
    points=[Vector{Tp}(undef,0) for iregion=1:nbregions]
    for iregion=1:nbregions
        for iprime=1:size(primepoints,2)
            @views push!(points[iregion],Tp(primepoints[:,iprime]))
        end
    end

    # remove some type instability here
    function collct(points,faces)
        trinodes=[1,2,3]
        for i=1:nbfaces
            iregion=bfaceregions[i]
            trinodes[1]=bfacenodes[1,i]
            trinodes[2]=bfacenodes[2,i]
            trinodes[3]=bfacenodes[3,i]
            if take(coord,trinodes,xyzcut)
                npts=size(points[iregion],1)
                @views push!(points[iregion],xcoord[trinodes[1]])
                @views push!(points[iregion],xcoord[trinodes[2]])
                @views push!(points[iregion],xcoord[trinodes[3]])
                @views push!(faces[iregion],(npts+1,npts+2,npts+3))
            end
        end
    end
    collct(points,faces)
    points,faces
end

# old version with function values
function extract_visible_bfaces3D(grid::ExtendableGrid,func,xyzcut)
    cutcoord=zeros(3)

    function take(coord,simplex,xyzcut)
        for idim=1:3
            for inode=1:3
                cutcoord[inode]=coord[idim,simplex[inode]]-xyzcut[idim]
            end
            if !mapreduce(a->a<=0,*,cutcoord)
                return false
            end
        end
        return true
    end

    coord=grid[Coordinates]
    nbfaces=num_bfaces(grid)
    bfacenodes=grid[BFaceNodes]

    pmark=zeros(UInt32,size(coord,2))
    faces=ElasticArray{UInt32}(undef,3,0)
    npoints=0

    for i=1:nbfaces
        tri=view(bfacenodes,:, i)
        if take(coord,tri,xyzcut)
            for inode=1:3
                if pmark[tri[inode]]==0
                    npoints+=1
                    pmark[tri[inode]]=npoints
                end
            end
            tri=map(i->pmark[i],tri)
            append!(faces,tri)
        end
    end

    points=Array{Float32,2}(undef,3,npoints)
    values=Vector{Float32}(undef,npoints)
    for i=1:size(coord,2)
        if pmark[i]>0
            @views points[:,pmark[i]].=(coord[1,i],coord[2,i],coord[3,i])
            values[pmark[i]]=func[i]
        end
    end
    points,faces,values
end

"""
  $(SIGNATURES)
  Calculate intersections between tetrahedron with given piecewise linear
  function data and plane

  Adapted from https://github.com/j-fu/gltools/blob/master/glm-3d.c#L341

  A non-empty intersection is either a triangle or a planar quadrilateral,
  define by either 3 or 4 intersection points between tetrahedron edges
  and the plane.

  Input:
  -       pointlist: 3xN array of grid point coordinates
  -    node_indices: 4 element array of node indices (pointing into pointlist and function_values)
  -   planeq_values: 4 element array of plane equation evaluated at the node coordinates
  - function_values: N element array of function values

  Mutates:
  -  ixcoord: 3x4 array of plane - tetedge intersection coordinates
  - ixvalues: 4 element array of fuction values at plane - tetdedge intersections

  Returns:
  - nxs,ixcoord,ixvalues

  This method can be used both for the evaluation of plane sections and for
  the evaluation of function isosurfaces.
"""
function tet_x_plane!(ixcoord,ixvalues,pointlist,node_indices,planeq_values,function_values; tol=0.0)

    # If all nodes lie on one side of the plane, no intersection
    if (mapreduce(a->a< -tol,*,planeq_values) || mapreduce(a->a>tol,*,planeq_values))
        return 0
    end
    # Interpolate coordinates and function_values according to
    # evaluation of the plane equation
    nxs=0
    @inbounds @simd for n1=1:3
        N1=node_indices[n1]
        @inbounds @fastmath @simd for n2=n1+1:4
            N2=node_indices[n2]
            if planeq_values[n1]*planeq_values[n2]<tol
                nxs+=1
                t= planeq_values[n1]/(planeq_values[n1]-planeq_values[n2])
                ixcoord[1,nxs]=pointlist[1,N1]+t*(pointlist[1,N2]-pointlist[1,N1])
                ixcoord[2,nxs]=pointlist[2,N1]+t*(pointlist[2,N2]-pointlist[2,N1])
                ixcoord[3,nxs]=pointlist[3,N1]+t*(pointlist[3,N2]-pointlist[3,N1])
                ixvalues[nxs]=function_values[N1]+t*(function_values[N2]-function_values[N1])
            end
        end
    end
    return nxs
end

"""
 We should be able to parametrize this
 with a pushdata function which will remove one copy
 step for GeometryBasics.mesh creation - perhaps a meshcollector struct we
 can dispatch on.
 flevel could be flevels
 xyzcut could be a vector of plane data
 perhaps we can also collect isolines.
 Just an optional collector parameter, defaulting to somethig makie independent.

    Better yet:

 struct TetrahedronMarcher
  ...
 end
 tm=TetrahedronMarcher(planes,levels)

 foreach tet
       collect!(tm, tet_node_coord, node_function_values)
 end
 tm.colors=AbstractPlotting.interpolated_getindex.((cmap,), mcoll.vals, (fminmax,))
 mesh!(collect(mcoll),backlight=1f0)
"""

"""
   $(SIGNATURES)

"""
function marching_tetrahedra(grid::ExtendableGrid,func,planes,flevels; kwargs...)
    coord=grid[Coordinates]
    cellnodes=grid[CellNodes]
    marching_tetrahedra(coord,cellnodes,func,planes,flevels;kwargs...)
end

"""
   $(SIGNATURES)

"""
function marching_tetrahedra(coord,cellnodes,func,planes,flevels;
                             tol=1.0e-12,
                             primepoints=zeros(0,0),
                             primevalues=zeros(0),
                             Tv=Float32,
                             Tp=SVector{3,Float32},
                             Tf=SVector{3,Int32})
    # We could rewrite this for Meshing.jl
    # CellNodes::Vector{Ttet}, Coord::Vector{Tpt}
    nplanes=length(planes)
    nlevels=length(flevels)
    nnodes=size(coord,2)
    ntet=size(cellnodes,2)

    all_planeq=Vector{Float32}(undef,nnodes)

    all_ixfaces=Vector{Tf}(undef,0)
    all_ixcoord=Vector{Tp}(undef,0)
    all_ixvalues=Vector{Tv}(undef,0)

    @assert(length(primevalues)==size(primepoints,2))
    for iprime=1:size(primepoints,2)
        @views push!(all_ixcoord,primepoints[:,iprime])
        @views push!(all_ixvalues,primevalues[iprime])
    end

    planeq=zeros(4)
    ixcoord=zeros(3,6)
    ixvalues=zeros(6)
    cn=zeros(4)
    node_indices=zeros(Int32,4)

    @inbounds @fastmath plane_equation(plane,coord)= coord[1]*plane[1]+coord[2]*plane[2]+coord[3]*plane[3]+plane[4]

    function pushtris(ns,ixcoord,ixvalues)
        # number of intersection points can be 3 or 4
        if ns>=3
            last_i=length(all_ixvalues)
            for is=1:ns
                @views push!(all_ixcoord,ixcoord[:,is])
                push!(all_ixvalues,ixvalues[is])
            end
            push!(all_ixfaces,(last_i+1,last_i+2,last_i+3))
            if ns==4
                push!(all_ixfaces,(last_i+3,last_i+2,last_i+4))
            end
        end
    end

    function calcxs()
        @inbounds for itet=1:ntet
            node_indices[1]=cellnodes[1,itet]
            node_indices[2]=cellnodes[2,itet]
            node_indices[3]=cellnodes[3,itet]
            node_indices[4]=cellnodes[4,itet]
            planeq[1]=all_planeq[node_indices[1]]
            planeq[2]=all_planeq[node_indices[2]]
            planeq[3]=all_planeq[node_indices[3]]
            planeq[4]=all_planeq[node_indices[4]]
            nxs=tet_x_plane!(ixcoord,ixvalues,coord,node_indices,planeq,func,tol=tol)
            pushtris(nxs,ixcoord,ixvalues)
        end
    end

    @inbounds for iplane=1:nplanes
        @views @inbounds map!(inode->plane_equation(planes[iplane],coord[:,inode]),all_planeq,1:nnodes)
        calcxs()
    end

    # allocation free (besides push!)
    @inbounds for ilevel=1:nlevels
        @views @inbounds @fastmath map!(inode->(func[inode]-flevels[ilevel]),all_planeq,1:nnodes)
        calcxs()
    end

    all_ixcoord, all_ixfaces, all_ixvalues
end

########################################################################################
"""
    $(SIGNATURES)

    Collect isoline snippets on triangles ready for linesegments!

"""

function marching_triangles(grid::ExtendableGrid,func,levels)
    coord::Matrix{Float64}=grid[Coordinates]
    cellnodes::Matrix{Int32}=grid[CellNodes]
    marching_triangles(coord,cellnodes,func,levels)
end

function marching_triangles(coord,cellnodes,func,levels)
    points=Vector{Point2f}(undef,0)
    function isect(nodes)
        (i1,i2,i3)=(1,2,3)

        f=(func[nodes[1]],func[nodes[2]],func[nodes[3]])

        f[1]  <= f[2]  ?  (i1,i2) = (1,2)   : (i1,i2) = (2,1)
        f[i2] <= f[3]  ?  i3=3              : (i2,i3) = (3,i2)
        f[i1] >  f[i2] ?  (i1,i2) = (i2,i1) : nothing

        (n1,n2,n3)=(nodes[i1],nodes[i2],nodes[i3])

        dx31=coord[1,n3]-coord[1,n1]
        dx21=coord[1,n2]-coord[1,n1]
        dx32=coord[1,n3]-coord[1,n2]

        dy31=coord[2,n3]-coord[2,n1]
        dy21=coord[2,n2]-coord[2,n1]
        dy32=coord[2,n3]-coord[2,n2]

        df31 = f[i3]!=f[i1] ? 1/(f[i3]-f[i1]) : 0.0
        df21 = f[i2]!=f[i1] ? 1/(f[i2]-f[i1]) : 0.0
        df32 = f[i3]!=f[i2] ? 1/(f[i3]-f[i2]) : 0.0

        for level ∈ levels
            if  (f[i1]<=level) && (level<f[i3])
	        α=(level-f[i1])*df31
	        x1=coord[1,n1]+α*dx31
	        y1=coord[2,n1]+α*dy31

	        if (level<f[i2])
	            α=(level-f[i1])*df21
	            x2=coord[1,n1]+α*dx21
		    y2=coord[2,n1]+α*dy21
                else
	            α=(level-f[i2])*df32
	            x2=coord[1,n2]+α*dx32
	            y2=coord[2,n2]+α*dy32
                end
                push!(points,Point2f(x1,y1))
                push!(points,Point2f(x2,y2))
            end
        end
    end

    for itri=1:size(cellnodes,2)
        @views isect(cellnodes[:,itri])
    end

    points
end

##############################################
# Create meshes from grid data
function regionmesh(grid,iregion)
    coord=grid[Coordinates]
    cn=grid[CellNodes]
    cr=grid[CellRegions]
    @views points=[Point2f(coord[:,i]) for i=1:size(coord,2)]
    faces=Vector{GLTriangleFace}(undef,0)
    for i=1:length(cr)
        if cr[i]==iregion
            @views push!(faces,cn[:,i])
        end
    end
    Mesh(points,faces)
end

function bfacesegments(grid,ibreg)
    coord=grid[Coordinates]
    nbfaces=num_bfaces(grid)
    bfacenodes=grid[BFaceNodes]
    bfaceregions=grid[BFaceRegions]
    points=Vector{Point2f}(undef,0)
    for ibface=1:nbfaces
        if bfaceregions[ibface]==ibreg
            push!(points,Point2f(coord[1,bfacenodes[1,ibface]],coord[2,bfacenodes[1,ibface]]))
            push!(points,Point2f(coord[1,bfacenodes[2,ibface]],coord[2,bfacenodes[2,ibface]]))
        end
    end
    points
end

function bfacesegments3(grid,ibreg)
    coord=grid[Coordinates]
    nbfaces=num_bfaces(grid)
    bfacenodes=grid[BFaceNodes]
    bfaceregions=grid[BFaceRegions]
    points=Vector{Point3f}(undef,0)
    for ibface=1:nbfaces
        if bfaceregions[ibface]==ibreg
            push!(points,Point3f(coord[1,bfacenodes[1,ibface]],coord[2,bfacenodes[1,ibface]],0.0))
            push!(points,Point3f(coord[1,bfacenodes[2,ibface]],coord[2,bfacenodes[2,ibface]],0.0))
        end
    end
    points
end

############################################
"""
$(SIGNATURES)

Assume that `points` are nodes of a polyline.
Place `nmarkers` equidistant markers  at the polyline, under
the assumption that the points are transformed via the transformation
matrix M vor visualization.
"""
function markerpoints(points,nmarkers,transform)
    dist(p1,p2)=norm(transform*(p1-p2))

    llen=0.0
    for i=2:length(points)
       llen+=dist(points[i],points[i-1])
    end

    mdist=llen/(nmarkers-1)

    mpoints=[points[1]]

    i=2
    l=0.0
    lnext=l+mdist
    while i<length(points)
        d=dist(points[i],points[i-1])
        while l+d <= lnext && i<length(points)
            i=i+1
            l=l+d
            d=dist(points[i],points[i-1])
        end

        while lnext <=l+d &&  length(mpoints)<nmarkers-1
            α=(lnext-l)/d
            push!(mpoints,Point2f(α*points[i]+ (1-α)*points[i-1]))
            lnext=lnext+mdist
        end
    end
    push!(mpoints,points[end])
end

function makeplanes(mmin,mmax,n)
    if isa(n,Number)
        if n==0
            return [Inf]
        end
        p=collect(range(mmin,mmax,length=ceil(n)+2))
        p[2:end-1]
    else
        n
    end
end

function makeplanes(xyzmin,xyzmax,x,y,z)
    planes=Vector{Vector{Float64}}(undef,0)
#    ε=1.0e-1*(xyzmax.-xyzmin)

    X=makeplanes(xyzmin[1],xyzmax[1],x)
    Y=makeplanes(xyzmin[2],xyzmax[2],y)
    Z=makeplanes(xyzmin[3],xyzmax[3],z)

    for i=1:length(X)
        x=X[i]
        x>xyzmin[1] && x<xyzmax[1]  && push!(planes,[1,0,0,-x])
    end

    for i=1:length(Y)
        y=Y[i]
        y>xyzmin[2] && y<xyzmax[2]  && push!(planes,[0,1,0,-y])
    end

    for i=1:length(Z)
        z=Z[i]
        z>xyzmin[3] && z<xyzmax[3]  && push!(planes,[0,0,1,-z])
    end
    planes
end

# Calculate isolevel values and function limits
function isolevels(ctx,func)

    limits=ctx[:limits]
    if limits==:auto || limits[1]>limits[2]
        limits=extrema(func)
    end

    if isa(ctx[:levels],Number)
        levels=collect(LinRange(limits[1],limits[2],ctx[:levels]+2))
    else
        levels=ctx[:levels]
    end

    levels,limits
end

"""
$(TYPEDSIGNATURES)

Extract values of given vector field (either nodal values of a piecewise linear vector field or a callback function
providing evaluation of the vector field for given generalized barycentric coordinates).
at all sampling points on  `offset+ i*spacing` for i in Z^d  defined by the tuples offset and spacing.

By default, offset is at the minimum of grid coordinates, and spacing is defined
the largest grid extend divided by 10.

The intermediate `rasterflux` in future versions can be used to calculate
streamlines.

The code is 3D ready.
"""
function vectorsample(grid::ExtendableGrid{Tv,Ti},v; offset=:default, spacing=:default,reltol=1.0e-10) where {Tv,Ti}

    coord=grid[Coordinates]
    cn=grid[CellNodes]
    ncells::Int=num_cells(grid)

    dim::Int=dim_space(grid)

    eltype= dim ==2 ? Triangle2D : Tetrahedron3D

    L2G = ExtendableGrids.L2GTransformer(eltype, grid, ON_CELLS)

    # memory for  inverse of local transformation matrix
    invA=zeros(dim+1,dim+1)

    # barycentric coordinates
    λ=zeros(dim+1)

    # coordinate window
    cminmax=extrema(coord, dims=(2,))

    if offset==:default
        offset=[cminmax[i][1] for i=1:dim]
    end

    # extent  of domain
    extent=maximum([cminmax[i][2]-cminmax[i][1] for i=1:dim])

    # scale tolerance
    tol=reltol*extent

    # point spacing
    if spacing==:default
        spacing=[extent/15 for i=1:dim]
    elseif isa(spacing,Number)
        spacing=[spacing for i=1:dim]
    end
    # else assume spacing vector has been given

    # index range
    ijkmax=ones(Int,3)
    for idim=1:dim
        ijkmax[idim]=ceil(Int64,(cminmax[idim][2]-offset[idim])/spacing[idim])+1
    end

    # The ijk raster corresponds to a  tensorproduct grid
    # spanned by x,y and z coordinate vectors. Here, we build them
    # in order to avoid to calculate them from the raster indices
    rastercoord=[zeros(Float32,ijkmax[idim]) for idim=1:dim]
    for idim=1:dim
        rastercoord[idim]=collect(range(offset[idim],step=spacing[idim],length=ijkmax[idim]))
    end

    # Memory for flux vectors on ijk grid
    rasterflux=zeros(Float32,dim,ijkmax[1], ijkmax[2], ijkmax[3])

    # type stable versions of offset and spacing
    O=zeros(dim)
    O.=offset
    S=zeros(dim)
    S.=spacing

    # memory for point, vector to be investigated
    X=zeros(Float32,dim)
    V=zeros(Float32,dim)

    tmin=zeros(Float64,dim)
    tmax=zeros(Float64,dim)

    # index vector
    I=ones(Int,3)

    # cell extent
    Imin=ones(Int,3)
    Imax=ones(Int,3)

    for icell::Int=1:ncells
        update_trafo!(L2G, icell) # 1 alloc: the only one left in this cell loop

        # Cell coordinate window
        inode=cn[1,icell]
        @views  tmin.=coord[:,inode]
        @views  tmax.=coord[:,inode]
        for i=2:dim+1
            inode=cn[i,icell]
            for idim=1:dim
                tmin[idim]=min(tmin[idim],coord[idim,inode])
                tmax[idim]=max(tmax[idim],coord[idim,inode])
            end
        end

        # min and max of raster indices falling into cell coordinate window
        for idim=1:dim
            Imin[idim]=floor(Int64,(tmin[idim]-O[idim])/S[idim])+1
            Imax[idim]=ceil(Int64,(tmax[idim]-O[idim])/S[idim])+1
        end

        # For raster indices falling into cell coordinate window,
        # check if raster points are in the cell
        # If so, obtain P1 interpolated raster data and
        # assign them to rasterflux
        for I[1] ∈ Imin[1]:Imax[1] # 0 alloc
            for I[2] ∈ Imin[2]:Imax[2]
                for I[3] ∈ Imin[3]:Imax[3]

                    # Fill raster point to be tested
                    for idim=1:dim
                        X[idim]=rastercoord[idim][I[idim]]
                    end

                    # Get barycentric coordinates
                    bary!(λ,invA,L2G,X)

                    # Check positivity of bc coordinates with some slack for
                    # round-off errors. Therefore a point may be found in two
                    # neigboring triangles. Constraining points to the raster ensures
                    # that only the last of them is taken.
                    if all(x-> x>-tol,λ)
                        # Interpolate vector value
                        if typeof(v) <: Function
                            v(V,λ,icell)
                        else
                            fill!(V,0.0)
                            for inode=1:dim+1
                                for idim=1:dim
                                    V[idim]+=λ[inode]*v[idim,cn[inode,icell]]
                                end
                            end
                        end

                        for idim=1:dim
                            rasterflux[idim,I[1],I[2],I[3]]=V[idim]
                        end
                    end
                end
            end
        end
    end
    rastercoord, rasterflux
end

"""
$(TYPEDSIGNATURES)

Extract  nonzero fluxes for quiver plots from rastergrid.

Returns qc, qv -  `d x nquiver` matrices.

If vnormalize is true, the vector field is normalized to vscale*min(spacing), otherwise, it
is scaled by vscale
Result data are meant to  be ready for being passed to calls to `quiver`.

"""
function quiverdata(rastercoord, rasterflux;   vscale=1.0, vnormalize=true)
    dim=length(rastercoord)

    imax=length(rastercoord[1])
    jmax=length(rastercoord[2])
    spacing=(rastercoord[1][2]-rastercoord[1][1],rastercoord[2][2]-rastercoord[2][1])
    kmax=1
    if dim>2
        spacing=(spacing...,rastercoord[3][2]-rastercoord[3][1])
        kmax=length(rastercoord[3])
    end

    # memory for point, vector to be investigated
    X=zeros(Float32,dim)
    V=zeros(Float32,dim)

    qvcoord=Vector{SVector{dim,Float32}}(undef,0)
    qvvalues=Vector{SVector{dim,Float32}}(undef,0)

    I=ones(Int,3)
    for I[1]=1:imax # 0 allocs besides of push
        for I[2]=1:jmax
            for I[3]=1:kmax
                for idim=1:dim
                    V[idim]=rasterflux[idim,I[1],I[2],I[3]]
                    X[idim]=rastercoord[idim][I[idim]]
                end
                if !iszero(V)
                    push!(qvcoord,X)
                    push!(qvvalues,V)
                end
            end
        end
    end

    # Reshape into matrices
    qc=reshape(reinterpret(Float32,qvcoord),(2,length(qvcoord)))
    qv=reshape(reinterpret(Float32,qvvalues),(2,length(qvvalues)))

    # Normalize vectors to raster point spacing
    if vnormalize
        @views vmax=maximum(norm,(qv[:,i] for i=1:length(qvvalues)))
        vscale=vscale*min(spacing...)/vmax
    end

    # Scale vectors with user input
    qv.*=vscale
    qc,qv
end

function bary!(λ,invA,L2G,x)
    mapderiv!(invA,L2G,nothing)
    fill!(λ ,0)
    for j = 1 : length(x)
        dj=x[j] - L2G.b[j]
        for  k = 1 : length(x)
            λ[k] += invA[j,k] * dj
        end
    end
    ExtendableGrids.postprocess_xreftest!(λ,Triangle2D)
end


text:
module ColonizationTests
    using MetacommunityDynamics
    using Test

    # ------------------------------------------------------------------------------------
    #
    # RandomColonization constructors
    #
    # ------------------------------------------------------------------------------------

    # ------------------------------------------------------------------------------------
    # RandomColonization constructor with no arguments
    # ------------------------------------------------------------------------------------
    rc_no_kw = RandomColonization()
    @test typeof(rc_no_kw) <: Colonization

    # ------------------------------------------------------------------------------------
    # RandomColonization constructor with keyword arguments
    # ------------------------------------------------------------------------------------
    v = rand()
    rc_kw = RandomColonization(probability=v)
    @test typeof(rc_kw) <: Colonization
    @test rc_kw.probability == v

    # ------------------------------------------------------------------------------------
    #
    # LevinsColonization constructors
    #
    # ------------------------------------------------------------------------------------

    # ------------------------------------------------------------------------------------
    # LevinsColonization constructor with no arguments
    # ------------------------------------------------------------------------------------
    rc_kw = LevinsColonization()
    @test typeof(rc_kw) <: Colonization

    # ------------------------------------------------------------------------------------
    # LevinsColonization constructor with no arguments
    # ------------------------------------------------------------------------------------
    v = rand()
    lc_kw = LevinsColonization(probability=v)
    @test typeof(lc_kw) <: Colonization
    @test lc_kw.probability == v

    # ------------------------------------------------------------------------------------
    #
    # SpatiallyExplicitLevinsColonization constructors
    #
    # ------------------------------------------------------------------------------------

end


text:
struct Liouvillian{N,F,D} <: AbstractParameterizedFunction{true}
    L₀::SparseMatrixCSC{ComplexF64,Int}
    Lₙ::NTuple{N,SparseMatrixCSC{ComplexF64,Int}}
    fₙ::F
    pₙ::NTuple{N,Vector{Float64}}
    dims::Dims{D}
    mass_matrix::UniformScaling{Bool}
    analytic::Nothing
end

Liouvillian(dims::Dims{D},L₀,Lₙ=(),fₙ=(),pₙ=()) where {D} =
    Liouvillian{length(Lₙ),typeof(fₙ),D}(L₀,Lₙ,fₙ,pₙ,dims,UniformScaling{Bool}(true),nothing)

dims(L::Liouvillian) = L.dims

dimsmatch(L::Liouvillian,A::QuObject) = dims(L)==dims(A) || throw(DimensionMismatch("subspace dimensions must match"))

# Liouvillian
function (L::Liouvillian)(dψ,ψ,p,t)
    mul!(dψ, L.L₀, ψ, 1, 0)
    applyfun!(dψ, L.Lₙ, L.fₙ, L.pₙ, t, ψ)
end
@inline function applyfun!(dψ,Lₙ,fₙ,pₙ,t,ψ)
    mul!(dψ, first(Lₙ), ψ, first(fₙ)(t,first(pₙ)), 1)
    applyfun!(dψ, tail(Lₙ), tail(fₙ), tail(pₙ), t, ψ)
end
@inline applyfun!(dψ,Lₙ::Tuple{},fₙ,pₙ,t,ψ) = nothing

# Jacobian
function (L::Liouvillian)(::Type{Val{:jac}},J,ψ,p,t)
    copy!(J, L.L₀)
    applyjac!(J, L.Lₙ, L.fₙ, L.pₙ, t)
end
@inline function applyjac!(J,Lₙ,fₙ,pₙ,t)
    J .+= first(fₙ)(t,first(pₙ)) .* first(Lₙ)
    applyjac!(J, tail(Lₙ), tail(fₙ), tail(pₙ), t)
end
@inline applyjac!(J,Lₙ::Tuple{},fₙ,pₙ,t) = nothing


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiRbacV1Subject <: SwaggerModel
    apiGroup::Any # spec type: Union{ Nothing, String } # spec name: apiGroup
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    namespace::Any # spec type: Union{ Nothing, String } # spec name: namespace

    function IoK8sApiRbacV1Subject(;apiGroup=nothing, kind=nothing, name=nothing, namespace=nothing)
        o = new()
        validate_property(IoK8sApiRbacV1Subject, Symbol("apiGroup"), apiGroup)
        setfield!(o, Symbol("apiGroup"), apiGroup)
        validate_property(IoK8sApiRbacV1Subject, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiRbacV1Subject, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiRbacV1Subject, Symbol("namespace"), namespace)
        setfield!(o, Symbol("namespace"), namespace)
        o
    end
end # type IoK8sApiRbacV1Subject

const _property_map_IoK8sApiRbacV1Subject = Dict{Symbol,Symbol}(Symbol("apiGroup")=>Symbol("apiGroup"), Symbol("kind")=>Symbol("kind"), Symbol("name")=>Symbol("name"), Symbol("namespace")=>Symbol("namespace"))
const _property_types_IoK8sApiRbacV1Subject = Dict{Symbol,String}(Symbol("apiGroup")=>"String", Symbol("kind")=>"String", Symbol("name")=>"String", Symbol("namespace")=>"String")
Base.propertynames(::Type{ IoK8sApiRbacV1Subject }) = collect(keys(_property_map_IoK8sApiRbacV1Subject))
Swagger.property_type(::Type{ IoK8sApiRbacV1Subject }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiRbacV1Subject[name]))}
Swagger.field_name(::Type{ IoK8sApiRbacV1Subject }, property_name::Symbol) =  _property_map_IoK8sApiRbacV1Subject[property_name]

function check_required(o::IoK8sApiRbacV1Subject)
    (getproperty(o, Symbol("kind")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiRbacV1Subject }, name::Symbol, val)
end


text:
module TestConfig

using ..Fixtures
using FactCheck
using Log4jl
import Log4jl: Configuration, name, source, logger, logger!, loggers, default!, root,
               appender, appenders, appender!, state, configure,
               LoggerConfig, level, level!, isadditive, references, reference!

testapnd = Fixtures.TestAppender()

facts("Configurations") do
    context("should have required methods implemented") do
        cfg = Fixtures.InvalidConfiguration()
        @fact_throws AssertionError name(cfg)
        @fact_throws AssertionError source(cfg)
        @fact_throws AssertionError root(cfg)
        @fact_throws AssertionError logger(cfg, Fixtures.TESTNAME)
        @fact_throws AssertionError loggers(cfg)
        @fact_throws AssertionError appender(cfg, Fixtures.TESTNAME)
        @fact_throws AssertionError appenders(cfg)
        @fact appender!(cfg, Fixtures.TESTNAME, testapnd) --> nothing
    end

    context("there exists empty configuration") do
        cfg = Log4jl.NullConfiguration()
        @fact name(cfg) --> "Null"
        @fact source(cfg) --> isempty
        @fact isa(root(cfg), Log4jl.LoggerConfig) --> true
        @fact logger(cfg, Fixtures.TESTNAME) --> root(cfg)
        @fact loggers(cfg) --> isempty
        @fact appender(cfg, Fixtures.TESTNAME) --> nothing
        @fact appenders(cfg) --> isempty
        @fact appender!(cfg, Fixtures.TESTNAME, testapnd) --> nothing
        @fact logger!(cfg, Fixtures.TESTNAME, LoggerConfig()) --> nothing
        @fact state(cfg) --> Log4jl.LifeCycle.STOPPED
        @fact configure(cfg) --> nothing
    end

    context("should be referenced to logger with `LoggerConfig` object") do
        root = Log4jl.LoggerConfig()
        @fact name(root) --> isempty
        @fact level(root) --> Log4jl.LOG4JL_DEFAULT_STATUS_LEVEL
        @fact level(Nullable(root)) --> Log4jl.LOG4JL_DEFAULT_STATUS_LEVEL
        level!(root, Fixtures.TESTLEVEL)
        @fact level(root) --> Fixtures.TESTLEVEL
        @fact level(Nullable(root)) --> Fixtures.TESTLEVEL
        @fact root --> not(isadditive)
        context("which references appenders") do
            @fact length(references(root)) --> 0
            reference!(root, testapnd)
            @fact length(references(root)) --> greater_than(0)
        end
    end

    context("there exists default configuration") do
        cfg = Log4jl.DefaultConfiguration()
        @fact name(cfg) --> "Default"
        @fact source(cfg) --> isempty
        @fact isa(root(cfg), Log4jl.LoggerConfig) --> true
        context("with default parameters when not started") do
            @fact logger(cfg, Fixtures.TESTNAME) --> root(cfg)
            @fact loggers(cfg) --> isempty
            @fact appenders(cfg) --> isempty
            @fact state(cfg) --> Log4jl.LifeCycle.INITIALIZED
        end
        context("it should be initialized by `configure`") do
            # No appender and references
            @fact length(appenders(cfg)) --> 0
            @fact length(loggers(cfg)) --> 0
            @fact length(references(root(cfg))) --> 0
            configure(cfg)
            apnds = appenders(cfg) |> values
            @fact length(apnds) --> greater_than(0)
            apnd = appender(cfg, "Default")
            @fact isa(apnd, Log4jl.Appenders.Console) --> true
            apndrefs = references(root(cfg))
            @fact length(apndrefs) --> greater_than(0)
            @fact isa(first(apndrefs), Log4jl.Appenders.Reference) --> true
            @fact first(apnds) --> first(apndrefs).appender

        end
    end
    context("have a root logger which can be setup by calling `default!` function") do
        cfg = Log4jl.DefaultConfiguration()
        @fact length(appenders(cfg)) --> 0
        @fact length(loggers(cfg)) --> 0
        apndref = default!(cfg)
        @fact length(appenders(cfg)) --> greater_than(0)
        @fact isa(apndref, Log4jl.Appenders.Reference) --> true
        @fact appenders(cfg)["Default"] --> apndref.appender
    end
end

end


text:
A = zeros(3, 4, 5);
A[:] = 1:60
Ap = @view(PermutedDimsArray(A,(3,1,2))[:,1:2,1])';

ap_index = ArrayInterface.StrideIndex(Ap)
for x_i in axes(Ap, 1)
    for y_i in axes(Ap, 2)
        @test ap_index[x_i, y_i] == ap_index[x_i, y_i]
    end
end
@test @inferred(ArrayInterface.known_offsets(ap_index)) === ArrayInterface.known_offsets(Ap)
@test @inferred(ArrayInterface.known_offset1(ap_index)) === ArrayInterface.known_offset1(Ap)
@test @inferred(ArrayInterface.offsets(ap_index, 1)) === ArrayInterface.offset1(Ap)
@test @inferred(ArrayInterface.offsets(ap_index, static(1))) === ArrayInterface.offset1(Ap)
@test @inferred(ArrayInterface.known_strides(ap_index)) === ArrayInterface.known_strides(Ap)
@test @inferred(ArrayInterface.contiguous_axis(ap_index)) == 1
@test @inferred(ArrayInterface.contiguous_axis(ArrayInterface.StrideIndex{2,(1,2),nothing,NTuple{2,Int},NTuple{2,Int}})) == nothing
@test @inferred(ArrayInterface.stride_rank(ap_index)) == (1, 3)

let v = Float64.(1:10)', v2 = transpose(parent(v))
  sv = @view(v[1:5])'
  sv2 = @view(v2[1:5])'
  @test @inferred(ArrayInterface.StrideIndex(sv)) === @inferred(ArrayInterface.StrideIndex(sv2)) === ArrayInterface.StrideIndex{2, (2, 1), 2}((StaticInt(1), StaticInt(1)), (StaticInt(1), StaticInt(1)))
  @test @inferred(ArrayInterface.stride_rank(parent(sv))) === @inferred(ArrayInterface.stride_rank(parent(sv2))) === (StaticInt(1),)
end