text:
using Base.Ryu

abstract type FmtType end

struct FmtDec <: FmtType end
struct FmtOct <: FmtType end
struct FmtHex <: FmtType end
struct FmtChr <: FmtType end
struct FmtStr <: FmtType end
struct FmtPtr <: FmtType end

struct FmtFltE <: FmtType end
struct FmtFltF <: FmtType end
struct FmtFltG <: FmtType end
struct FmtFltA <: FmtType end

struct FmtN <: FmtType end

const VALID_FMTS = b"duoxefgacsipnDUOXEFGACS"

const FMT_TYPES = [FmtDec, FmtDec, FmtOct, FmtHex, FmtFltE, FmtFltF, FmtFltG, FmtFltA, FmtChr, FmtStr,
                   FmtDec, FmtPtr, FmtN]

# format specifier categories
const FmtInts = Union{FmtDec, FmtOct, FmtHex}
const FmtFlts = Union{FmtFltE, FmtFltF, FmtFltG, FmtFltA}

"""
Typed representation of a format specifier.

Fields are the various modifiers allowed for various format specifiers.
"""
struct FmtSpec{T<:FmtType}
    char::UInt8
    tsep::UInt8 # thousands separator (default ',')
    leftalign::Bool
    plus::Bool
    space::Bool
    zero::Bool
    altf::Bool # alternate formatting ('#' flag)
    width::Int # negative means *, 0 means no width specified
    prec::Int # -2 means *, -1 means no precision specified
end

# recreate the format specifier string from a typed FmtSpec
Base.string(f::FmtSpec; modifier::String="") =
    string("%",
           f.leftalign ? "-" : "",
           f.plus ? "+" : "",
           f.space ? " " : "",
           f.zero ? "0" : "",
           f.altf ? "#" : "",
           f.tsep == 0 ? "" : "'",
           f.width > 0 ? "$(f.width)" : f.width < 0 ? "*" : "",
           f.prec < -1 ? "*" : f.prec >= 0 ? ".$(f.prec)" : "",
           modifier,
           Char(f.char))
Base.show(io::IO, f::FmtSpec) = print(io, string(f))

floatfmt(s::FmtSpec) =
    FmtSpec{FmtFltF}(s.char, s.tsep, s.leftalign, s.plus, s.space, s.zero, s.altf, s.width, 0)
ptrfmt(s::FmtSpec, x) =
    FmtSpec{FmtHex}(s.char, s.tsep, s.leftalign, s.plus, s.space, s.zero, true, s.width,
                    ifelse(sizeof(x) == 8, 16, 8))

function FmtSpec(s::FmtSpec{F}; width=s.width, prec=s.prec) where {F}
    s.width < 0 || throw(ArgumentError("Already has width $(s.width)"))
    s.prec < -1 || throw(ArgumentError("Already has precision $(s.prec)"))
    width < 0 && throw(ArgumentError("Width $(width) must be >= 0"))
    prec < -1 && throw(ArgumentError("Precision $(prec) must be >= -1"))
    FmtSpec{F}(s.char, s.tsep, s.leftalign, s.plus, s.space, s.zero, s.altf, width, prec)
end

# parse format string
function FmtSpec(f::AbstractString)
    isempty(f) && throw(ArgumentError("empty format string"))
    bytes = codeunits(f) # Note: codeunits are not necessarily *bytes*!
    len = length(bytes)
    bytes[1] === UInt8('%') || throw(ArgumentError("Format string doesn't start with %"))
    b = bytes[2]
    pos = 3
    # positioned at start of first format str %
    # parse flags
    leftalign = plus = space = zero = altf = false
    tsep = UInt8(0)
    while true
        if b == UInt8('-')
            leftalign = true
        elseif b == UInt8('+')
            plus = true
        elseif b == UInt8(' ')
            space = true
        elseif b == UInt8('0')
            zero = true
        elseif b == UInt8('#')
            altf = true
        elseif b == UInt8('\'')
            tsep = UInt8(',')
        else
            break
        end
        pos > len && throw(ArgumentError("incomplete format string: '$f'"))
        b = bytes[pos]
        pos += 1
    end
    leftalign && (zero = false)
    # parse width
    if b == UInt('*')
        width = -1
    else
        width = 0
        while b - UInt8('0') < 0x0a
            width = 10 * width + (b - UInt8('0'))
            b = bytes[pos]
            pos += 1
            pos > len && break
        end
    end
    # parse prec
    prec = -1
    parsedprecdigits = false
    if b == UInt8('.')
        pos > len && throw(ArgumentError("incomplete format string: '$f'"))
        parsedprecdigits = true
        b = bytes[pos]
        pos += 1
        if pos <= len
            if b == UInt('*')
                prec == -2
            else
                prec = 0
                while b - UInt8('0') < 0x0a
                    prec = 10 * prec + (b - UInt8('0'))
                    b = bytes[pos]
                    pos += 1
                    pos > len && break
                end
            end
        end
    end
    # parse length modifier (ignored)
    if b == UInt8('h') || b == UInt8('l')
        prev = b
        b = bytes[pos]
        pos += 1
        if b == prev
            pos > len && throw(ArgumentError("invalid format string: '$f'"))
            b = bytes[pos]
            pos += 1
        end
    elseif b in b"Ljqtz"
        b = bytes[pos]
        pos += 1
    end
    # parse type
    fmtind = findfirst(isequal(b), VALID_FMTS)
    fmtind === nothing &&
        throw(ArgumentError(string("invalid format string: '$f', ",
                                   "invalid type specifier: '$(Char(b))'")))
    # Check for uppercase variants
    fmtind > 13 && (fmtind -= 13)
    fmttyp = FMT_TYPES[fmtind]
    if fmttyp <: FmtInts && prec > 0
        zero = false
    elseif !parsedprecdigits
        if (fmttyp === FmtStr || fmttyp === FmtChr || fmttyp === FmtFltA)
            prec = -1
        elseif fmttyp <: FmtFlts
            prec = 6
        end
    end
    FmtSpec{fmttyp}(b, tsep, leftalign, plus, space, zero, altf, width, prec)
end

@inline gethexbase(spec) = spec.char < UInt8('a') ? b"0123456789ABCDEF" : b"0123456789abcdef"

@inline upchar(spec, ch) = (spec.char & 0x20) | UInt8(ch)

# write out a single arg according to format options
# char
@inline function writechar(buf, pos, c::Char)
    u = bswap(reinterpret(UInt32, c))
    while true
        @inbounds buf[pos] = u % UInt8
        pos += 1
        (u >>= 8) == 0 && break
    end
    return pos
end

@inline function padn(buf, pos, cnt)
    @inbounds for _ = 1:cnt
        buf[pos] = UInt8(' ')
        pos += 1
    end
    pos
end

@inline function padzero(buf, pos, cnt)
    @inbounds for _ = 1:cnt
        buf[pos] = UInt8('0')
        pos += 1
    end
    pos
end

function _fmt(buf, pos, spec::FmtSpec{FmtChr}, arg)
    ch = arg isa String ? arg[1] : Char(arg)
    width = spec.width - 1
    width <= 0 && return writechar(buf, pos, ch)
    if spec.leftalign
        padn(buf, writechar(buf, pos, ch), width)
    else
        writechar(buf, padn(buf, pos, width), ch)
    end
end

@inline function outch(buf, pos, ch)
    @inbounds buf[pos] = UInt8(ch)
    pos + 1
end
@inline function outch(buf, pos, c1, c2)
    @inbounds buf[pos] = UInt8(c1)
    @inbounds buf[pos + 1] = UInt8(c2)
    pos + 2
end
@inline function outch(buf, pos, c1, c2, c3)
    @inbounds buf[pos] = UInt8(c1)
    @inbounds buf[pos + 1] = UInt8(c2)
    @inbounds buf[pos + 2] = UInt8(c3)
    pos + 3
end

# strings
function _fmt(buf, pos, spec::FmtSpec{FmtStr}, arg)
    altf, width, prec = spec.altf, spec.width, spec.prec
    str = altf && (arg isa Symbol || arg isa AbstractString) ? repr(arg) : string(arg)
    slen = length(str)
    op = p = prec == -1 ? slen : min(slen, prec)
    # Make sure there is enough room in buffer
    nlen = max(width, p)
    buflen = sizeof(buf) - pos
    buflen < nlen && resize!(buf, nlen + pos)
    !spec.leftalign && width > p && (pos = padn(buf, pos, width - p))
    for c in str
        p == 0 && break
        pos = writechar(buf, pos, c)
        p -= 1
    end
    spec.leftalign && width > op && return padn(buf, pos, width - op)
    return pos
end

const BaseInt = Union{Int8, Int16, Int32, Int64, Int128}
const BaseUns = Union{UInt8, UInt16, UInt32, UInt64, UInt128}

base(::Type{FmtOct}) = 8
base(::Type{FmtDec}) = 10
base(::Type{FmtHex}) = 16

# integers

_fmt(buf, pos, spec::FmtSpec{<:FmtInts}, arg::AbstractFloat) =
    _fmt(buf, pos, floatfmt(spec), arg)

_fmt(buf, pos, spec::FmtSpec{T}, arg::Real) where {T<:FmtInts} =
    _fmt(buf, pos, spec, arg < 0, string(abs(arg); base=base(T)))
_fmt(buf, pos, spec::FmtSpec{<:FmtInts}, arg::BaseUns) =
    _fmt(buf, pos, spec, false, arg)
_fmt(buf, pos, spec::FmtSpec{<:FmtInts}, arg::BaseInt) =
    _fmt(buf, pos, spec, arg < 0, unsigned(abs(arg)))

function _fmt(buf, pos, spec::FmtSpec{F}, neg, x::T) where {F<:FmtInts,T<:Union{String,BaseUns}}
    n = T === String ? sizeof(x) :
        F === FmtDec ? dec_len(x) :
        F === FmtHex ? (sizeof(x)<<1) - (leading_zeros(x)>>2) :
        div((sizeof(x)<<3) - leading_zeros(x)+2, 3)
    i = n
    arglen = n + (neg || (spec.plus | spec.space)) +
        ((spec.altf && (F !== FmtDec)) ? ifelse(F === FmtOct, 1, 2) : 0)
    width, prec = spec.width, spec.prec
    arglen2 = arglen < width && prec > 0 ? arglen + min(max(0, prec - n), width - arglen) : arglen

    # Make sure that remaining output buffer is large enough
    # This means that it isn't necessary to preallocate for cases that usually will never happen
    buflen = sizeof(buf) - pos
    buflen < arglen2 && resize!(buf, arglen2 + pos)

    !spec.leftalign && !spec.zero && arglen2 < width &&
        (pos = padn(buf, pos, width - arglen2))
    if neg
        pos = outch(buf, pos, '-')
    elseif spec.plus # plus overrides space
        pos = outch(buf, pos, '+')
    elseif spec.space
        pos = outch(buf, pos, ' ')
    end
    if spec.altf
        if F === FmtOct
            pos = outch(buf, pos, '0')
        elseif F === FmtHex
            pos = outch(buf, pos, '0', upchar(spec, 'X'))
        end
    end
    if spec.zero && arglen2 < width
        pos = padzero(buf, pos, width - arglen2)
    elseif n < prec
        pos = padzero(buf, pos, prec - n)
    elseif arglen < arglen2
        pos = padzero(buf, pos, arglen2 - arglen)
    end
    if T === String
        GC.@preserve buf x begin
            unsafe_copyto!(pointer(buf, pos), pointer(x), n)
        end
    elseif F === FmtDec
        while i > 0
            d, r = divrem(x, 10)
            @inbounds buf[pos + i - 1] = UInt8('0') + r
            x = oftype(x, d)
            i -= 1
        end
    elseif F === FmtHex
        hexbase = gethexbase(spec)
        while i > 0
            @inbounds buf[pos + i - 1] = hexbase[(x & 0x0f) + 1]
            x >>= 4
            i -= 1
        end
    else # F === FmtOct
        while i > 0
            @inbounds buf[pos + i - 1] = UInt8('0') + (x & 0x07)
            x >>= 3
            i -= 1
        end
    end
    (spec.leftalign && arglen2 < width) ? padn(buf, pos + n, width - arglen2) : (pos + n)
end

# floats
"""
    tofloat(x)

Convert an argument to a Base float type for printf formatting.
By default, arguments are converted to `Float64` via `Float64(x)`.
Custom numeric types that have a conversion to a Base float type
that wish to hook into printf formatting can extend this method like:

```julia
Printf.tofloat(x::MyCustomType) = convert_my_custom_type_to_float(x)
```

For arbitrary precision numerics, you might extend the method like:

```julia
Printf.tofloat(x::MyArbitraryPrecisionType) = BigFloat(x)
```
"""
tofloat(x) = Float64(x)
tofloat(x::Base.IEEEFloat) = x

function output_fmt_a(buf, pos, spec, neg, x)
    if neg
        pos = outch(buf, pos, '-')
    elseif spec.plus
        pos = outch(buf, pos, '+')
    elseif spec.space
        pos = outch(buf, pos, ' ')
    end
    isnan(x) && return outch(buf, pos, 'N', 'a', 'N')
    isfinite(x) || return outch(buf, pos, 'I', 'n', 'f')
    prec = spec.prec
    pos = outch(buf, pos, '0', upchar(spec, 'X'))
    if x == 0
        pos = outch(buf, pos, '0')
        prec > 0 && (pos = padzero(buf, pos, prec))
        return outch(buf, pos, upchar(spec, 'P'), '+', '0')
    end
    s, p = frexp(x)
    prec = spec.prec
    if prec > -1
        sigbits = 4 * min(prec, 13)
        s = 0.25 * round(ldexp(s, 1 + sigbits))
        # ensure last 2 exponent bits either 01 or 10
        u = (reinterpret(UInt64, s) & 0x003f_ffff_ffff_ffff) >> (52 - sigbits)
        i = n = (sizeof(u) << 1) - (leading_zeros(u) >> 2)
    else
        s *= 2.0
        u = (reinterpret(UInt64, s) & 0x001f_ffff_ffff_ffff)
        t = (trailing_zeros(u) >> 2)
        u >>= (t << 2)
        i = n = 14 - t
    end
    frac = u > 9 || spec.altf || prec > 0
    hexbase = gethexbase(spec)
    while i > 1
        buf[pos + i] = hexbase[(u & 0x0f) + 1]
        u >>= 4
        i -= 1
        prec -= 1
    end
    frac && (buf[pos + 1] = UInt8('.'))
    buf[pos] = hexbase[(u & 0x0f) + 1]
    pos += n + frac
    while prec > 0
        pos = outch(buf, pos, '0')
        prec -= 1
    end
    pos = outch(buf, pos, upchar(spec, 'P'))
    p -= 1
    pos = outch(buf, pos, p < 0 ? UInt8('-') : UInt8('+'))
    p = p < 0 ? -p : p
    n = i = ndigits(p, base=10, pad=1)
    while i > 0
        d, r = divrem(p, 10)
        buf[pos + i - 1] = UInt8('0') + r
        p = oftype(p, d)
        i -= 1
    end
    return pos + n
end

const _strspec = Dict{FmtSpec,ASCIIStr}()

"""Get the format specification, prepared for use with MPFR for BigFloats"""
function _get_strspec(spec)
    global _strspec
    chkspec = get(_strspec, spec, nothing)
    chkspec === nothing || return chkspec
    _strspec[spec] = string(spec; modifier="R")
end

@static if VERSION < v"1.5"
_snprintf(ptr, siz, spec, arg) =
    ccall((:mpfr_snprintf,:libmpfr), Int32,
          (Ptr{UInt8}, Culong, Ptr{UInt8}, Ref{BigFloat}...),
          ptr, siz, spec, arg)
else
_snprintf(ptr, siz, spec, arg) =
    @ccall "libmpfr".mpfr_snprintf(ptr::Ptr{UInt8}, siz::Csize_t, spec::Ptr{UInt8};
                                   arg::Ref{BigFloat})::Cint
end

function _fmt(buf, pos, spec::FmtSpec{<:FmtFlts}, arg::BigFloat)
    isfinite(arg) || return _fmt(buf, pos, spec, Float64(arg))
    strspec = _get_strspec(spec)
    siz = length(buf) - pos + 1
    # Calculate size needed for most common outputs
    len = max(spec.width, ceil(Int, precision(arg) * log(2)/log(10)) + 24)
    if len > siz
        resize!(buf, pos + len + 1)
        siz = length(buf) - pos + 1
    end
    GC.@preserve buf begin
        len = _snprintf(pointer(buf, pos), siz, strspec, arg)
        if len > siz
            resize!(buf, pos + len + 1)
            len = _snprintf(pointer(buf, pos), len + 1, strspec, arg)
        end
    end
    len > 0 || error("invalid printf formatting for BigFloat")
    pos + len
end

function _fmt(buf, pos, spec::FmtSpec{T}, arg) where {T <: FmtFlts}
    # Make sure there is enough room
    width = spec.width
    buflen = sizeof(buf) - pos
    needed = max(width, 309 + 17 + 5)
    buflen < needed && resize!(buf, pos + needed)

    x = tofloat(arg)
    if T === FmtFltE
        newpos = Ryu.writeexp(buf, pos, x, spec.prec, spec.plus, spec.space,
                              spec.altf, upchar(spec, 'E'), UInt8('.'))
    elseif T === FmtFltF
        newpos = Ryu.writefixed(buf, pos, x, spec.prec, spec.plus, spec.space, spec.altf,
                                UInt8('.'))
    elseif T === FmtFltG
        prec = spec.prec
        prec = prec == 0 ? 1 : prec
        x = round(x, sigdigits=prec)
        newpos = Ryu.writeshortest(buf, pos, x, spec.plus, spec.space, spec.altf, prec,
                                   upchar(spec, 'E'), true, UInt8('.'))
    elseif T === FmtFltA
        newpos = output_fmt_a(buf, pos, spec, x < 0, abs(x))
    end
    if newpos - pos < width
        # need to pad
        if spec.leftalign
            # easy case, just pad spaces after number
            newpos = padn(buf, newpos, width - (newpos - pos))
        else
            # right aligned
            n = width - (newpos - pos)
            if spec.zero
                ex = (arg < 0 || (spec.plus | spec.space)) + ifelse(T === FmtFltA, 2, 0)
                so = pos + ex
                len = (newpos - pos) - ex
                copyto!(buf, so + n, buf, so, len)
                for i = so:(so + n - 1)
                    buf[i] = UInt8('0')
                end
                newpos += n
            else
                copyto!(buf, pos + n, buf, pos, newpos - pos)
                for i = pos:(pos + n - 1)
                    buf[i] = UInt8(' ')
                end
                newpos += n
            end
        end
    end
    return newpos
end

# pointers
_fmt(buf, pos, spec::FmtSpec{FmtPtr}, arg) =
    _fmt(buf, pos, ptrfmt(spec, arg), UInt(arg))

@inline _dec_len1(v) = ifelse(v < 100, ifelse(v < 10, 1, 2), 3)
@inline _dec_len2(v) = v < 1_000 ? _dec_len1(v) : ifelse(v < 10_000, 4, 5)
@inline _dec_len4(v) = v < 100_000 ? _dec_len2(v) :
    (v < 10_000_000 ? ifelse(v < 1_000_000, 6, 7) : ifelse(v < 100_000_000, 8, 9))
@inline function _dec_len8(v)
    if v < 1_000_000_000 # 1 - 9 digits
        _dec_len4(v)
    elseif v < 10_000_000_000_000 # 10 - 13 digits
        if v < 100_000_000_000
            ifelse(v < 10_000_000_000, 10, 11)
        else
            ifelse(v < 1_000_000_000_000, 12, 13)
        end
    elseif v < 100_000_000_000_000_000 # 14 - 17 digits
        if v < 1_000_000_000_000_000
            ifelse(v < 100_000_000_000_000, 14, 15)
        else
            ifelse(v < 10_000_000_000_000_000, 16, 17)
        end
    else
        v < 10_000_000_000_000_000_000 ? ifelse(v < 1_000_000_000_000_000_000, 18, 19) : 20
    end
end

@inline dec_len(v::UInt8)  = _dec_len1(v)
@inline dec_len(v::UInt16) = _dec_len2(v)
@inline dec_len(v::UInt32) = _dec_len4(v)
@inline dec_len(v::UInt64) = _dec_len8(v)
@inline function dec_len(v::UInt128)
    v <= typemax(UInt64) && return _dec_len8(v%UInt64)
    v = div(v, 0x8ac7230489e80000) # 10^19
    v <= typemax(UInt64) ? _dec_len8(v%UInt64) + 19 : 39
end


text:
function get_prediction_variances(model_matrix::Array{Float64, 2})
    information_matrix = model_matrix' * model_matrix

    if det(information_matrix) != 0.0
        dispersion_matrix = inv(information_matrix)
        rows              = size(dispersion_matrix, 1)

        prediction_variances = [dispersion_matrix[i, :]' * dispersion_matrix *
                                dispersion_matrix[i, :] for i = 1:rows]

        return prediction_variances
    else
        return 0.0
    end
end

"""
    d_optimality(model_matrix::Array{Float64, 2})

Compute the ``D``-optimality of a given design's model matrix.

# Examples

```jldoctest
julia> using ExperimentalDesign, DataStructures, StatsModels

julia> A = plackett_burman(4)
4×3 Array{Int64,2}:
  1   1   1
  1  -1  -1
 -1  -1   1
 -1   1  -1

julia> factors = OrderedDict([(:f1, [-1., 1.]), (:f2, [-1., 1.]), (:f3, [-1., 1.])])
DataStructures.OrderedDict{Symbol,Array{Float64,1}} with 3 entries:
  :f1 => [-1.0, 1.0]
  :f2 => [-1.0, 1.0]
  :f3 => [-1.0, 1.0]

julia> m = generate_model_matrix(@formula(y ~ f1 + f2 + f3), float(A), factors)
4×4 Array{Float64,2}:
 1.0   1.0   1.0   1.0
 1.0   1.0  -1.0  -1.0
 1.0  -1.0  -1.0   1.0
 1.0  -1.0   1.0  -1.0

julia> d_optimality(m)
256.0

```

# Formula

The ``D``-optimality of a design is the determinant of the information matrix
``\\mathbf{X}^{\\prime}\\mathbf{X}``, where ``\\mathbf{X}`` is the model matrix of
a design.

"""
function d_optimality(model_matrix::Array{Float64, 2})
    det_information_matrix = det(model_matrix' * model_matrix)

    return det_information_matrix < 0.0 ? 0.0 : det_information_matrix
end

"""
    d_efficiency_lower_bound(model_matrix::Array{Float64, 2})

Compute a lower bound for the ``D``-efficiency of a given design's model matrix
according to Castillo's "Process Optimization : A Statistical Approach".

```jldoctest
julia> using ExperimentalDesign, DataStructures, StatsModels

julia> A = plackett_burman(4)
4×3 Array{Int64,2}:
  1   1   1
  1  -1  -1
 -1  -1   1
 -1   1  -1

julia> factors = OrderedDict([(:f1, [-1., 1.]), (:f2, [-1., 1.]), (:f3, [-1., 1.])])
DataStructures.OrderedDict{Symbol,Array{Float64,1}} with 3 entries:
  :f1 => [-1.0, 1.0]
  :f2 => [-1.0, 1.0]
  :f3 => [-1.0, 1.0]

julia> m = generate_model_matrix(@formula(y ~ f1 + f2 + f3), float(A), factors)
4×4 Array{Float64,2}:
 1.0   1.0   1.0   1.0
 1.0   1.0  -1.0  -1.0
 1.0  -1.0  -1.0   1.0
 1.0  -1.0   1.0  -1.0

julia> d_efficiency_lower_bound(m)
1.0

```

# Formula

For a design ``A_{n,p}`` with ``n`` experiments or rows, ``p`` factors or
columns, and model matrix ``\\mathbf{X}``, the lower bound for the
``D``-efficiency of ``A`` is:

```math
D_{eff}^{(L)} = \\dfrac{|\\mathbf{X}^{\\prime}\\mathbf{X}|^{1/p}}{n}
```
"""
function d_efficiency_lower_bound(model_matrix::Array{Float64, 2})
    return ^(d_optimality(model_matrix), 1 / size(model_matrix, 2)) /
           size(model_matrix, 1)
end

function a_optimality(model_matrix::Array{Float64, 2})
    information_matrix = model_matrix' * model_matrix

    if det(information_matrix) != 0.0
        return trace(inv(information_matrix)) / size(model_matrix, 2)
    else
        return 0.0
    end

end

function v_optimality(model_matrix::Array{Float64, 2})
    prediction_variances = get_prediction_variances(model_matrix)
    rows                 = size(model_matrix, 1)

    return sum(prediction_variances) / rows
end

function g_optimality(model_matrix::Array{Float64, 2})
    prediction_variances = get_prediction_variances(model_matrix)
    return max(prediction_variances...)
end

function g_efficiency(model_matrix::Array{Float64, 2})
    prediction_variances = get_prediction_variances(model_matrix)
    max_variance         = max(prediction_variances...)

    g_e = size(model_matrix, 2) / max_variance

    if g_e == Inf
        return 0.0
    else
        return g_e
    end
end

function d_efficiency_lower_bound_algdesign(model_matrix::Array{Float64, 2})
    g_e   = g_efficiency(model_matrix)
    d_elb = exp(1 - (1 / g_e))

    if d_elb == Inf
        return 0.0
    else
        return d_elb
    end
end

function condition_number(model_matrix::Array{Float64, 2})
    condition_number = cond(model_matrix)

    if condition_number == Inf
        return 0.0
    else
        return condition_number
    end
end

ext: jl


text:
# TODO: we could get away with removing samplerate from here and storing
# the A and B signals as SampleBufs
struct GolaySequence{AT<:AbstractVector, SR} <: IRMeasurement
    A::AT
    B::AT
    decay::Int
    gain::Float64
    bandlimit::Int
    samplerate::SR
end

GolaySequence(A::AT, B::AT, decay, gain, bandlimit, sr::SR) where {AT, SR} =
    GolaySequence{AT, SR}(A, B, decay, gain, bandlimit, sr)

"""
    golay(length; options...)

Create an IR measurement using a complimentary Golay sequence following the
method from [1], assuming that the system being measured has a response less
than length `samples`. The actual length of the generated stimulus sequence
might be greater than `samples`.

You can also specify the length with a duration in seconds along with a
sampling rate. The sampling rate can either be specified as a raw number or
a unitful frequency.

TODO: update length docs

## Options

$optiondoc_gain
- `bandlimit::Integer`: Defaults to 1, which creates a standard binary golay code
  with energy up to nyquist. Setting this to a higher number will create a
  bandlimited version, so `bandlimit=2` will only have energy up to 1/2 the
  nyquist frequency. Often the upper frequencies of the impulse response are
  time-variant, so they may not be accurately measurable anyways. Decreasing the
  bandwidth can also help avoid nonlinearities due to slew-rate limiting in the
  playback system.
- `decay`: Specify how long you expect the system decay length to be. The
  simulus will wait this long between the A and B parts. Note that this should
  include any delay introduced by the system or measurement process.

## Example

```julia
using Plots: plot
using MeasureIR: golay, stimulus, analyze
using Unitful: s, Hz

meas = golay(2s, 48kHz)

# generate the test stimuli suitable for playback
stim = stimulus(meas)

# create a synthetic impulse response and simulate the system. This is where
# you'd normally play the stimuli through your system and record the response
irsim = 1./exp.(0:0.1:9.9) .* (rand(100) .- 0.5)
output = conv(stim, irsim)

# analyze to reconstruct the impulse response
ir = analyze(meas, output)

plot([irsim[1:100], ir[1:100]], labels=["Convolved IR", "Measured IR"])
```

[1]: Berdahl, Edgar J. and Julius O. Smith III, "Transfer Function Measurement
Toolbox", REALSIMPLE Project. Released 2008-06-05 under the Creative Commons
License (Attribution 2.5) Center for Computer Research in Music and Acoustics
(CCRMA), Stanford University. Accessed 2018 from
https://ccrma.stanford.edu/realsimple/imp_meas/Golay_Code_Theory.html
"""
function golay end

function golay(length; decaylength=length, samplerate=1, gain=golay_gain, bandlimit=1)
    L = ceil(Int, uconvert(NoUnits, length*samplerate/bandlimit))
    dec = ceil(Int, uconvert(NoUnits, decaylength*samplerate))
    L = nextpow(2, L)
    A, B = _golay(L)
    GolaySequence(A, B, dec, gain, bandlimit, samplerate)
end

Base.:(==)(g1::GolaySequence, g2::GolaySequence) = g1.A == g2.A && g1.B == g2.B

function stimulus(sig::GolaySequence)
    stim = [sig.gain .* sig.A;
            zeros(sig.decay);
            sig.gain .* sig.B]
    if sig.bandlimit != 1
        # force upsampling ratio to be a rational (DSP.jl issue #211)
        # TODO: this is a bit broken because it doesn't handle the warm-up
        # in the convolution. but `resample` doesn't give you exactly the
        # right number of samples, so this needs to be looked into.
        # TODO: also I need to think about how we're handling the decay time
        # here, i.e. whether it shuold be before of after the resampling
        stim = filt(FIRFilter(sig.bandlimit//1), stim)
    end

    stim
end

# naive implementation, allocates a lot, but runs reasonably fast
function _golay(L)
    if L == 1
        [1.0], [1.0]
    else
        a, b = _golay(L÷2)
        [a; b], [a; -b]
    end
end

function _analyze(sig::GolaySequence, response::AbstractArray)
    # stimuli has 2 measurements, each with a length-L golay sequence. They are
    # separated by `sig.decay` samples of silence.
    L = Int(length(sig.A) * sig.bandlimit)
    # np = noisefloor(sig, response)
    # we chop off the silence part - there's nothing causal for us there
    respA = timeslice(response, 1 : L+sig.decay)
    respB = timeslice(response, L+sig.decay+1 : size(response, 1))
    # @views respB = truncresponse(
    #         timeslice(response, L+sig.decay+1 : size(response, 1)), 2L, np)
    if sig.bandlimit != 1
        # we can use the full-bandwidth zero-stuffed signal for the
        # cross-correlation, so that any weirdness caused by our upsampling
        # procedure only affects us once, and we get a nice bandlimited
        # impulse in the pure loopback case. There's an argument for using the
        # bandlimited version in that frequency content above the cutoff must
        # have been caused by nonlinearity, but for now we'll go for the clean
        # impulse. There's also some kind of strange issue where there ends up
        # being a bunch of pre-ringing garbage when using the bandlimited
        # version. Not sure what that's about. Also this way it doesn't cut off
        # the beginning of the bandlimited impulse when there's no delay.
        A = zerostuff(sig.A, Int(sig.bandlimit)) * sig.bandlimit
        B = zerostuff(sig.B, Int(sig.bandlimit)) * sig.bandlimit
        # A = filt(FIRFilter(sig.bandlimit), sig.A)
        # B = filt(FIRFilter(sig.bandlimit), sig.B)
    else
        A = sig.A
        B = sig.B
    end
    # run the cross-correlation, chopping off the non-causal part and the part
    # that would be corrupted if the IR is too long
    irA = mapslices(respA, dims=1) do v
        # warning - xcorr pads the inputs to equal length, so the total result
        # length should be 2*maxlen-1
        xc = xcorr(v, A; padmode=:longest) / (2L * sig.gain)
        time0 = (length(xc)+1)÷2-1
        xc[time0 .+ (1:sig.decay)]
    end
    irB = mapslices(respB, dims=1) do v
        xc = xcorr(v, B; padmode=:longest) / (2L * sig.gain)
        time0 = (length(xc)+1)÷2-1
        xc[time0 .+ (1:sig.decay)]
    end
    irA + irB
end

# this is no longer used because the brickwall cutoff causes unacceptable
# ringing in the time-domain signal, causing aliasing in the deconvolution
"""
upsamples/interpolates `A` treating it as periodic.
"""
function circ_upsample(A, ratio)
    L = length(A)
    L2 = L*ratio
    zpad = L2÷2 - L÷2
    spec = [rfft(A); zeros(Complex128, zpad)]
    if iseven(L)
        # in the original spectrum the nyquist bin was doubled up because of
        # periodicity. Now the energy is split between two bins.
        spec[L÷2+1] *= 0.5
    end

    irfft(spec, L2) * ratio
end

"""
upsample by zero-stuffing in between samples
"""
function zerostuff(A, ratio)
    A2 = zeros(length(A)*ratio)
    A2[1:ratio:end] .= A
    A2
end


text:
#=
[OpenBCI.jl]
Julia = 0.7+
Author = "William Herrera"
Module Version = 0.025
Copyright = "Copyright William Herrera, 2018"
Created = "16 Jan 2018"
Purpose = "EEG WiFi routines using OpenBCI Arduino hardware"
=#

module OpenBCI

export makeganglionbdfplus, makecyton8bdfplus, makecyton16bdfplus

using Dates
using EDFPlus
using HTTP
import HTTP: get, post
using JSON
using Memento
using Sockets

const logger = getlogger(@__MODULE__)

"""
See [http://docs.openbci.com/Hardware/03-Cyton_Data_Format#cyton-data-format-binary-format]
The ganglion data is same as cyton when hooked to wifi board except that only the
first 4 channels have data, with others specified to be 0
"""
const SAMPLERATE = 250.0         # Hz, the default fs for ganglion
const STARTBYTE = 0xA0           # start of 33 byte data packet
const ENDBYTE = 0xC0             # end of 33 byte data packet
const DEFAULT_STREAM_PORT = 5020 # streaming port, obscure zenginkyo-1 is 5020

# Microvolts per digital unit for Ganglion is 15686 / 8388607 = 0.001869917138805
const INT12MINIMUM = -2048
const INT12MAXIMUM =  2047
const INT24MINIMUM = -8388608    # digital minimum
const INT24MAXIMUM =  8388607    # digital maximum
const CYTONPHYSICALMINIMUM = -187500
const CYTONPHYSICALMAXIMUM =  187500
const GANGLIONPHYSICALMINIMUM = -15686  # microvolts
const GANGLIONPHYSICALMAXIMUM =  15686  # microvolts

# accelerometer conversion factor
const scale_fac_accel_G_per_count = 0.032  # 0.032 G or 32 mG per accelerometer unit

# Commands for board as in SDK http://docs.openbci.com/software/01-Open BCI_SDK
const command_stop = "s"
const command_startText = "x"
const command_startBinary = "b"
const command_startBinary_wAux = "n"
const command_startBinary_4chan = "v"
const command_activateFilters = "F"
const command_deactivateFilters = "g"
const command_deactivate_channel = ["1", "2", "3", "4", "5", "6", "7", "8"]
const command_activate_channel = ["q", "w", "e", "r", "t", "y", "u", "i"]
const command_activate_leadoffP_channel = ["!", "@", "#", "\$", "%", "^", "&", "*"]    # shift + 1-8
const command_deactivate_leadoffP_channel = ["Q", "W", "E", "R", "T", "Y", "U", "I"]  # letters (plus shift) right below 1-8
const command_activate_leadoffN_channel = ["A", "S", "D", "F", "G", "H", "J", "K"]    # letters (plus shift) below the letters below 1-8
const command_deactivate_leadoffN_channel = ["Z", "X", "C", "V", "B", "N", "M", "<"]  # letters (plus shift) below the letters below the letters below 1-8
const command_biasAuto = "`"
const command_biasFixed = "~"
const sratecommands = Dict(25600 =>'0', 12800 =>'1', 6400 =>'2', 3200 =>'3', 1600=>'4', 800=>'5', 400=>'6', 200=>'7')

"""
    postwrite(server, command, saywarn=true)

Send a command to the OpenBCI hardware via the WiFI http server POST /command JSON interface
"""
function postwrite(server, command, saywarn=true)
    resp = post("$server/command",
                "Content-Type"=>"application/json",
                JSON.json(Dict("command"=>command)))
    if resp.status == 200
        info(logger, "command was $command, response was $(String(resp.body))")
    elseif saywarn
        warn(
            logger, "Got status $resp.status, " *
            "failed to get proper response from $server after command $command"
        )
    end
end

""" The functions below all are specific calls to the postwrite function above """
stop(server) = postwrite(server,command_stop)           # stop streaming "s"
softreset(server) = postwrite(server, "v")              # reset peripherals only
start(server) = postwrite(server, command_startBinary)  # 'b', binary streaming should now start
getfs(server) = postwrite(server, "~~")                 # reply is sample rate
setnonstandardfs(server, srate=200) = postwrite(server, sratecommands[srate])
getregisters(server) = postwrite(server, "?")
enablechannels(server, chan=[1,2,3,4]) = for i in chan postwrite(server, b"qwer"[i]) end
disablechannel(server, chan) = postwrite(server,chan)
startsquarewave(server) = postwrite(server, "[")
stopsquarewave(server) = postwrite(server, "]", false)
startimpedancetest(server) = postwrite(server, "z")
stopimpedancetest(server) = postwrite(server, "Z", false)
startaccelerometer(server) = postwrite(server, "n")
stopaccelerometer(server) = postwrite(server, "N", false)
startSDlogging(server) = postwrite(server, "a")         # up to 14 seconds of SD card logging
stopSDlogging(server) = postwrite(server, "j", false)   # stop logging
attachshield(server) = postwrite(server, "{")
detachshield(server) = postwrite(server, "}")
resetshield(server) = postwrite(server,";")

"""
    asyncsocketserver(serveraddress, portnum, packetchannel)

Server task, run in separate task, gets stream of packets from the OpenBCI Wifi hardware
and sends back to main process via a channel
# Arguments
- serveraddress: in form of "http://111.222.133/" (used to restart if error)
- portnum: our port number the OpenBCI Wifi will connect to
- packetchannel: the channel over which we send data to the parent task
"""
function asyncsocketserver(serveraddress, portnum, packetchannel)
    numberofgets = 1
    wifisocket = TCPSocket()
    try
        info(logger, "Entered server async code")
        while true
            server = listen(IPv4(0), portnum)
            wifisocket = accept(server)
            info(logger, "socket service connected to ganglion board")
            bytes = b""
            top = 1
            while isopen(wifisocket)
                bytes = vcat(bytes,read(wifisocket,33))
                if length(bytes) >= 33
                    if bytes[1] == 0xA0 # in sync?
                        put!(packetchannel, bytes[1:33])
                        bytes = bytes[34:end]
                    elseif (top = something(findfirst(x->x==0xA0, bytes), 0)) > 0
                        info(logger, "sync: dropping bytes above position $top")
                        bytes = bytes[top:end]
                    else
                        info(logger, "sync: dumping buffer")
                        bytes = b""
                    end
                else
                    yield()
                    sleep(0.01)
                end
            end
            # if we got here we may need to do a restart, but try a new get also
            numberofgets += 1
            info(logger, "Redoing get for binary stream, will now have done get $numberofgets times")
            get("$serveraddress/stream/start")
        end
    catch y
        info(logger, "Caught exception $y")
        # either error or the channel to process the packets has been closed
        info(logger, "Exiting WiFi streaming task")
        close(wifisocket)
    end
end

"""
    rawOpenBCIboard(ip_board, ip_ours; portnum,fs,latency,locallogging,logSD,useaccelerometer,impedancetest,maketestwave)

Set up the raw OpenBCI WiFI shield connection with the OpenBCI ganglion board.
# Arguments
-ip_board         the ip number of the WiFi shield
-ip_ours          the ip to get the board's stream, usually this computer
--- optional named arguments ---
-portnum          the port number to which the shield will stream data
-fs               sampling rate, usually 250 == SAMPLERATE
-latency          latency in microseconds, time between packets, default 15 msec
-locallogging     true if loglevel is to be info rather than warn
-logSD            true if should log to SD card
-useaccelerometer true if accelerometer data to be sent
-impedancetest    true if impedance check to be done
-maketestwave     true if test square-wave signal to be generated
"""
function rawOpenBCIboard(ip_board, ip_ours; portnum=DEFAULT_STREAM_PORT,
                         fs=250, latency=15000, locallogging=false,
                         logSD=false, useaccelerometer=true,
                         impedancetest=false, maketestwave=false)
    if locallogging
        setlevel!(logger, "info")
        info(logger, "--- Logging starting session ---")
    else
        setlevel!(logger, "warn")
    end
    serveraddress = "http://$ip_board"
    resp = get("$serveraddress/board")
    if resp.status == 200
        # expect: {"board_connected": true, "board_type": "string",
        # "gains": [ null ], "num_channels": 0}
        sysinfo = JSON.parse(String(resp.body))
        info(logger, "board reports $sysinfo")
        num_signals = sysinfo["num_channels"]
        if !(num_signals in (4, 8, 16))
            warn(logger, "board reports $num_signals channels")
        end
    end
    # It's important to start the server process before we set up the TCP port
    # so that the board will find the socket for the connection right away.
    packetchannel = Channel(2400)
    @async(asyncsocketserver(serveraddress, portnum, packetchannel))
    # Now we set up the TCP port connection from the board to our service
    jso = Dict("ip"=>ip_ours, "port"=>portnum, "output"=>"raw", "latency"=>latency)
    resp = post("$serveraddress/tcp", "Content-Type"=>"application/json", JSON.json(jso))
    info(logger, "sending json")
    if resp.status == 200
        tcpinfo = JSON.parse(String(resp.body))
        if haskey(tcpinfo, "connected") && tcpinfo["connected"]
            info(
                logger,
                "Wifi shield TCP server, command connection established, info is $tcpinfo"
            )
        else
            throw("TCP connection failure with $serveraddress")
        end
    else
        warn(logger, "tcp config error status code: $(statuscode(resp))")
    end
    if fs != 250 && fs in [1600, 800, 400, 200]
        setnonstandardfs(serveraddress, fs)
    end
    if logSD
        startSDlogging(serveraddress)
    end
    if useaccelerometer
        startaccelerometer(serveraddress)
    else
        stopaccelerometer(serveraddress)
    end
    info(logger, "completed config of ganglion")
    if maketestwave
        startsquarewave(serveraddress)
    else
        stopsquarewave(serveraddress)
    end
    # now do data collection as a task that will terminate only when we tell it later
    # stop via exception when the channel is closed
    info(logger, "asking for stream")
    get("$serveraddress/stream/start")
    sleep(1)
    if impedancetest
        info(logger, "Impedance check will be done for one second.")
        startimpedancetest(serveraddress)
        sleep(1)
        stopimpedancetest(serveraddress)
    end
    # return the channel where the data packets will be stacked and the number of channels
    packetchannel, num_signals
end

"""
    startBDFPluswritefile

Create the header and file descriptions of the future BDF+ file
"""
function startBDFPluswritefile(
    signalchannels::Int,
    patientID="",
    recording="",
    patientcode="",
    gender="",
    birthdate="",
    patientname="",
    patient_additional="",
    admincode="",
    technician="",
    equipment="",
    recording_additional=""
)
    bdfh = BEDFPlus()

    bdfh.writemode = true
    bdfh.version = "$(EDFPlus.version())"
    bdfh.edf = false
    bdfh.edfplus = false
    bdfh.bdf = false
    bdfh.bdfplus = true
    bdfh.discontinuous = false
    bdfh.filetype = EDFPlus.BDFPLUS
    bdfh.channelcount = signalchannels + 1  # add 1 for the annotation channel
    bdfh.patient = patientID
    bdfh.recording = recording
    bdfh.patientcode = patientcode
    bdfh.gender = gender
    bdfh.birthdate = birthdate
    bdfh.patientname = patientname
    bdfh.patient_additional = patient_additional
    bdfh.admincode = admincode
    bdfh.technician = technician
    bdfh.equipment = equipment
    bdfh.recording_additional = recording_additional

    bdfh.annotationchannel = bdfh.channelcount
    return bdfh
end

"""
    startBDFPluswritefile(json_idfile, signalcount)

json_idfile is a file containing JSON formatted patient information.
"""
function startBDFPluswritefile(json_idfile::String, signalcount::Int)
    # if the json file fails as a data source, we use anonymous defaults
    patientID=""
    recording=""
    patientcode=""
    gender=""
    birthdate=""
    patientname=""
    patient_additional=""
    admincode=""
    technician=""
    equipment=""
    recording_additional=""
    try
        jfh = open(json_idfile, "r")
        dict = JSON.parse(read(jfh, String))
        close(jfh)
        info(logger, "Using patient file ID information: $dict")
        if haskey(dict, "patientID") patientID = dict["patientID"] end
        if haskey(dict, "recording") recording = dict["recording"] end
        if haskey(dict, "patientcode") patientcode = dict["patientcode"] end
        if haskey(dict, "gender") gender = dict["gender"] end
        if haskey(dict, "birthdate") birthdate = uppercase(dict["birthdate"]) end
        if haskey(dict, "patientname") patientname = dict["patientname"] end
        if haskey(dict, "patient_additional") patient_additional = dict["patient_additional"] end
        if haskey(dict, "admincode") admincode = dict["admincode"] end
        if haskey(dict, "technician") technician = dict["technician"] end
        if haskey(dict, "equipment") equipment = dict["equipment"] end
        if haskey(dict, "recording_additional") recording_additional = dict["recording_additional"] end
    catch y
        warn(logger, "Error reading startBDFPluswritefile ID file $json_idfile: $y")
    end
    return startBDFPluswritefile(
        signalcount,
        patientID,
        recording,
        patientcode,
        gender,
        birthdate,
        patientname,
        patient_additional,
        admincode,
        technician,
        equipment,
        recording_additional
    )
end

"""
    setplustimenow(bdfh)

Set time of the BDF+ file data being acquired to current time
"""
function setplustimenow(bdfh)
    datetime = now()
    bdfh.startdate_day = Dates.day(datetime)
    bdfh.startdate_month = Dates.month(datetime)
    bdfh.startdate_year = Dates.year(datetime)
    bdfh.starttime_hour = Dates.hour(datetime)
    bdfh.starttime_minute = Dates.minute(datetime)
    bdfh.starttime_second = Dates.second(datetime)
    bdfh.startdatestring = uppercase(Dates.format(datetime, "dd-u-yyyy"))
    return datetime
end

""" function to help set up BDF+ file header """
function makechannelsignalparam(
    bdfh, records, size, interval, num_signals;
    fs=250,
    labels=[string(locali) for locali in 1:num_signals+1],
    transtype="active electrode",
    physdim="uV",
    prefilter="None"
)
    bdfh.signalparam::Array{ChannelParam,1} = []
    for i in 1:num_signals+1
        parm = ChannelParam()
        parm.label = labels[i]
        parm.transducer = transtype
        parm.physdimension = physdim
        if num_signals < 5
            parm.physmin = GANGLIONPHYSICALMINIMUM
            parm.physmax = GANGLIONPHYSICALMAXIMUM
        else
            parm.physmin = CYTONPHYSICALMINIMUM
            parm.physmax = CYTONPHYSICALMAXIMUM
        end
        parm.digmin = INT24MINIMUM
        parm.digmax = INT24MAXIMUM
        parm.smp_per_record = fs
        parm.bufoffset = (i-1) * 3 * fs + 1
        parm.prefilter= prefilter
        if i == num_signals+1
            parm.annotation = true
            parm.transducer=""
            parm.prefilter = ""
        end
        push!(bdfh.signalparam, parm)
    end
    bdfh.datarecords = records
    bdfh.datarecord_duration = interval
    bdfh.recordsize = size
    bdfh.annotationchannel = num_signals+1
    bdfh.file_duration = Float64(records*interval)
    bdfh.headersize = 256 * (num_signals+2)
    return bdfh
end

"""
    makeBDFplusrecord(rectime, packetchannel, acceldata, reclen, num_channels, daisy=false)

Make a single record from signals at time rectime after start of recording.
Use the last (annotation) channel for the timestamp and accelerometer
annotation data. The record size defaults to be 1 second in duration.
So, if data rate is 250 Hz, for ganglion board's 4 data channels, five
channels total and 3 bytes per datapoint, equals 3 * 250 * 5 = 3750 bytes,
9 cyton channnels is 6750 bytes, and the 17 for the cyton with daisyboard
are then 12750 bytes. The wifi server sends packet data via a driver task,
which fills a Channel (packetchannel) with the data in raw packet form to
read by the main thread. We write one averaged accelerometer reading per
packet. Reclen must be a multiple of 3 * (number of signals + 1).
Returns: one record of length reclen bytes.
# Arguments
-rectime time in seconds since start of recording
-packetchannel the Channel that the server task uses to send data
-acceldata true if acceleromenter data is to be used
-reclen record size, see above for best defaults
-num_channels total number of channels in the record including annotation channel
-daisy true if cyton has a daisyboard (take two packets per 16 channels if true)
"""
function makeBDFplusrecord(
    rectime, packetchannel, acceldata, reclen, num_channels, daisy=false
)
    if reclen % (3*num_channels) != 0
        throw("makeganglionrecord record length $reclen is not a multiple of 3*num_channels")
    end
    siglen = div(reclen, num_channels)
    chan = zeros(UInt8, (num_channels, siglen))
    annotpos = 1
    xaccel = yaccel = zaccel = numaccelpackets = 0
    for sigpos in 1:3:siglen-1
        data = take!(packetchannel)
        if daisy # Cyton board using 16 channel daisy
            if data[2] & 1 == 0  # even packet number, maybe 0
                data = take!(packetchannel) # get an odd numbered packet
            end
            data2 = take!(packetchannel) # second packet for daisy
        end
        # the bigendians and the littleendians are clashing again...
        # OpenBSD boards send bigendian, BDF and EDF files are littleendian
        for i in 1:num_channels-1
            if daisy && i > 8
                j = i - 8  # daisy packet channels 9 through 16
                chan[i, sigpos:sigpos+2] .= reverse(data2[j*3:j*3+2])
            else # all but daisy channels
                chan[i, sigpos:sigpos+2] .= reverse(data[i*3:i*3+2])
            end
        end
        if sigpos == 1
            timestamp = EDFPlus.trimrightzeros(string(rectime))
            chan[num_channels, 1:length(timestamp)+4] =
                unsafe_wrap(Array{UInt8,1}, "+$timestamp\x14\x14\x00")
            annotpos += (length(timestamp) + 4)
        end
        # TODO: we can set the ganglion board to send button press data instead of accel data
        # this would be logged as a button press annotation in that record
        # We assume accelerometer data here is the standard non-time-stamped version
        if acceldata && data[33] == 0xC0 && something(findfirst(data[27:32]), 0) > 0
            xaccel += data[27] >> 8 + data[28]
            yaccel += data[29] >> 8 + data[30]
            zaccel += data[31] >> 8 + data[32]
            numaccelpackets += 1
        end
    end
    if acceldata
        xax = round(32.0 * xaccel / numaccelpackets, 4)
        yax = round(32.0 * yaccel / numaccelpackets, 4)
        zax = round(32.0 * zaccel / numaccelpackets, 4)
        atime = EDFPlus.trimrightzeros(string(rectime + reclen/(SAMPLERATE*15* 2)))
        annot = "+" * atime * "\x14$xax $yax $zax (x,y,z) accelerometer data in 1/1000 G units\x14\x00"
        chan[num_channels, annotpos:annotpos+length(annot)-1] .= unsafe_wrap(Array{UInt8,1}, annot)
    end
    recbytes = b""
    for i in 1:num_channels
        recbytes = vcat(recbytes, chan[i,:])
    end
    rec = Array{Int32,1}(undef, div(reclen,3))
    for i in 1:3:reclen-1
        rec[div(i,3)+1] = Int(reinterpret(EDFPlus.Int24, recbytes[i:i+2])[1])
    end
    return rec
end

"""
    nilfunc(bdfh, pcount, maxrecords)

# Arguments
- bdfh handle to the EEG BDF format data file
- pcount the record number, a counting number (1, 2, ...)
- maxrecords the number of records planned for the file

This default function simply logs the records as received.
This function, a callback, can in practice be used to do realtime fft,
spike detection, BCI, etc. Substitute your function taking those arguments
for this function in order to process the EEG data in real time.
This function is the ``inspector`` function argument in functions below.
"""
nilfunc(bdfh, pcount, maxrecords) = info(logger, "Record $pcount of $maxrecords received.")

"""
    makeganglionbdfplus(path, ip_board, ip_ours, records; idfile,inspector,portnum,recordsize,fs,latency,locallogging,logSD,accelannotations,impedancetest,maketestwave)

Set up and sample streaming ganglion board and write a BDF+ file as output.
# Required arguments
- path             pathname of BDF+ file to be written at end of recording run
- ip_board         the ip number of the WiFi shield
- ip_ours          the ip number of the machine getting the stream, generally this computer
- records          number of records to write, is same as seconds in file with defaults

# Optional named arguments
- idfile           optional JSON file for patient and machine data
- inspector        logging or detection function, called once every BDF+ record
- portnum          the port number to which the shield will stream data
- recordsize       size of each record to write in bytes
- fs               sampling rate, usually 250 == SAMPLERATE
- latency          latency in microseconds, time between packets, default 15 msec
- locallogging     true if loglevel is to be info rather than warn
- logSD            true if should log to SD card for 14 sec
- accelannotations true if accelerometer data to be recorded
- impedancetest    true if impedance check to be done
- maketestwave     true if squarewave test signal to be generated
"""
function makeganglionbdfplus(
    path, ip_board, ip_ours, records=60;
    idfile="",
    inspector=nilfunc,
    portnum=DEFAULT_STREAM_PORT,
    recordsize=3750,
    fs=SAMPLERATE,
    latency=15000,
    locallogging=true,
    logSD=false,
    accelannotations=false,
    impedancetest=false,
    maketestwave=false
)
    bdfh = (idfile == "") ? startBDFPluswritefile(4) : startBDFPluswritefile(idfile, 4)
    packetinterval = recordsize / 3750.0
    makechannelsignalparam(bdfh, records, recordsize, packetinterval, 4)
    bdfh.BDFsignals = zeros(Int32,(records,div(recordsize,3)))
    packetchannel, numsig = rawOpenBCIboard(
        ip_board,
        ip_ours,
        portnum=portnum,
        fs=fs,
        latency=latency,
        locallogging=locallogging,
        logSD=logSD,
        useaccelerometer=accelannotations,
        impedancetest=impedancetest,
        maketestwave=maketestwave
    )
    setplustimenow(bdfh)
    pcount = 0
    packettime = 0.0
    while pcount < records
        rec = makeBDFplusrecord(packettime, packetchannel, accelannotations, recordsize, 5)
        pcount += 1
        inspector(bdfh, pcount, records)
        bdfh.BDFsignals[pcount,:] = rec
        packettime += packetinterval
    end
    EDFPlus.writefile!(bdfh, path)
end

"""
    makecyton8bdfplus(path, ip_board, ip_ours, records; idfile,inspector,portnum,recordsize,fs,latency,locallogging,logSD,accelannotations,impedancetest,maketestwave)

Set up and sample streaming cyton 8-channel board and write a BDF+ file as output.
# Required arguments
- path             pathname of BDF+ file to be written at end of recording run
- ip_board         the ip number of the WiFi shield
- ip_ours          the ip number of the machine getting the stream, generally this computer
- records          number of records to write, is same as seconds in file with defaults

# Optional named arguments
- idfile           optional JSON file for patient and machine data
- inspector        logging or detection function, called once every BDF+ record
- portnum          the port number to which the shield will stream data
- recordsize       size of each record to write in bytes
- fs               sampling rate, usually 250 == SAMPLERATE
- latency          latency in microseconds, time between packets, default 15 msec
- locallogging     true if loglevel is to be info rather than warn
- logSD            true if should log to SD card for 14 sec
- accelannotations true if accelerometer data to be recorded
- impedancetest    true if impedance check to be done
- maketestwave     true if squarewave test signal to be generated
"""
function makecyton8bdfplus(
    path, ip_board, ip_ours, records=60;
    idfile="",
    inspector=nilfunc,
    portnum=DEFAULT_STREAM_PORT,
    recordsize=6750,
    fs=SAMPLERATE,
    latency=15000,
    locallogging=true,
    logSD=false,
    accelannotations=false,
    impedancetest=false,
    maketestwave=false
)
    bdfh = (idfile == "") ? startBDFPluswritefile(8) : startBDFPluswritefile(idfile, 8)
    packetinterval = recordsize / 6750.0
    makechannelsignalparam(bdfh, records, recordsize, packetinterval, 8)
    bdfh.BDFsignals = zeros(Int32,(records,div(recordsize,3)))
    packetchannel, numsig = rawOpenBCIboard(
        ip_board,
        ip_ours,
        portnum=portnum,
        fs=fs,
        latency=latency,
        locallogging=locallogging,
        logSD=logSD,
        useaccelerometer=accelannotations,
        impedancetest=impedancetest,
        maketestwave=maketestwave
    )
    setplustimenow(bdfh)
    pcount = 0
    packettime = 0.0
    while pcount < records
        rec = makeBDFplusrecord(packettime, packetchannel, accelannotations, recordsize, 9)
        pcount += 1
        inspector(bdfh, pcount, records)
        bdfh.BDFsignals[pcount,:] = rec
        packettime += packetinterval
    end
    EDFPlus.writefile!(bdfh, path)
end

"""
    makecyton16bdfplus(path, ip_board, ip_ours, records; idfile,inspector,portnum,recordsize,fs,latency,locallogging,logSD,accelannotations,impedancetest,maketestwave)

Set up cyton 16-channel board with daisy board for 16 signal channels and write a BDF+ file as output.
# Required arguments
- path             pathname of BDF+ file to be written at end of recording run
- ip_board         the ip number of the WiFi shield
- ip_ours          the ip number of the machine getting the stream, generally this computer
- records          number of records to write, is same as seconds in file with defaults

# Optional named arguments
- idfile           optional JSON file for patient and machine data
- inspector        logging or detection function, called once every BDF+ record
- portnum          the port number to which the shield will stream data
- recordsize       size of each record to write in bytes
- fs               sampling rate, usually 250 == SAMPLERATE
- latency          latency in microseconds, time between packets, default 15 msec
- locallogging     true if loglevel is to be info rather than warn
- logSD            true if should log to SD card for 14 sec
- accelannotations true if accelerometer data to be recorded
- impedancetest    true if impedance check to be done
- maketestwave     true if squarewave test signal to be generated
"""
function makecyton16bdfplus(
    path, ip_board, ip_ours, records=60;
    idfile="",
    inspector=nilfunc,
    portnum=DEFAULT_STREAM_PORT,
    recordsize=12750,
    fs=SAMPLERATE,
    latency=15000,
    locallogging=true,
    logSD=false,
    accelannotations=false,
    impedancetest=false,
    maketestwave=false
)
    bdfh = (idfile == "") ? startBDFPluswritefile(16) : startBDFPluswritefile(idfile, 16)
    packetinterval = recordsize / 12750.0
    makechannelsignalparam(bdfh, records, recordsize, packetinterval, 16)
    bdfh.BDFsignals = zeros(Int32,(records,div(recordsize,3)))
    packetchannel, numsig = rawOpenBCIboard(
        ip_board,
        ip_ours,
        portnum=portnum,
        fs=fs,
        latency=latency,
        locallogging=locallogging,
        logSD=logSD,
        useaccelerometer=accelannotations,
        impedancetest=impedancetest,
        maketestwave=maketestwave
    )
    setplustimenow(bdfh)
    pcount = 0
    packettime = 0.0
    while pcount < records
        rec = makeBDFplusrecord(packettime, packetchannel, accelannotations, recordsize, 17)
        pcount += 1
        inspector(bdfh, pcount, records)
        bdfh.BDFsignals[pcount,:] = rec
        packettime += packetinterval
    end
    EDFPlus.writefile!(bdfh, path)
end

end # module


text:
__precompile__()
module SingularIntegralEquations
    using Base, Compat, BandedMatrices, ApproxFun, DualNumbers

export cauchy, cauchyintegral, stieltjes, logkernel,
       stieltjesintegral, hilbert, pseudohilbert, pseudocauchy,
       stieltjesjacobimoment, logjacobimoment

import Base: values,getindex,setindex!,*,.*,+,.+,-,.-,==,<,<=,>,
                >=,./,/,.^,^,\,∪,transpose

import BandedMatrices: bzeros

import Compat: view

import ApproxFun
import ApproxFun: bandinds,SpaceOperator,bilinearform,linebilinearform,dotu,
                  plan_transform,plan_itransform,transform,itransform,transform!,itransform!,
                  rangespace, domainspace, promotespaces, InterlaceOperator, coefficientmatrix,
                  canonicalspace, domain, space, Space, promotedomainspace, promoterangespace, AnyDomain, CalculusOperator,
                  SumSpace,PiecewiseSpace, interlace,Multiplication,ArraySpace,DiagonalArrayOperator,
                  BandedMatrix,ChebyshevDirichlet,PolynomialSpace,AbstractProductSpace,evaluate,order,
                  RealBasis,ComplexBasis,AnyBasis,UnsetSpace,BivariateFun,linesum,complexlength,
                  Fun, ProductFun, LowRankFun, mappoint, JacobiZ,
                  real, UnivariateSpace, RealUnivariateSpace, setdomain, eps, choosedomainspace, isapproxinteger, BlockOperator,
                  ConstantSpace,ReOperator,DirectSumSpace,TupleSpace, ZeroSpace,
                  DiagonalInterlaceOperator, LowRankPertOperator, LaurentDirichlet, setcanonicaldomain,
                  IntervalCurve,PeriodicCurve, reverseorientation, op_eltype, @wrapper, mobius,
                  defaultgetindex, WeightSpace, pochhammer, spacescompatible, ∞, LowRankMatrix, refactorsvd!

# we don't override for Bool and Function to make overriding below easier
# TODO: change when cauchy(f,z,s) calls cauchy(f.coefficients,space(f),z,s)

for OP in (:stieltjes,:stieltjesintegral,:pseudostieltjes)
    @eval begin
        $OP(f::Fun)=$OP(space(f),coefficients(f))
        $OP(f::Fun,z,s...)=$OP(space(f),coefficients(f),z,s...)
        $OP(f::Fun,z,s::Function)=$OP(f,z,s==+)
    end
end

hilbert(f)=Hilbert()*f
hilbert(S,f,z)=hilbert(Fun(f,S))(z)
hilbert(f::Fun,z)=hilbert(space(f),coefficients(f),z)

logkernel(f::Fun,z)=logkernel(space(f),coefficients(f),z)

cauchy(f...)=stieltjes(f...)*(im/(2π))
pseudocauchy(f...)=pseudostieltjes(f...)*(im/(2π))
cauchyintegral(u...)=stieltjesintegral(u...)*(im/(2π))

include("LinearAlgebra/LinearAlgebra.jl")
include("Operators/Operators.jl")
include("FundamentalSolutions/FundamentalSolutions.jl")
include("HypergeometricFunctions/HypergeometricFunctions.jl")

include("JacobiQWeight.jl")
include("JacobiQ.jl")
include("stieltjesmoment.jl")

include("circlecauchy.jl")
include("intervalcauchy.jl")
include("singfuncauchy.jl")

include("vectorcauchy.jl")

include("GreensFun/GreensFun.jl")

include("periodicline.jl")
include("arc.jl")
include("curve.jl")
include("asymptotics.jl")

include("fractals.jl")
include("clustertree.jl")

if isdir(Pkg.dir("TikzGraphs"))
    include("introspect.jl")
end

include("Extras/Extras.jl")

end #module


text:
#= List of Variables

tr      trace dataset
wv      waveform dataset
sta     station dataset
idₛ     station id
idₜ     trace id
stas    list of station names
traces  list of trace names
trₛ     names of traces in sta

idⱼ     S[j] :id
iⱼ      S[j] first index in :x to write
jⱼ      S[j] last index in :x to write
lxⱼ     S[j] value of :t[end,1]
ntⱼ     S[j] number of times in :t[:,1]
nⱼ      S[j] length(:x)
tⱼ      S[j] :t

eᵣ      request window end time [ns from epoch]
nᵣ      request window length in samples
sᵣ      request window start time [ns from epoch]

eₜ      trace end time [ns from epoch]
fs      trace sampling frequency [Hz]
iₜ      trace first index to read
jₜ      trace last index to read
nₜ      trace total number of samples
sₜ      trace start time [ns from epoch]

nₒ      overlap window number of samples
sₒ      overlap window start time [ns from epoch]

j         index to channel in S
v         verbosity
Δ         sampling interval [ns]

because I managed to confuse even myself with the bookkeeping, and that's not
happened since coding SeisIO began
=#

function read_asdf!(  S::GphysData,
                      hdf::String,
                      id::String,
                      s::TimeSpec,
                      t::TimeSpec,
                      msr::Bool,
                      v::Integer  )

  SX = SeisData() # for XML
  idr = isa(id, String) ? id_to_regex(id) : id
  (v > 2) && println("Reading IDs that match ", idr)

  if typeof(s) == String && typeof(t) == String
    d0 = s
    d1 = t
    sᵣ = DateTime(s).instant.periods.value*1000 - dtconst
    eᵣ = DateTime(t).instant.periods.value*1000 - dtconst
  else
    (d0, d1) = parsetimewin(s, t)
    sᵣ = DateTime(d0).instant.periods.value*1000 - dtconst
    eᵣ = DateTime(d1).instant.periods.value*1000 - dtconst
  end
  sᵣ *= 1000
  eᵣ *= 1000
  Δ = 0
  fs = 0.0

  # this nesting is a mess
  netsta = netsta_to_regex(id)
  idr = id_to_regex(id)
  f = h5open(hdf, "r")
  wv = f["Waveforms"]
  stas = names(wv)
  sort!(stas)
  (v > 2) && println("Net.sta found: ", stas)

  for idₛ in stas
    if occursin(netsta, idₛ)
      sta = wv[idₛ]
      traces = names(sta)
      sort!(traces)
      (v > 2) && println("Traces found: ", traces)
      for idₜ in traces
        if idₜ == "StationXML"
          sxml = String(UInt8.(read(sta[idₜ])))
          read_station_xml!(SX, sxml, d0, d1, msr, v)
        elseif occursin(idr, idₜ)
          tr = sta[idₜ]
          nₜ = length(tr)
          sₜ = read(tr["starttime"])
          fs = read(tr["sampling_rate"])

          # convert fs to sampling interval in ns
          Δ = round(Int64, 1.0e9/fs)
          eₜ = sₜ + (nₜ-1)*Δ
          (v > 2) && println("sₜ = ", sₜ,"; eₜ = ", eₜ, "; sᵣ = ", sᵣ, "; eᵣ = ", eᵣ)

          if (sᵣ ≤ eₜ) && (eᵣ ≥ sₜ)
            lxⱼ = 0
            iₜ, jₜ, sₒ = get_trace_bounds(sᵣ, eᵣ, sₜ, eₜ, Δ, nₜ)
            nₒ = jₜ-iₜ+1
            idⱼ = String(split(idₜ, "_", limit=2, keepempty=true)[1])
            j = findid(idⱼ, S.id)
            (v > 2) && println("idⱼ = ", idⱼ, " (found at index in S = ", j, ")")
            nᵣ = div(eᵣ-sᵣ, Δ)+1
            if j == 0
              T = eltype(tr)
              push!(S, SeisChannel(id = idⱼ,
                                   fs = fs,
                                   x = Array{T,1}(undef, nᵣ)))
              j = S.n
              if has(tr, "event_id")
                S.misc[j]["event_id"] = read(tr["event_id"])
              end
            else
              tⱼ = getindex(getfield(S, :t), j)
              ntⱼ = div(lastindex(tⱼ), 2)
              nⱼ = lastindex(S.x[j])
              if ntⱼ > 0
                lxⱼ = getindex(tⱼ, ntⱼ)
                check_for_gap!(S, j, div(sₒ, 1000), nₒ, v)
              end
              if lxⱼ + nₒ > nⱼ
                resize!(S.x[j], lxⱼ + max(nₒ, nᵣ))
              end
              if S.fs[j] == 0.0
                S.fs[j] = fs
              end
            end

            if lxⱼ == 0
              tⱼ = mk_t(nₒ, div(sₒ, 1000))
              setindex!(getfield(S, :t), tⱼ, j)
            end

            iⱼ = lxⱼ+1
            jⱼ = iⱼ+nₒ-1
            load_data!(S.x[j], tr, iₜ:jₜ, iⱼ:jⱼ)
          end
          HDF5.h5d_close(tr)
        end
      end
      HDF5.h5g_close(sta)
    end
  end

  # merge in the XML that we read
  sxml_mergehdr!(S, SX, false, true, v)
  trunc_x!(S)

  # Done
  HDF5.h5g_close(wv)
  close(f)
  return S
end

function read_asdf(hdf::String, id::Union{String,Regex}, s::TimeSpec, t::TimeSpec, msr::Bool, v::Integer)
  S = SeisData()
  read_asdf!(S, hdf, id, s, t, msr, v)
  return S
end


text:
export BrickletLEDStripV2LEDValuesLowLevel
struct BrickletLEDStripV2LEDValuesLowLevel
    value_length::Integer
    value_chunk_offset::Integer
    value_chunk_data::Vector{Integer}
end

export BrickletLEDStripV2SPITFPErrorCount
struct BrickletLEDStripV2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletLEDStripV2Identity
struct BrickletLEDStripV2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletLEDStripV2
"""
Controls up to 2048 RGB(W) LEDs
"""
mutable struct BrickletLEDStripV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletLEDStripV2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_led_strip_v2")
        deviceInternal = package.BrickletLEDStripV2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_led_values_low_level
"""
    $(SIGNATURES)

Sets the RGB(W) values for the LEDs starting from *index*.
You can set at most 2048 RGB values or 1536 RGBW values (6144 byte each).

To make the colors show correctly you need to configure the chip type
(see :func:`Set Chip Type`) and a channel mapping (see :func:`Set Channel Mapping`)
according to the connected LEDs.

If the channel mapping has 3 colors, you need to give the data in the sequence
RGBRGBRGB... if the channel mapping has 4 colors you need to give data in the
sequence RGBWRGBWRGBW...

The data is double buffered and the colors will be transfered to the
LEDs when the next frame duration ends (see :func:`Set Frame Duration`).

Generic approach:

* Set the frame duration to a value that represents the number of frames per
  second you want to achieve.
* Set all of the LED colors for one frame.
* Wait for the :cb:`Frame Started` callback.
* Set all of the LED colors for next frame.
* Wait for the :cb:`Frame Started` callback.
* And so on.

This approach ensures that you can change the LED colors with a fixed frame rate.
"""
function set_led_values_low_level(device::BrickletLEDStripV2, index, value_length, value_chunk_offset, value_chunk_data)
    device.deviceInternal.set_led_values_low_level(index, value_length, value_chunk_offset, value_chunk_data)
end

export get_led_values_low_level
"""
    $(SIGNATURES)

Returns *length* RGB(W) values starting from the
given *index*.

If the channel mapping has 3 colors, you will get the data in the sequence
RGBRGBRGB... if the channel mapping has 4 colors you will get the data in the
sequence RGBWRGBWRGBW...
(assuming you start at an index divisible by 3 (RGB) or 4 (RGBW)).
"""
function get_led_values_low_level(device::BrickletLEDStripV2, index, length)
    return device.deviceInternal.get_led_values_low_level(index, length)
end

export set_frame_duration
"""
    $(SIGNATURES)

Sets the frame duration.

Example: If you want to achieve 20 frames per second, you should
set the frame duration to 50ms (50ms * 20 = 1 second).

For an explanation of the general approach see :func:`Set LED Values`.

Default value: 100ms (10 frames per second).
"""
function set_frame_duration(device::BrickletLEDStripV2, duration)
    device.deviceInternal.set_frame_duration(duration)
end

export get_frame_duration
"""
    $(SIGNATURES)

Returns the frame duration as set by :func:`Set Frame Duration`.
"""
function get_frame_duration(device::BrickletLEDStripV2)
    return device.deviceInternal.get_frame_duration()
end

export get_supply_voltage
"""
    $(SIGNATURES)

Returns the current supply voltage of the LEDs.
"""
function get_supply_voltage(device::BrickletLEDStripV2)
    return device.deviceInternal.get_supply_voltage()
end

export set_clock_frequency
"""
    $(SIGNATURES)

Sets the frequency of the clock.

The Bricklet will choose the nearest achievable frequency, which may
be off by a few Hz. You can get the exact frequency that is used by
calling :func:`Get Clock Frequency`.

If you have problems with flickering LEDs, they may be bits flipping. You
can fix this by either making the connection between the LEDs and the
Bricklet shorter or by reducing the frequency.

With a decreasing frequency your maximum frames per second will decrease
too.
"""
function set_clock_frequency(device::BrickletLEDStripV2, frequency)
    device.deviceInternal.set_clock_frequency(frequency)
end

export get_clock_frequency
"""
    $(SIGNATURES)

Returns the currently used clock frequency as set by :func:`Set Clock Frequency`.
"""
function get_clock_frequency(device::BrickletLEDStripV2)
    return device.deviceInternal.get_clock_frequency()
end

export set_chip_type
"""
    $(SIGNATURES)

Sets the type of the LED driver chip. We currently support the chips

* WS2801,
* WS2811,
* WS2812 / SK6812 / NeoPixel RGB,
* SK6812RGBW / NeoPixel RGBW (Chip Type = WS2812),
* WS2813 / WS2815 (Chip Type = WS2812)
* LPD8806 and
* APA102 / DotStar.
"""
function set_chip_type(device::BrickletLEDStripV2, chip)
    device.deviceInternal.set_chip_type(chip)
end

export get_chip_type
"""
    $(SIGNATURES)

Returns the currently used chip type as set by :func:`Set Chip Type`.
"""
function get_chip_type(device::BrickletLEDStripV2)
    return device.deviceInternal.get_chip_type()
end

export set_channel_mapping
"""
    $(SIGNATURES)

Sets the channel mapping for the connected LEDs.

If the mapping has 4 colors, the function :func:`Set LED Values` expects 4
values per pixel and if the mapping has 3 colors it expects 3 values per pixel.

The function always expects the order RGB(W). The connected LED driver chips
might have their 3 or 4 channels in a different order. For example, the WS2801
chips typically use BGR order, then WS2812 chips typically use GRB order and
the APA102 chips typically use WBGR order.

The APA102 chips are special. They have three 8-bit channels for RGB
and an additional 5-bit channel for the overall brightness of the RGB LED
making them 4-channel chips. Internally the brightness channel is the first
channel, therefore one of the Wxyz channel mappings should be used. Then
the W channel controls the brightness.
"""
function set_channel_mapping(device::BrickletLEDStripV2, mapping)
    device.deviceInternal.set_channel_mapping(mapping)
end

export get_channel_mapping
"""
    $(SIGNATURES)

Returns the currently used channel mapping as set by :func:`Set Channel Mapping`.
"""
function get_channel_mapping(device::BrickletLEDStripV2)
    return device.deviceInternal.get_channel_mapping()
end

export set_frame_started_callback_configuration
"""
    $(SIGNATURES)

Enables/disables the :cb:`Frame Started` callback.
"""
function set_frame_started_callback_configuration(device::BrickletLEDStripV2, enable)
    device.deviceInternal.set_frame_started_callback_configuration(enable)
end

export get_frame_started_callback_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by
:func:`Set Frame Started Callback Configuration`.
"""
function get_frame_started_callback_configuration(device::BrickletLEDStripV2)
    return device.deviceInternal.get_frame_started_callback_configuration()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletLEDStripV2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletLEDStripV2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletLEDStripV2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletLEDStripV2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletLEDStripV2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletLEDStripV2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletLEDStripV2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletLEDStripV2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletLEDStripV2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletLEDStripV2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletLEDStripV2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletLEDStripV2)
    return device.deviceInternal.get_identity()
end

export set_led_values
"""
Sets the RGB(W) values for the LEDs starting from *index*.
You can set at most 2048 RGB values or 1536 RGBW values (6144 byte each).

To make the colors show correctly you need to configure the chip type
(see :func:`Set Chip Type`) and a channel mapping (see :func:`Set Channel Mapping`)
according to the connected LEDs.

If the channel mapping has 3 colors, you need to give the data in the sequence
RGBRGBRGB... if the channel mapping has 4 colors you need to give data in the
sequence RGBWRGBWRGBW...

The data is double buffered and the colors will be transfered to the
LEDs when the next frame duration ends (see :func:`Set Frame Duration`).

Generic approach:

* Set the frame duration to a value that represents the number of frames per
  second you want to achieve.
* Set all of the LED colors for one frame.
* Wait for the :cb:`Frame Started` callback.
* Set all of the LED colors for next frame.
* Wait for the :cb:`Frame Started` callback.
* And so on.

This approach ensures that you can change the LED colors with a fixed frame rate.
"""
function set_led_values(device::BrickletLEDStripV2, index, value)

        index = UInt16(index)
    value = convert(Vector{UInt8}, value)

    if length(value) > 65535
        throw(TinkerforgeInvalidParameterError("Value can be at most 65535 items long"))
    end

    value_length = length(value)
    value_chunk_offset = 0

    if value_length == 0
        value_chunk_data = [0] * 58
        ret = set_led_values_low_level(device, index, value_length, value_chunk_offset, value_chunk_data)
    else
        lock(device.stream_lock) do
            while value_chunk_offset < value_length
                value_chunk_data = create_chunk_data(value, value_chunk_offset, 58, 0)
                ret = set_led_values_low_level(device, index, value_length, value_chunk_offset, value_chunk_data)
                value_chunk_offset += 58
            end
        end
    end

        return ret
end

export get_led_values
"""
Returns *length* RGB(W) values starting from the
given *index*.

If the channel mapping has 3 colors, you will get the data in the sequence
RGBRGBRGB... if the channel mapping has 4 colors you will get the data in the
sequence RGBWRGBWRGBW...
(assuming you start at an index divisible by 3 (RGB) or 4 (RGBW)).
"""
function get_led_values(device::BrickletLEDStripV2, index, length)

        index = UInt16(index)
    length = UInt16(length)

    lock(device.stream_lock) do
        ret = get_led_values_low_level(device, index, length)
            value_length = ret.value_length
        value_out_of_sync = ret.value_chunk_offset != 0

        value_data = ret.value_chunk_data

        while !value_out_of_sync && length(value_data) < value_length
            ret = get_led_values_low_level(device, index, length)
                value_length = ret.value_length
            value_out_of_sync = ret.value_chunk_offset != length(value_data)
            value_data += ret.value_chunk_data
        end

        if value_out_of_sync # discard remaining stream to bring it back in-sync
            while ret.value_chunk_offset + 60 < value_length
                ret = get_led_values_low_level(device, index, length)
                    value_length = ret.value_length
            end

            throw(TinkerforgeStreamOutOfSyncError("Value stream is out-of-sync"))
        end
    end

        return value_data[:value_length]
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletLEDStripV2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using Base: @deprecate, @deprecate_binding

# BEGIN TimeZones 0.9 deprecations

@deprecate build(version::AbstractString, regions; kwargs...) build(version; kwargs...) false

# END TimeZones 0.9 deprecations


text:
using StructJuMP
using Base.Test

@testset "printhook" begin

    numScen = 2
    m = StructuredModel(num_scenarios=numScen)

    @variable(m, 0 <= x <= 1)
    @variable(m, 0 <= y <= 1)

    @constraint(m, x + y == 1)
    @objective(m, :Min, x*x + y)

    for i in 1:numScen
        bl = StructuredModel(parent=m, id=i)
        @variable(bl, w >= 0)
        @constraint(bl, w - x - y <= 1)
        @objective(bl, :Min, w*w + w)
    end

    str = string(m)
    @test contains(str, "Child")
end


text:
module TestExtras
using Test
using CategoricalArrays

const ≅ = isequal

@testset "categorical" begin
    for ord in (false, true)
        x = categorical(["a"], true; ordered=ord)
        @test x isa CategoricalVector{String, UInt8}
        @test x == ["a"]
        @test isordered(x) === ord

        x = categorical(["a"], false; ordered=ord)
        @test x isa CategoricalVector{String, UInt32}
        @test x == ["a"]
        @test isordered(x) === ord
    end
end

@testset "allow_missing argument" begin
    x = categorical(["a", "b", missing])
    levels!(x, ["a"], allow_missing=true)
    @test x ≅ ["a", missing, missing]

    x = cut([1, missing, 100], [1, 2], allow_missing=true)
    @test x ≅ ["[1, 2)", missing, missing]
end

end


text:
"""
    li6(z::ComplexF64)::ComplexF64

Returns the complex 6th order polylogarithm of a complex number `z` of type `ComplexF64`.

Author: Alexander Voigt

License: MIT

# Example
```julia
li6(1.0 + 1.0im)
```
"""
function li6(z::ComplexF64)::ComplexF64
    function clog(z)
        az::Float64 = angle(z)
        return 0.5*log(abs2(z)) + (imag(z) == 0.0 && az < 0.0 ? -az : az)*1.0im
    end

    z6::Float64 = 1.0173430619844491
    bf = (
        1.0                   , -31.0/64.0             ,
        1.5241340877914952e-01, -3.4365555877057613e-02,
        5.7174797239368999e-03, -6.8180453746570645e-04,
        4.9960361948734493e-05, -4.9166051196039048e-07,
       -3.0632975161302164e-07,  1.4414599270849095e-08,
        3.7272438230924107e-09, -3.7300867345487607e-10,
       -5.1246526816085832e-11,  9.0541930956636683e-12,
        6.7381882615512517e-13, -2.1215831150303135e-13,
       -6.8408811719011698e-15,  4.8691178462005581e-15
    )

    if imag(z) == 0.0
        if real(z) == 0.0
            return 0.0 + 0.0im
        end
        if real(z) == 1.0
            return z6 + 0.0im
        end
        if real(z) == -1.0
            return -31.0/32.0*z6 + 0.0im
        end
    end

    nz::Float64  = abs2(z)
    pz::Float64  = angle(z)
    lnz::Float64 = 0.5*log(nz)

    if lnz*lnz + pz*pz < 1.0 # |log(z)| < 1
        v::ComplexF64  = lnz + pz*im
        v2::ComplexF64 = v*v
        c1::Float64 = 1.0369277551433699 # zeta(5)
        c2::Float64 = 0.54116161685556910
        c3::Float64 = 0.20034281719326571
        c4::Float64 = 0.068538919452009435
        c5::ComplexF64 = (137.0/60.0 - clog(-v))/120.0
        c6::Float64 = -1.0/1440.0

        cs = (
            -1.6534391534391534e-05, 2.2964432686654909e-08,
            -9.9413128513657614e-11, 6.6912682653423394e-13,
            -5.7933058574392549e-15
        )

        return z6 + v * c1 +
            v2 * (c2 + v * c3 +
            v2 * (c4 + v * c5 +
            v2 * (c6 +
            v * (cs[1] +
            v2 * (cs[2] +
            v2 * (cs[3] +
            v2 * (cs[4] +
            v2 * (cs[5]))))))))
    end

    (u::ComplexF64, rest::ComplexF64, sgn::Float64) = if nz <= 1.0
        (-clog(1.0 - z), 0.0 + 0.0im, 1.0)
    else # nz > 1.0
        arg::Float64 = pz > 0.0 ? pz - pi : pz + pi
        lmz::ComplexF64 = lnz + arg*im # clog(z)
        lmz2::ComplexF64 = lmz*lmz
        (-clog(1.0 - 1.0/z), -31.0*pi^6/15120.0 + lmz2*(-7.0/720.0*pi^4 + lmz2*(-1.0/144.0*pi^2 - 1.0/720.0*lmz2)), -1.0)
    end

    u2::ComplexF64 = u*u
    u4::ComplexF64 = u2*u2
    u8::ComplexF64 = u4*u4

    rest + sgn * (
       u*bf[1] +
       u2*(bf[2] + u*bf[3]) +
       u4*(bf[4] + u*bf[5] + u2*(bf[6] + u*bf[7])) +
       u8*(bf[8] + u*bf[9] + u2*(bf[10] + u*bf[11]) +
           u4*(bf[12] + u*bf[13] + u2*(bf[14] + u*bf[15]))) +
       u8*u8*(bf[16] + u*bf[17] + u2*bf[18])
    )
end


text:
getline(source::IO, ::UInt8) = Vector{UInt8}(readline(source))
getline(source::IO, ::Base.Chars) = readline(source)

"""
`scan(source, blank; skip, nrow)

Reads fixed wdith format file or stream `source`.
Returns `Vector{UnitRange{Int}}` with autotetected fields in `source`.
Detects only fields that exist in all checked lines.

Parameters:
* `source::Union{IO, AbstractString}`: stream or filename to read from
* `blank::Union{Base.Chars, UInt8}=Base._default_delims`: which characters are considered non-data
* `skip::Int=0`: number of lines to skip at the beginning of the file
* `nrow::Int=0`: number of rows containing data to read (possibly including header);
  `0` means to read all data

  If `blank` is `UInt8` then byte indexing is used, otherwise character indexinx is assumed
"""
function scan(source::IO, blank::Union{Base.Chars, UInt8}=Base._default_delims;
              skip::Int=0, nrow::Int=0)
    for i in 1:skip
        line = readline(source)
    end

    allblank = Int[]
    maxwidth = 0
    firstline = true
    row = 0
    while (row < nrow || nrow == 0) && !eof(source)
        line = getline(source, blank)
        isempty(line) && eof(source) && break
        thisblank = Int[]
        for (i, c) in enumerate(line)
            c in blank && push!(thisblank, i)
        end
        if firstline
            allblank = thisblank
            firstline = false
        else
            allblank = intersect(thisblank, allblank)
        end
        maxwidth = max(maxwidth, length(line))
        row += 1
    end
    # if character at maxwidth character index was not blank
    # add a virtual blank at maxwidth+1
    (isempty(allblank) || allblank[end] < maxwidth) && push!(allblank, maxwidth+1)
    last_blank = 0
    range = UnitRange{Int}[]
    maxwidth == 0 && return range
    for this_blank in allblank
        # do not create zero width columns
        if this_blank > last_blank + 1
            push!(range, (last_blank+1):(this_blank-1))
        end
        last_blank = this_blank
    end
    range
end

function scan(source::AbstractString, blank::Union{Base.Chars, UInt8}=Base._default_delims;
              skip::Int=0, nrow::Int=0)
    open(source) do handle
        scan(handle, blank, skip=skip, nrow=nrow)
    end
end


text:
mutable struct MazeDirected{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Wall, Goal}}
    agent_pos::CartesianIndex{2}
    agent_dir::AbstractDirection
    reward::T
    rng::R
    terminal_reward::T
    goal_pos::CartesianIndex{2}
    done::Bool
end

@generate_getters(MazeDirected)
@generate_setters(MazeDirected)

mutable struct MazeUndirected{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Wall, Goal}}
    agent_pos::CartesianIndex{2}
    reward::T
    rng::R
    terminal_reward::T
    goal_pos::CartesianIndex{2}
    done::Bool
end

@generate_getters(MazeUndirected)
@generate_setters(MazeUndirected)

#####
# Directed
#####

"""
Maze generation uses the iterative implementation of randomized depth-first search from [Wikipedia](https://en.wikipedia.org/wiki/Maze_generation_algorithm#Iterative_implementation).
"""
function MazeDirected(; T = Float32, height = 9, width = 9, rng = Random.GLOBAL_RNG)
    @assert isodd(height) && isodd(width) "height and width must be odd numbers"
    vertical_range = 2:2:height-1
    horizontal_range = 2:2:width-1

    objects = (AGENT, WALL, GOAL)
    world = GridWorldBase(objects, height, width)

    world[WALL, :, :] .= true
    for i in vertical_range, j in horizontal_range
        world[WALL, i, j] = false
    end

    goal_pos = CartesianIndex(height - 1, width - 1)
    world[GOAL, goal_pos] = true

    agent_pos = CartesianIndex(2, 2)
    world[AGENT, agent_pos] = true

    agent_dir = RIGHT
    reward = zero(T)
    terminal_reward = one(T)
    done = false

    env = MazeDirected(world, agent_pos, agent_dir, reward, rng, terminal_reward, goal_pos, done)

    RLBase.reset!(env)

    return env
end

RLBase.state_space(env::MazeDirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const MAZE_DIRECTED_LAYERS = SA.SVector(2, 3)
RLBase.state(env::MazeDirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_agent_dir(env), get_half_size(env), MAZE_DIRECTED_LAYERS)

RLBase.state_space(env::MazeDirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::MazeDirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = (copy(get_grid(env)), get_agent_dir(env))

RLBase.action_space(env::MazeDirected, ::RLBase.DefaultPlayer) = DIRECTED_NAVIGATION_ACTIONS
RLBase.reward(env::MazeDirected, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::MazeDirected) = get_done(env)

function RLBase.reset!(env::MazeDirected{T}) where {T}
    world = get_world(env)
    rng = get_rng(env)
    height = get_height(env)
    width = get_width(env)

    vertical_range = 2:2:height-1
    horizontal_range = 2:2:width-1

    world[AGENT, get_agent_pos(env)] = false
    world[GOAL, get_goal_pos(env)] = false

    world[WALL, :, :] .= true
    for i in vertical_range, j in horizontal_range
        world[WALL, i, j] = false
    end

    generate_maze!(env)

    new_goal_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_goal_pos!(env, new_goal_pos)
    world[GOAL, new_goal_pos] = true

    agent_start_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_agent_pos!(env, agent_start_pos)
    world[AGENT, agent_start_pos] = true

    set_agent_dir!(env, rand(rng, DIRECTIONS))

    set_reward!(env, zero(T))
    set_done!(env, false)

    return nothing
end

function (env::MazeDirected{T})(action::AbstractTurnAction) where {T}
    new_dir = turn(action, get_agent_dir(env))
    set_agent_dir!(env, new_dir)
    world = get_world(env)

    if world[GOAL, get_agent_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

function (env::MazeDirected{T})(action::AbstractMoveAction) where {T}
    world = get_world(env)

    agent_pos = get_agent_pos(env)
    dest = move(action, get_agent_dir(env), agent_pos)

    if !world[WALL, dest]
        world[AGENT, agent_pos] = false
        set_agent_pos!(env, dest)
        world[AGENT, dest] = true
    end

    if world[GOAL, get_agent_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

#####
# Undirected
#####

"""
Maze generation uses the iterative implementation of randomized depth-first search from [Wikipedia](https://en.wikipedia.org/wiki/Maze_generation_algorithm#Iterative_implementation).
"""
function MazeUndirected(; T = Float32, height = 9, width = 9, rng = Random.GLOBAL_RNG)
    @assert isodd(height) && isodd(width) "height and width must be odd numbers"
    vertical_range = 2:2:height-1
    horizontal_range = 2:2:width-1

    objects = (AGENT, WALL, GOAL)
    world = GridWorldBase(objects, height, width)

    world[WALL, :, :] .= true
    for i in vertical_range, j in horizontal_range
        world[WALL, i, j] = false
    end

    goal_pos = CartesianIndex(height - 1, width - 1)
    world[GOAL, goal_pos] = true

    agent_pos = CartesianIndex(2, 2)
    world[AGENT, agent_pos] = true

    reward = zero(T)
    terminal_reward = one(T)
    done = false

    env = MazeUndirected(world, agent_pos, reward, rng, terminal_reward, goal_pos, done)

    RLBase.reset!(env)

    return env
end

RLBase.state_space(env::MazeUndirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const MAZE_UNDIRECTED_LAYERS = SA.SVector(2, 3)
RLBase.state(env::MazeUndirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_half_size(env), MAZE_UNDIRECTED_LAYERS)

RLBase.state_space(env::MazeUndirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::MazeUndirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = copy(get_grid(env))

RLBase.action_space(env::MazeUndirected, player::RLBase.DefaultPlayer) = UNDIRECTED_NAVIGATION_ACTIONS
RLBase.reward(env::MazeUndirected, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::MazeUndirected) = get_done(env)

function RLBase.reset!(env::MazeUndirected{T}) where {T}
    world = get_world(env)
    rng = get_rng(env)
    height = get_height(env)
    width = get_width(env)

    vertical_range = 2:2:height-1
    horizontal_range = 2:2:width-1

    world[AGENT, get_agent_pos(env)] = false
    world[GOAL, get_goal_pos(env)] = false

    world[WALL, :, :] .= true
    for i in vertical_range, j in horizontal_range
        world[WALL, i, j] = false
    end

    generate_maze!(env)

    new_goal_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_goal_pos!(env, new_goal_pos)
    world[GOAL, new_goal_pos] = true

    agent_start_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_agent_pos!(env, agent_start_pos)
    world[AGENT, agent_start_pos] = true

    set_reward!(env, zero(T))
    set_done!(env, false)

    return nothing
end

function (env::MazeUndirected{T})(action::AbstractMoveAction) where {T}
    world = get_world(env)

    agent_pos = get_agent_pos(env)
    dest = move(action, agent_pos)

    if !world[WALL, dest]
        world[AGENT, agent_pos] = false
        set_agent_pos!(env, dest)
        world[AGENT, dest] = true
    end

    if world[GOAL, get_agent_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

#####
# Common
#####

function get_candidate_neighbors(pos::CartesianIndex{2})
    shifts = ((-2, 0), (0, -2), (0, 2), (2, 0))
    return map(shift -> CartesianIndex(pos.I .+ shift), shifts)
end

function generate_maze!(env::Union{MazeDirected, MazeUndirected})
    world = get_world(env)
    rng = get_rng(env)
    height = get_height(env)
    width = get_width(env)

    vertical_range = 2:2:height-1
    horizontal_range = 2:2:width-1

    visited = falses(height, width)
    stack = DS.Stack{CartesianIndex{2}}()

    first_cell = CartesianIndex(rand(rng, vertical_range), rand(rng, horizontal_range))
    visited[first_cell] = true
    push!(stack, first_cell)

    while !isempty(stack)
        current_cell = pop!(stack)
        candidate_neighbors = get_candidate_neighbors(current_cell)
        neighbors = filter(pos -> (pos.I[1] in vertical_range) && (pos.I[2] in horizontal_range), candidate_neighbors)
        unvisited_neighbors = filter(pos -> !visited[pos], neighbors)
        if length(unvisited_neighbors) > 0
            push!(stack, current_cell)

            next_cell = rand(rng, unvisited_neighbors)

            mid_pos = CartesianIndex((current_cell.I .+ next_cell.I) .÷ 2)
            world[WALL, mid_pos] = false

            visited[next_cell] = true
            push!(stack, next_cell)
        end
    end

    return nothing
end