text:
struct Delim{T} <: Dates.AbstractDateToken
    d::T
end

charactercode(::Dates.DatePart{c}) where {c} = c

function Format(f::AbstractString, locale::Dates.DateLocale=Dates.ENGLISH)
    tokens = Dates.AbstractDateToken[]
    prev = ()
    prev_offset = 1

    letters = String(collect(keys(Dates.CONVERSION_SPECIFIERS)))
    for m in eachmatch(Regex("(?<!\\\\)([\\Q$letters\\E])\\1*"), f)
        tran = replace(f[prev_offset:prevind(f, m.offset)], r"\\(.)" => s"\1")

        if !isempty(prev)
            letter, width = prev
            typ = Dates.CONVERSION_SPECIFIERS[letter]

            push!(tokens, Dates.DatePart{letter}(width, isempty(tran)))
        end

        if !isempty(tran)
            push!(tokens, length(tran) == 1 ? Delim(first(tran)) : Delim(tran))
        end

        letter = f[m.offset]
        width = length(m.match)

        prev = (letter, width)
        prev_offset = m.offset + width
    end

    tran = replace(f[prev_offset:lastindex(f)], r"\\(.)" => s"\1")

    if !isempty(prev)
        letter, width = prev
        typ = Dates.CONVERSION_SPECIFIERS[letter]

        push!(tokens, Dates.DatePart{letter}(width, false))
    end

    if !isempty(tran)
        push!(tokens, length(tran) == 1 ? Delim(first(tran)) : Delim(tran))
    end

    return Format(tokens, locale)
end

function Format(f::AbstractString, locale::AbstractString)
    Format(f, Dates.LOCALES[locale])
end

function Format(df::Dates.DateFormat{S, T}) where {S, T}
    N = fieldcount(T)
    dftokens = df.tokens
    tokens = Vector{Dates.AbstractDateToken}(undef, fieldcount(T))
    Base.@nexprs 15 i -> begin
        if i <= N
            @inbounds tok = dftokens[i]
            @inbounds tokens[i] = tok isa Dates.Delim ? Delim(tok.d) : dftokens[i]
        end
    end
    if N > 15
        for i = 16:N
            @inbounds tok = dftokens[i]
            @inbounds tokens[i] = tok isa Dates.Delim ? Delim(tok.d) : dftokens[i]
        end
    end
    return Format(tokens, df.locale)
end

function Base.show(io::IO, df::Format)
    print(io, "Parsers.dateformat\"")
    for t in df.tokens
        _show_content(io, t)
    end
    print(io, '"')
end
Base.Broadcast.broadcastable(x::Format) = Ref(x)

function _show_content(io::IO, d::Dates.DatePart{c}) where c
    for i = 1:d.width
        print(io, c)
    end
end

function _show_content(io::IO, d::Delim{<:AbstractChar})
    if d.d in keys(Dates.CONVERSION_SPECIFIERS)
        for i = 1:1
            print(io, '\\', d.d)
        end
    else
        for i = 1:1
            print(io, d.d)
        end
    end
end

function _show_content(io::IO, d::Delim)
    for c in d.d
        if c in keys(Dates.CONVERSION_SPECIFIERS)
            print(io, '\\')
        end
        print(io, c)
    end
end

function Base.show(io::IO, d::Delim)
    print(io, "Delim(")
    _show_content(io, d)
    print(io, ")")
end

macro dateformat_str(str)
    Format(str)
end

# Standard formats
const ISODateTimeFormat = Format("yyyy-mm-dd\\THH:MM:SS.s")
const ISODateFormat = Format("yyyy-mm-dd")
const ISOTimeFormat = Format("HH:MM:SS.s")
const RFC1123Format = Format("e, dd u yyyy HH:MM:SS")

default_format(::Type{DateTime}) = ISODateTimeFormat
default_format(::Type{Date}) = ISODateFormat
default_format(::Type{Time}) = ISOTimeFormat

maxdigits(d::Dates.DatePart) = d.fixed ? d.width : typemax(Int64)

for c in "yYmdHIMS"
    @eval begin
        @inline function tryparsenext(d::Dates.DatePart{$c}, source, pos, len, b, code)
            return tryparsenext_base10(source, pos, len, b, code, maxdigits(d))
        end
    end
end

@inline function tryparsenext_base10(source, pos, len, b, code, maxdigits)
    x::Int64 = 0
    b -= UInt8('0')
    if b > 0x09
        # character isn't a digit, INVALID value
        code |= INVALID
        b += UInt8('0')
        @goto done
    end
    ndigits = 0
    @inbounds while true
        x = Int64(10) * x + Int64(b)
        ndigits += 1
        pos += 1
        incr!(source)
        if eof(source, pos, len)
            code |= EOF
            @goto done
        end
        b = peekbyte(source, pos) - UInt8('0')
        if b > 0x09 || ndigits == maxdigits
            b += UInt8('0')
            @goto done
        end
    end
    @label done
    return x, pos, b, code
end

ascii_lc(c::UInt8) = c in UInt8('A'):UInt8('Z') ? c + 0x20 : c

function tryparsenext(d::Dates.DatePart{'p'}, source, pos, len, b, code)
    ap = ascii_lc(b)
    pos += 1
    incr!(source)
    if !(ap == UInt8('a') || ap == UInt8('p')) || eof(source, pos, len)
        code |= INVALID_TOKEN
    else
        b = peekbyte(source, pos)
        if ascii_lc(b) != UInt8('m')
            code |= INVALID_TOKEN
        end
        pos += 1
        incr!(source)
        if eof(source, pos, len)
            code |= EOF
        else
            b = peekbyte(source, pos)
        end
    end
    return ap == UInt8('a') ? Dates.AM : Dates.PM, pos, b, code
end

function nextchar(source, pos, len, b)
    u = UInt32(b) << 24
    if !Base.between(b, 0x80, 0xf7)
        pos += 1
        incr!(source)
        return reinterpret(Char, u), pos
    end
    return nextchar_continued(source, pos, len, u)
end

function nextchar_continued(source, pos, len, u)
    u < 0xc0000000 && (pos += 1; incr!(source); @goto ret)
    # first continuation byte
    pos += 1
    incr!(source)
    eof(source, pos, len) && @goto ret
    b = peekbyte(source, pos)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 16
    # second continuation byte
    pos += 1
    incr!(source)
    (eof(source, pos, len)) | (u < 0xe0000000) && @goto ret
    b = peekbyte(source, pos)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 8
    # third continuation byte
    pos += 1
    incr!(source)
    (eof(source, pos, len)) | (u < 0xf0000000) && @goto ret
    b = peekbyte(source, pos)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b); pos += 1; incr!(source)
@label ret
    return reinterpret(Char, u), pos
end

for (tok, fn) in zip("uUeE", Any[Dates.monthabbr_to_value, Dates.monthname_to_value, Dates.dayabbr_to_value, Dates.dayname_to_value])
    @eval function tryparsenext(d::Dates.DatePart{$tok}, source, pos, len, b, code, locale)
        startpos = pos
        while true
            c, pos = nextchar(source, pos, len, b)
            if !isletter(c) || eof(source, pos, len)
                pos -= 1
                break
            end
            b = peekbyte(source, pos)
        end
        val = 0
        if startpos == pos
            code |= INVALID_TOKEN
        else
            if source isa AbstractVector{UInt8}
                word = unsafe_string(pointer(source, startpos), pos - startpos)
            else # source isa IO
                fastseek!(source, startpos - 1)
                word = String(read(source, pos - startpos))
            end
            val = $fn(word, locale)
            if val == 0
                code |= INVALID_TOKEN
            end
        end
        return Int64(val), pos, b, code
    end
end

@inline function tryparsenext(d::Dates.DatePart{'s'}, source, pos, len, b, code)
    ms0, newpos, b, code = tryparsenext_base10(source, pos, len, b, code, maxdigits(d))
    invalid(code) && return ms0, newpos, b, code
    len = newpos - pos
    if len > 3
        ms, r = divrem(ms0, Int64(10) ^ (len - 3))
        if r != 0
            code |= INVALID
        end
    else
        ms = ms0 * Int64(10) ^ (3 - len)
    end
    return ms, newpos, b, code
end

@inline function tryparsenext(d::Delim{<:AbstractChar}, source, pos, len, b, code)
    u = bswap(reinterpret(UInt32, d.d))
    while true
        if b != (u & 0x000000ff)
            code |= INVALID_TOKEN
            break
        end
        u >>= 8
        pos += 1
        incr!(source)
        if eof(source, pos, len)
            code |= EOF | (u == UInt32(0) ? SUCCESS : INVALID_TOKEN)
            break
        end
        b = peekbyte(source, pos)
        u == UInt32(0) && break
    end
    return pos, b, code
end

function tryparsenext(d::Delim{String}, source, pos, len, b, code)
    bytes = codeunits(d.d)
    tlen = length(bytes)
    for dpos = 1:tlen
        @inbounds c = bytes[dpos]
        if b != c
            code |= INVALID_TOKEN
            break
        end
        pos += 1
        incr!(source)
        if eof(source, pos, len)
            code |= EOF | (dpos == tlen ? SUCCESS : INVALID_TOKEN)
            break
        end
        b = peekbyte(source, pos)
    end
    return pos, b, code
end

# fallback that would call custom DatePart overloads that are expecting a string
function tryparsenext(tok, source, pos, len, b, code)::Tuple{Any, Int, UInt8, ReturnCode}
    strlen = min(len - pos + 1, 64)
    str = getstring(source, PosLen(pos, strlen), 0x00)
    res = Dates.tryparsenext(tok, str, 1, strlen)
    if res === nothing
        val = nothing
        code |= INVALID_TOKEN
    else
        val, i = res
        pos += i - 1
        if eof(source, pos, len)
            code |= EOF
        else
            b = peekbyte(source, pos)
        end
    end
    return val, pos, b, code
end

@inline function typeparser(::Type{T}, source, pos, len, b, code, options) where {T <: Dates.TimeType}
    fmt = options.dateformat
    df = fmt === nothing ? default_format(T) : fmt
    tokens = df.tokens
    locale::Dates.DateLocale = df.locale
    year = month = day = Int64(1)
    hour = minute = second = millisecond = Int64(0)
    tz = ""
    @static if VERSION >= v"1.3-DEV"
        ampm = Dates.TWENTYFOURHOUR
    else
        ampm = 0
    end
    extras = nothing
    for tok in tokens
        # @show pos, Char(b), code, typeof(tok)
        eof(code) && break
        if tok isa Delim{Char}
            pos, b, code = tryparsenext(tok, source, pos, len, b, code)
        elseif tok isa Delim{String}
            pos, b, code = tryparsenext(tok, source, pos, len, b, code)
        elseif T !== Time && tok isa Dates.DatePart{'y'}
            year, pos, b, code = tryparsenext(tok, source, pos, len, b, code)
        elseif T !== Time && tok isa Dates.DatePart{'Y'}
            year, pos, b, code = tryparsenext(tok, source, pos, len, b, code)
        elseif T !== Time && tok isa Dates.DatePart{'m'}
            month, pos, b, code = tryparsenext(tok, source, pos, len, b, code)
        elseif T !== Time && tok isa Dates.DatePart{'u'}
            month, pos, b, code = tryparsenext(tok, source, pos, len, b, code, locale)
        elseif T !== Time && tok isa Dates.DatePart{'U'}
            month, pos, b, code = tryparsenext(tok, source, pos, len, b, code, locale)
        elseif T !== Time && tok isa Dates.DatePart{'d'}
            day, pos, b, code = tryparsenext(tok, source, pos, len, b, code)
        elseif T !== Time && tok isa Dates.DatePart{'e'}
            _, pos, b, code = tryparsenext(tok, source, pos, len, b, code, locale)
        elseif T !== Time && tok isa Dates.DatePart{'E'}
            _, pos, b, code = tryparsenext(tok, source, pos, len, b, code, locale)
        elseif T !== Date && tok isa Dates.DatePart{'H'}
            hour, pos, b, code = tryparsenext(tok, source, pos, len, b, code)
        elseif T !== Date && tok isa Dates.DatePart{'I'}
            hour, pos, b, code = tryparsenext(tok, source, pos, len, b, code)
        elseif T !== Date && tok isa Dates.DatePart{'M'}
            minute, pos, b, code = tryparsenext(tok, source, pos, len, b, code)
        elseif T !== Date && tok isa Dates.DatePart{'S'}
            second, pos, b, code = tryparsenext(tok, source, pos, len, b, code)
        elseif T !== Date && tok isa Dates.DatePart{'p'}
            ampm, pos, b, code = tryparsenext(tok, source, pos, len, b, code)
        elseif T !== Date && tok isa Dates.DatePart{'s'}
            millisecond, pos, b, code = tryparsenext(tok, source, pos, len, b, code)
        elseif tok isa Dates.DatePart{'z'}
            tz, pos, b, code = tryparsenext(tok, source, pos, len, b, code)
        elseif tok isa Dates.DatePart{'Z'}
            tz, pos, b, code = tryparsenext(tok, source, pos, len, b, code)
        else
            # non-Dates defined character code
            # allocate extras if not already and parse
            if extras === nothing
                extras = Dict{Type, Any}()
            end
            extraval, pos, b, code = tryparsenext(tok, source, pos, len, b, code)::Tuple{Any, Int, UInt8, ReturnCode}
            extras[Dates.CONVERSION_SPECIFIERS[charactercode(tok)]] = extraval
        end
        if invalid(code)
            if invalidtoken(code)
                code &= ~INVALID_TOKEN
            end
            break
        end
        # @show pos, Char(b), code
    end

    if T === Time
        @static if VERSION >= v"1.3-DEV"
            valid = Dates.validargs(T, hour, minute, second, millisecond, Int64(0), Int64(0), ampm)
        else
            valid = Dates.validargs(T, hour, minute, second, millisecond, Int64(0), Int64(0))
        end
    elseif T === Date
        valid = Dates.validargs(T, year, month, day)
    elseif T === DateTime
        @static if VERSION >= v"1.3-DEV"
            valid = Dates.validargs(T, year, month, day, hour, minute, second, millisecond, ampm)
        else
            valid = Dates.validargs(T, year, month, day, hour, minute, second, millisecond)
        end
    elseif T.name.name === :ZonedDateTime
        valid = Dates.validargs(T, year, month, day, hour, minute, second, millisecond, tz)
    else
        # custom TimeType
        if extras === nothing
            extras = Dict{Type, Any}()
        end
        extras[Year] = year; extras[Month] = month; extras[Day] = day;
        extras[Hour] = hour; extras[Minute] = minute; extras[Second] = second; extras[Millisecond] = millisecond;
        types = Dates.CONVERSION_TRANSLATIONS[T]
        vals = Vector{Any}(undef, length(types))
        for (i, type) in enumerate(types)
            vals[i] = get(extras, type) do
                Dates.CONVERSION_DEFAULTS[type]
            end
        end
        valid = Dates.validargs(T, vals...)
    end
    if invalid(code) || valid !== nothing
        if T.name.name === :ZonedDateTime
            x = T(0, TimeZone("UTC"))
        else
            x = T(0)
        end
        code |= INVALID
    else
        if T === Time
            @static if VERSION >= v"1.3-DEV"
                x = Time(Nanosecond(1000000 * millisecond + 1000000000 * second + 60000000000 * minute + 3600000000000 * (Dates.adjusthour(hour, ampm))))
            else
                x = Time(Nanosecond(1000000 * millisecond + 1000000000 * second + 60000000000 * minute + 3600000000000 * hour))
            end
        elseif T === Date
            x = Date(Dates.UTD(Dates.totaldays(year, month, day)))
        elseif T === DateTime
            @static if VERSION >= v"1.3-DEV"
                x = DateTime(Dates.UTM(millisecond + 1000 * (second + 60 * minute + 3600 * (Dates.adjusthour(hour, ampm)) + 86400 * Dates.totaldays(year, month, day))))
            else
                x = DateTime(Dates.UTM(millisecond + 1000 * (second + 60 * minute + 3600 * hour + 86400 * Dates.totaldays(year, month, day))))
            end
        elseif T.name.name === :ZonedDateTime
            x = T(year, month, day, hour, minute, second, millisecond, tz)
        else
            # custom TimeType
            x = T(vals...)
        end
        code |= OK
    end
    if eof(source, pos, len)
        code |= EOF
    end
    return x, code, pos
end


text:
using LinearAlgebra
using FastLapackInterface
using PolynomialMatrixEquations
#using SolveEyePlusMinusAkronB: EyePlusAtKronBWs, generalized_sylvester_solver!

using LinearAlgebra.BLAS

struct LinearRationalExpectationsWs
    algo::String
    endogenous_nbr::Int64
    exogenous_nbr::Int64
    exogenous_deterministic_nbr::Int64
    forward_indices::Vector{Int64}
    purely_forward_indices::Vector{Int64}
    current_indices::Vector{Int64}
    backward_indices::Vector{Int64}
    both_indices::Vector{Int64}
    static_indices::Vector{Int64}
    dynamic_indices::Vector{Int64}
    non_backward_indices::Vector{Int64}
    current_dynamic_indices::Vector{Int64}
    forward_indices_d::Vector{Int64}
    backward_indices_d::Vector{Int64}
    current_dynamic_indices_d::Vector{Int64}
    current_dynamic_indices_d_j::Vector{Int64}
    exogenous_indices::Vector{Int64}
    static_indices_j::Vector{Int64}
    static_nbr::Int64
    dynamic_nbr::Int64
    forward_nbr::Int64
    backward_nbr::Int64
    both_nbr::Int64
    current_nbr::Int64
    jacobian_static::Matrix{Float64}
    qr_ws::QrWs
    solver_ws::Union{GsSolverWs, CyclicReductionWs}
    a::Matrix{Float64}
    b::Matrix{Float64}
    c::Matrix{Float64}
    d::Matrix{Float64}
    e::Matrix{Float64}
    x::Matrix{Float64}
    ghx::Matrix{Float64}
    gx::Matrix{Float64}
    hx::Matrix{Float64}
    temp1::Matrix{Float64}
    temp2::Matrix{Float64}
    temp3::Matrix{Float64}
    temp4::Matrix{Float64}
    temp5::Matrix{Float64}
    temp6::Matrix{Float64}
    temp7::Matrix{Float64}
    temp8::Matrix{Float64}
    temp9::Matrix{Float64}
    b10::Matrix{Float64}
    b11::Matrix{Float64}
    icolsD::Vector{Int64}
    icolsE::Vector{Int64}
    jcolsD::Vector{Int64}
    jcolsE::Vector{Int64}
    colsUD::Vector{Int64}
    colsUE::Vector{Int64}
    AGplusB::Matrix{Float64}
    linsolve_static_ws::LinSolveWs
    AGplusB_linsolve_ws::LinSolveWs
    #    eye_plus_at_kron_b_ws::EyePlusAtKronBWs

    function LinearRationalExpectationsWs(algo::String,
                                          endogenous_nbr::Int64,
                                          exogenous_nbr::Int64,
                                          exogenous_deterministic_nbr::Int64,
                                          forward_indices::Vector{Int64},
                                          current_indices::Vector{Int64},
                                          backward_indices::Vector{Int64},
                                          both_indices::Vector{Int64},
                                          static_indices::Vector{Int64})
        static_nbr = length(static_indices)
        forward_nbr = length(forward_indices)
        backward_nbr = length(backward_indices)
        both_nbr = length(both_indices)
        current_nbr = length(current_indices)
        dynamic_nbr = endogenous_nbr - static_nbr
        dynamic_indices = setdiff(collect(1:endogenous_nbr), static_indices)
        current_dynamic_indices = setdiff(current_indices, static_indices)
        purely_forward_indices = setdiff(forward_indices, both_indices)
        non_backward_indices = sort(union(purely_forward_indices, static_indices))
        forward_indices_d = findall(in(forward_indices), dynamic_indices)
        backward_indices_d = findall(in(backward_indices), dynamic_indices)
        current_dynamic_indices_d = findall(in(current_dynamic_indices), dynamic_indices)
        current_dynamic_indices_d_j = backward_nbr .+ findall(in(dynamic_indices), current_indices)
        exogenous_indices = backward_nbr + current_nbr + forward_nbr .+ (1:exogenous_nbr)
        jacobian_static = Matrix{Float64}(undef, endogenous_nbr, static_nbr)
        static_indices_j = backward_nbr .+ [findfirst(isequal(x), current_indices) for x in static_indices]
        qr_ws = QrWs(jacobian_static)
        if algo == "GS"
            de_order = forward_nbr + backward_nbr
            d = zeros(de_order, de_order)
            e = similar(d)
            solver_ws = GsSolverWs(d, e, backward_nbr)
            a = Matrix{Float64}(undef, 0, 0)
            b = similar(a)
            c = similar(a)
            x = similar(a)
        elseif algo == "CR"
            a = Matrix{Float64}(undef, dynamic_nbr, dynamic_nbr)
            b = similar(a)
            c = similar(a)
            x = similar(a)
            solver_ws = CyclicReductionWs(dynamic_nbr)
            d = Matrix{Float64}(undef, 0, 0)
            e = similar(d)
        end
        ghx = Matrix{Float64}(undef, endogenous_nbr, backward_nbr)
        gx = Matrix{Float64}(undef, forward_nbr, backward_nbr)
        hx = Matrix{Float64}(undef,  backward_nbr, backward_nbr)
        temp1 = Matrix{Float64}(undef, static_nbr, forward_nbr)
        temp2 = Matrix{Float64}(undef, static_nbr, backward_nbr)
        temp3 = Matrix{Float64}(undef, forward_nbr, backward_nbr)
        temp4 = Matrix{Float64}(undef, endogenous_nbr - static_nbr, backward_nbr)
        temp5 = Matrix{Float64}(undef, endogenous_nbr, exogenous_nbr)
        temp6 = Matrix{Float64}(undef, static_nbr, backward_nbr)
        temp7 = Matrix{Float64}(undef, endogenous_nbr, backward_nbr)
        temp8 = Matrix{Float64}(undef, endogenous_nbr, forward_nbr)
        temp9 = Matrix{Float64}(undef, endogenous_nbr, current_nbr)
        b10 = Matrix{Float64}(undef, static_nbr,static_nbr)
        b11 = Matrix{Float64}(undef, static_nbr, endogenous_nbr - static_nbr)
        current_backward_indices = findall(in(backward_indices), current_indices)
        current_forward_indices = findall(in(forward_indices), current_indices)
        # derivatives of current values of variables that are both
        # forward and backward are included in the D matrix
        k1 = findall(in(current_indices), backward_indices)
        icolsD = [k1;  backward_nbr .+ (1:forward_nbr)]
        jcolsD = [backward_nbr .+ current_backward_indices;
                  backward_nbr + current_nbr .+ (1:forward_nbr)]
        k2a = findall(in(purely_forward_indices), forward_indices)
        k2b = findall(in(purely_forward_indices), current_indices)
        icolsE = [1:backward_nbr; backward_nbr .+ k2a]
        jcolsE = [1:backward_nbr; backward_nbr .+ k2b]
        colsUD = findall(in(forward_indices), backward_indices)
        colsUE = backward_nbr .+ findall(in(backward_indices), forward_indices)
        linsolve_static_ws = LinSolveWs(static_nbr)
        AGplusB = Matrix{Float64}(undef, endogenous_nbr, endogenous_nbr)
        AGplusB_linsolve_ws = LinSolveWs(endogenous_nbr)
        #        if m.serially_correlated_exogenous
        #            eye_plus_at_kron_b_ws = EyePlusAtKronBWs(ma, mb, mc, 1)
        #        else
        #            eye_plus_at_kron_b_ws = EyePlusAtKronBWs(1, 1, 1, 1)
        # end
        new(algo, endogenous_nbr, exogenous_nbr,
            exogenous_deterministic_nbr, forward_indices,
            purely_forward_indices, current_indices, backward_indices,
            both_indices, static_indices, dynamic_indices,
            non_backward_indices, current_dynamic_indices,
            forward_indices_d, backward_indices_d,
            current_dynamic_indices_d, current_dynamic_indices_d_j,
            exogenous_indices, static_indices_j, static_nbr,
            dynamic_nbr, forward_nbr, backward_nbr, both_nbr,
            current_nbr, jacobian_static, qr_ws, solver_ws, a, b, c,
            d, e, x, ghx, gx, hx, temp1, temp2, temp3, temp4, temp5,
            temp6, temp7, temp8, temp9, b10, b11, icolsD, icolsE,
            jcolsD, jcolsE, colsUD, colsUE, AGplusB,
            linsolve_static_ws, AGplusB_linsolve_ws)
    end
end

Base.@kwdef struct CyclicReductionOptions
    maxiter::Int64 = 100
    tol::Float64   = 1e-8
end

Base.@kwdef struct GeneralizedSchurOptions
    # Near unit roots are considered stable roots
    criterium::Float64 = 1.0 + 1e-6
end

Base.@kwdef struct LinearRationalExpectationsOptions
    cyclic_reduction::CyclicReductionOptions = CyclicReductionOptions()
    generalized_schur::GeneralizedSchurOptions = GeneralizedSchurOptions()
end

struct LinearRationalExpectationsResults
    g1::Matrix{Float64}  # full approximation
    gs1::Matrix{Float64} # state transition matrices: states x states
    hs1::Matrix{Float64} # states x shocks
    gns1::Matrix{Float64} # non states x state
    hns1::Matrix{Float64} # non states x shocsk
    # solution first order derivatives w.r. to state variables
    g1_1::SubArray{Float64, 2, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true}
    # solution first order derivatives w.r. to current exogenous variables
    g1_2::SubArray{Float64, 2, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true}
    endogenous_variance::Matrix{Float64}
    #    g1_3::SubArray # solution first order derivatives w.r. to lagged exogenous variables
    stationary_variables::Vector{Bool}
    function LinearRationalExpectationsResults(endogenous_nbr::Int64,
                                               exogenous_nbr::Int64,
                                               backward_nbr::Int64)
        state_nbr = backward_nbr + exogenous_nbr
        non_backward_nbr = endogenous_nbr - backward_nbr
        g1 =  zeros(endogenous_nbr,(state_nbr + 1))
        gs1 = zeros(backward_nbr,backward_nbr)
        hs1 = zeros(backward_nbr,exogenous_nbr)
        gns1 = zeros(non_backward_nbr,backward_nbr)
        hns1 = zeros(non_backward_nbr,exogenous_nbr)
        g1_1 = view(g1, :, 1:backward_nbr)
        g1_2 = view(g1, :, backward_nbr .+ (1:exogenous_nbr))
        endogenous_variance = zeros(endogenous_nbr, endogenous_nbr)
        stationary_variables = Vector{Bool}(undef, endogenous_nbr)
        new(g1, gs1, hs1, gns1, hns1, g1_1, g1_2, endogenous_variance, stationary_variables)
#        g1_3 = view(g[1], :, backward_nbr + exogenous_nbr .+ lagged_exogenous_nbr)
#        new(g, gs, g1_1, g1_2, g1_3, AGplusB, AGplusB_linsolve_ws)
    end
end

"""
remove_static! removes a subset of variables (columns) and rows by QR decomposition
jacobian: on entry jacobian matrix of the original model
          on exit transformed jacobian. The rows corresponding to the dynamic part
                  are at the bottom
p_static: a vector of indices of static variables in jacobian matrix
ws: FirstOrderWs workspace. On exit contains the triangular part conrresponding
                                    to static variables in jacobian_static
"""
function remove_static!(jacobian::Matrix{Float64},
                        ws::LinearRationalExpectationsWs)
    ws.jacobian_static .= view(jacobian, :, ws.static_indices_j)
    geqrf_core!(ws.jacobian_static, ws.qr_ws)
    ormqr_core!('L', ws.jacobian_static', jacobian, ws.qr_ws)
end

"""
Computes the solution for the static variables:
G_y,static = -B_s,s^{-1}(A_s*Gy,fwrd*Gs + B_s,d*Gy,dynamic + C_s)
"""
function add_static!(results::LinearRationalExpectationsResults,
                     jacobian::Matrix{Float64},
                     ws::LinearRationalExpectationsWs)
    # static rows are at the top of the QR transformed Jacobian matrix
    # B_s,s
    fill!(ws.b10, 0.0)
    vj = view(jacobian, 1:ws.static_nbr, ws.static_indices_j)
    copy!(ws.b10, vj)
    # B_s,d
    vb = view(ws.b11, :, ws.current_dynamic_indices_d)
    vj = view(jacobian, 1:ws.static_nbr, ws.current_dynamic_indices_d_j)
    copy!(vb, vj)
    # A_s
    vj = view(jacobian, 1:ws.static_nbr, ws.backward_nbr + ws.current_nbr .+ (1:ws.forward_nbr))
    copy!(ws.temp1, vj)
    # C_s
    vj = view(jacobian, 1:ws.static_nbr, 1:ws.backward_nbr)
    copy!(ws.temp2, vj)
    # Gy.fwrd
    vg = view(results.g1_1, ws.forward_indices, 1:ws.backward_nbr)
    copy!(ws.temp3, vg)
    # Gy.dynamic
    vg = view(results.g1_1, ws.dynamic_indices, 1:ws.backward_nbr)
    copy!(ws.temp4, vg)
    # ws.temp2 = B_s,d*Gy.dynamic + C_s
    mul!(ws.temp2, ws.b11, ws.temp4, 1.0, 1.0)
    # ws.temp6 = A_s*Gy.fwrd*gs1
    mul!(ws.temp6, ws.temp1, ws.temp3)
    mul!(ws.temp2, ws.temp6, results.gs1, -1.0, -1.0)
    # ws.temp3 = B_s,s\ws.temp2
    linsolve_core!(ws.b10, ws.temp2, ws.linsolve_static_ws)
    vg = view(results.g1, ws.static_indices, 1:ws.backward_nbr)
    copy!(vg, ws.temp2)
end

function get_abc!(ws::LinearRationalExpectationsWs,
                  jacobian::AbstractMatrix{Float64})
    fill!(ws.a, 0.0)
    fill!(ws.b, 0.0)
    fill!(ws.c, 0.0)

    va = view(ws.a, :, ws.forward_indices_d)
    vj = view(jacobian, ws.static_nbr .+ (1:ws.dynamic_nbr),
              ws.backward_nbr + ws.current_nbr .+ (1:ws.forward_nbr))
    copy!(va, vj)
    vb = view(ws.b, :, ws.current_dynamic_indices_d)
    vj = view(jacobian, ws.static_nbr .+ (1:ws.dynamic_nbr),
              ws.current_dynamic_indices_d_j)
    copy!(vb, vj)
    vc = view(ws.c, :, ws.backward_indices_d)
    vj = view(jacobian, ws.static_nbr .+ (1:ws.dynamic_nbr),
              1:ws.backward_nbr)
    copy!(vc, vj)
end

function get_de!(ws::LinearRationalExpectationsWs,
                 jacobian::AbstractMatrix{Float64})
    n1 = ws.backward_nbr + ws.forward_nbr - ws.both_nbr
    fill!(ws.d, 0.0)
    fill!(ws.e, 0.0)
    vd = view(ws.d, 1:ws.dynamic_nbr, ws.icolsD)
    vj = view(jacobian, ws.static_nbr .+ (1:ws.dynamic_nbr),
              ws.jcolsD)
    copy!(vd, vj)
    ve = view(ws.e, 1:ws.dynamic_nbr, ws.icolsE)
    vj = view(jacobian, ws.static_nbr .+ (1:ws.dynamic_nbr),
              ws.jcolsE)
    ve .= .-vj
    for i = 1:ws.both_nbr
        k = ws.dynamic_nbr + i
        m = ws.colsUD[i]
        ws.d[k, m] = 1.0
    end
    for i = 1:ws.both_nbr
        k = ws.dynamic_nbr + i
        m = ws.colsUE[i]
        ws.e[k, m] = 1.0
    end
end

"""
Computes LU decomposition of A*G + B
"""
function make_lu_AGplusB!(AGplusB::AbstractMatrix{Float64},
                          A::AbstractMatrix{Float64},
                          G::AbstractMatrix{Float64},
                          B::AbstractMatrix{Float64},
                          ws::LinearRationalExpectationsWs)
    fill!(AGplusB, 0.0)
    vAG = view(AGplusB, :, ws.current_indices)
    copy!(vAG, B)
    vG = view(G, ws.forward_indices, :)
    copy!(ws.temp3, vG)
    mul!(ws.temp7, A, ws.temp3)
    vAG = view(AGplusB, :, ws.backward_indices)
    vAG .+= ws.temp7
    FastLapackInterface.lu!(AGplusB, ws.AGplusB_linsolve_ws)
end

function solve_for_derivatives_with_respect_to_shocks!(results::LinearRationalExpectationsResults,
                                                       jacobian::AbstractMatrix{Float64},
                                                       ws::LinearRationalExpectationsWs)
    #=
    if model.lagged_exogenous_nbr > 0
        f6 = view(jacobian,:,model.i_lagged_exogenous)
        for i = 1:model.current_exogenous_nbr
            for j = 1:model.endo_nbr
                results.g1_3[i,j] = -f6[i,j]
            end
        end
        linsolve_core_no_lu!(results.f1g1plusf2, results.g1_3, ws)
    end
    =#
    if ws.exogenous_nbr > 0
        results.g1_2 .= .-view(jacobian, :, ws.exogenous_indices)

#        if ws.serially_correlated_exogenous
            # TO BE DONE
        #        else
        linsolve_core_no_lu!(ws.AGplusB, results.g1_2, ws.AGplusB_linsolve_ws)
#        end
    end
end

function first_order_solver!(results::LinearRationalExpectationsResults,
                             algo::String,
                             jacobian::AbstractMatrix{Float64},
                             options::LinearRationalExpectationsOptions,
                             ws::LinearRationalExpectationsWs)
    remove_static!(jacobian, ws)
    if algo == "CR"
        get_abc!(ws, jacobian)
        cyclic_reduction!(ws.x, ws.c, ws.b, ws.a, ws.solver_ws,
                          options.cyclic_reduction.tol, options.cyclic_reduction.maxiter)
        vg = view(results.gs1, :, 1:ws.backward_nbr)
        vx = view(ws.x, ws.backward_indices_d, ws.backward_indices_d)
        copy!(vg, vx)
        vg = view(results.g1,ws.dynamic_indices, 1:ws.backward_nbr)
        vx = view(ws.x, :, ws.backward_indices_d)
        copy!(vg, vx)
    elseif algo == "GS"
        get_de!(ws, jacobian)
        gs_solver!(ws.solver_ws, ws.d, ws.e, ws.backward_nbr,
                   options.generalized_schur.criterium)
        results.gs1 .= ws.solver_ws.g1
        vs = view(ws.solver_ws.g1, 1:ws.backward_nbr, 1:ws.backward_nbr)
        vr = view(results.g1, ws.backward_indices,1:ws.backward_nbr)
        copy!(vr,vs)
        vs = view(ws.solver_ws.g2,
                  ws.icolsE[ws.backward_nbr .+ (1:(ws.forward_nbr - ws.both_nbr))]
                  .-ws.backward_nbr, :)
        vr = view(results.g1, ws.purely_forward_indices, 1:ws.backward_nbr)
        copy!(vr,vs)
    else
        error("Algorithm $algo not recognized")
    end
    if ws.static_nbr > 0
        add_static!(results, jacobian, ws)
    end
    #    A = view(jacobian, :, ws.backward_nbr + ws.current_nbr .+ (1:ws.forward_nbr))
    #    B = view(jacobian, :, ws.backward_nbr .+ ws.current_indices)
    vt = view(ws.temp8, :, 1:ws.forward_nbr)
    vj = view(jacobian, :, ws.backward_nbr + ws.current_nbr .+ (1:ws.forward_nbr))
    copy!(vt, vj)
    vt = view(ws.temp9, :, 1:ws.current_nbr)
    vj = view(jacobian, :, ws.backward_nbr .+ (1:ws.current_nbr))
    copy!(vt, vj)
    make_lu_AGplusB!(ws.AGplusB, ws.temp8, results.g1_1, ws.temp9, ws)
    solve_for_derivatives_with_respect_to_shocks!(results, jacobian, ws)
    vh = view(results.g1_2, ws.backward_indices, :)
    results.hs1 .= vh
    vgns = view(results.g1_1, ws.non_backward_indices, :)
    results.gns1 .= vgns
    vhns = view(results.g1_2, ws.non_backward_indices, :)
    results.hns1 .= vhns
end

ext: jl


text:
using Test

module A
    const y = 2
    function foo end
    export y, foo
end

module B
    using ImportAll
    import ..A
    @importall A
    foo(x) = x + 2
end

module C
    using ImportAll
    @importall Base
    struct A end
    iterate(::A) = nothing
    length(::A) = 0
end

@testset "import all" begin
    @test B.y == 2
    @test A.foo(1) == 3
    @test collect(C.A()) == []
end


text:
module pureSubstanceThermo
    import Printf.@sprintf
    include("data.jl")
    include("SolidsEnthalpy.jl")
    include("SolidsCp.jl")
    include("SolidsEntropy.jl")
    include("GasThermo.jl")
    include("MolarMass.jl")
    export GasEntropy, GasEnthalpy, GasCp
    export SolidsCp, SolidsEntropy, SolidsEnthalpy
    export Mm
end  # module SolidsThermo


text:
function LinearEvolution(arg0::RealMatrix, arg1::RealMatrix, arg2::RealVector, arg3::RealMatrix, arg4::RealMatrix)
    return LinearEvolution((RealMatrix, RealMatrix, RealVector, RealMatrix, RealMatrix), arg0, arg1, arg2, arg3, arg4)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_command(obj::LinearEvolution)
    return jcall(obj, "getCommand", RealVector, ())
end

function get_control_matrix(obj::LinearEvolution)
    return jcall(obj, "getControlMatrix", RealMatrix, ())
end

function get_measurement_jacobian(obj::LinearEvolution)
    return jcall(obj, "getMeasurementJacobian", RealMatrix, ())
end

function get_process_noise_matrix(obj::LinearEvolution)
    return jcall(obj, "getProcessNoiseMatrix", RealMatrix, ())
end

function get_state_transition_matrix(obj::LinearEvolution)
    return jcall(obj, "getStateTransitionMatrix", RealMatrix, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
struct Poisson2{T<:AbstractFloat, N, TD<:AbstractArray{Bool, N}, TR8, TR27, TRI}
    dx::NTuple{N, T}
    idx2::NTuple{N, T}
    D::T
    iD::T
    interior::TD
    boundary::TD
    R8::TR8
    R27::TR27
    RI::TRI
end

Poisson2(interior, dx) = Poisson2(Float64, interior, dx)

function Poisson2(
    ::Type{T},
    interior::AbstractArray{Bool, N},
    dx::NTuple{N, Real}
) where {T<:AbstractFloat, N}
    dxT = convert.(T, dx)

    idx2 = inv.(dxT.*dxT)
    D = -2*sum(idx2)

    boundary = boundary_mask(interior)

    ax = map(a -> first(a)+1:last(a)-1, axes(interior))

    tsz = padded_tilesize(T, ntuple(_ -> 2, Val(N)), 1)
    R8 = vec(collect(TileIterator(ax, tsz)))

    tsz = padded_tilesize(T, ntuple(_ -> 3, Val(N)), 1)
    R27 = vec(collect(TileIterator(ax, tsz)))

    i = 0
    RI = Vector{typeof(first(eachindex(interior)))}(undef, sum(interior))
    @inbounds for I in eachindex(interior)
        if interior[I]
            RI[i += 1] = I
        end
    end

    return Poisson2(dxT, idx2, D, inv(D), interior, boundary, R8, R27, RI)
end

Base.size(A::Poisson2) = size(A.interior)
Base.eltype(::Poisson2{T}) where {T} = T

function Base.:(*)(
    A::Poisson2{<:AbstractFloat, 3},
    x::AbstractArray{<:AbstractFloat, 3}
)
    return mul!(tzero(x), A, x)
end

function Base.:(*)(
    A::Poisson2{<:AbstractFloat, 3},
    x::AbstractArray{<:AbstractFloat, 4}
)
    return mul!(tzero(x), A, x)
end

function LinearAlgebra.mul!(
    d2x::AbstractArray{<:AbstractFloat, 3},
    A::Poisson2{<:AbstractFloat, 3},
    x::AbstractArray{T, 3},
) where {T<:AbstractFloat}
    D = convert(T, A.D)
    idx2 = convert(T, A.idx2[1])
    idy2 = convert(T, A.idx2[2])
    idz2 = convert(T, A.idx2[3])

    # @threads ~10% faster here
    # Polyester needs A.interior inside or the loop slows down by 100%???
    # ie. slow down when m = A.interior ... if m[i,j,k] ...
    @inbounds @batch for (I, J, K) in A.R8
        for k in K
            for j in J
                for i in I
                    z1 = x[i,j,k-1]
                    y1 = x[i,j-1,k]
                    x1 = x[i-1,j,k]
                    x0 = x[i,j,k]
                    x2 = x[i+1,j,k]
                    y2 = x[i,j+1,k]
                    z2 = x[i,j,k+1]

                    if A.interior[i,j,k]
                        Δz = z1 + z2
                        Δy = y1 + y2
                        Δx = x1 + x2

                        Δ = D * x0
                        Δ = muladd(idz2, Δz, Δ)
                        Δ = muladd(idy2, Δy, Δ)
                        Δ = muladd(idx2, Δx, Δ)

                        d2x[i,j,k] = Δ
                    end
                end
            end
        end
    end

    return d2x
end

function LinearAlgebra.mul!(
    d2x::AbstractArray{<:AbstractFloat, 4},
    A::Poisson2{<:AbstractFloat, 3},
    x::AbstractArray{T, 4},
) where {T<:AbstractFloat}
    D = convert(T, A.D)
    idx2 = convert(T, A.idx2[1])
    idy2 = convert(T, A.idx2[2])
    idz2 = convert(T, A.idx2[3])

    # @threads ~10% faster here
    # Polyester needs A.interior inside or the loop slows down by 100%???
    # ie. slow down when m = A.interior ... if m[i,j,k] ...
    @inbounds for t in axes(x, 4)
        xt = @view(x[:,:,:,t])
        d2xt = @view(d2x[:,:,:,t])
        @batch for (I, J, K) in A.R8
            for k in K
                for j in J
                    for i in I
                        z1 = xt[i,j,k-1]
                        y1 = xt[i,j-1,k]
                        x1 = xt[i-1,j,k]
                        x0 = xt[i,j,k]
                        x2 = xt[i+1,j,k]
                        y2 = xt[i,j+1,k]
                        z2 = xt[i,j,k+1]

                        if A.interior[i,j,k]
                            Δz = z1 + z2
                            Δy = y1 + y2
                            Δx = x1 + x2

                            Δ = D * x0
                            Δ = muladd(idz2, Δz, Δ)
                            Δ = muladd(idy2, Δy, Δ)
                            Δ = muladd(idx2, Δx, Δ)

                            d2xt[i,j,k] = Δ
                        end
                    end
                end
            end
        end
    end

    return d2x
end

function residual!(
    r::AbstractArray{<:AbstractFloat, 3},
    b::AbstractArray{<:AbstractFloat, 3},
    A::Poisson2{<:AbstractFloat, 3},
    x::AbstractArray{T, 3},
) where {T<:AbstractFloat}
    D = convert(T, -A.D)
    idx2 = convert(T, -A.idx2[1])
    idy2 = convert(T, -A.idx2[2])
    idz2 = convert(T, -A.idx2[3])

    # @threads ~10% faster here
    # Polyester needs A.interior inside or the loop slows down by 100%???
    # ie. slow down when m = A.interior ... if m[i,j,k] ...
    @inbounds @batch for (I, J, K) in A.R8
        for k in K
            for j in J
                for i in I
                    z1 = x[i,j,k-1]
                    y1 = x[i,j-1,k]
                    x1 = x[i-1,j,k]
                    x0 = x[i,j,k]
                    x2 = x[i+1,j,k]
                    y2 = x[i,j+1,k]
                    z2 = x[i,j,k+1]

                    if A.interior[i,j,k]
                        Δz = z1 + z2
                        Δy = y1 + y2
                        Δx = x1 + x2

                        Δ = muladd(D, x0, b[i,j,k])
                        Δ = muladd(idz2, Δz, Δ)
                        Δ = muladd(idy2, Δy, Δ)
                        Δ = muladd(idx2, Δx, Δ)

                        r[i,j,k] = Δ
                    end
                end
            end
        end
    end

    return r
end

function __norm_residual(
    b::AbstractArray{<:AbstractFloat, 3},
    A::Poisson2{<:AbstractFloat, 3},
    x::AbstractArray{T, 3},
) where {T<:AbstractFloat}
    m = A.interior

    D = convert(T, -A.D)
    idx2 = convert(T, -A.idx2[1])
    idy2 = convert(T, -A.idx2[2])
    idz2 = convert(T, -A.idx2[3])

    @inbounds @batch threadlocal=zero(T)::T for (I, J, K) in A.R8
        for k in K
            for j in J
                for i in I
                    z1 = x[i,j,k-1]
                    y1 = x[i,j-1,k]
                    x1 = x[i-1,j,k]
                    x0 = x[i,j,k]
                    x2 = x[i+1,j,k]
                    y2 = x[i,j+1,k]
                    z2 = x[i,j,k+1]

                    if m[i,j,k]
                        Δz = z1 + z2
                        Δy = y1 + y2
                        Δx = x1 + x2

                        Δ = muladd(D, x0, b[i,j,k])
                        Δ = muladd(idz2, Δz, Δ)
                        Δ = muladd(idy2, Δy, Δ)
                        Δ = muladd(idx2, Δx, Δ)

                        threadlocal = muladd(Δ, Δ, threadlocal)
                    end
                end
            end
        end
    end

    return sqrt(sum(threadlocal::Vector{T}))
end

boundary_mask(m::AbstractArray{Bool, 3}) = boundary_mask!(tzero(m), m)

function boundary_mask!(
    mb::AbstractArray{Bool, 3},
    m::AbstractArray{Bool, 3}
)
    sz = size(m)
    nx, ny, nz = sz
    nxc, nyc, nzc = restrict_size(sz)

    b = tzero(mb)

    # array boundary
    outer = CartesianIndices(ntuple(n -> 2:sz[n]-1, Val(3)))
    inner = CartesianIndices(ntuple(n -> 3:sz[n]-2-Int(iseven(sz[n])), Val(3)))
    if !isempty(inner)
        @inbounds for I in EdgeIterator(outer, inner)
            b[I] = m[I]
        end
    end

    # interior
    @inbounds @batch minbatch=8 for kc in 2:nzc-1
        k = (kc << 1) - 1
        for jc in 2:nyc-1
            j = (jc << 1) - 1
            for ic in 2:nxc-1
                i = (ic << 1) - 1

                if !m[i,j,k]     ||
                   !m[i+1,j,k]   ||
                   !m[i,j+1,k]   ||
                   !m[i+1,j+1,k] ||
                   !m[i,j,k+1]   ||
                   !m[i+1,j,k+1] ||
                   !m[i,j+1,k+1] ||
                   !m[i+1,j+1,k+1]

                    b[i,j,k]       = m[i,j,k]
                    b[i+1,j,k]     = m[i+1,j,k]
                    b[i,j+1,k]     = m[i,j+1,k]
                    b[i+1,j+1,k]   = m[i+1,j+1,k]
                    b[i,j,k+1]     = m[i,j,k+1]
                    b[i+1,j,k+1]   = m[i+1,j,k+1]
                    b[i,j+1,k+1]   = m[i,j+1,k+1]
                    b[i+1,j+1,k+1] = m[i+1,j+1,k+1]
                end
            end
        end
    end

    @inbounds @batch minbatch=8 for k in 2:nz-1
        for j in 2:ny-1
            for i in 2:nx-1
                if m[i,j,k]
                    mb[i,j,k] =
                        b[i,j,k-1] ||
                        b[i,j-1,k] ||
                        b[i-1,j,k] ||
                        b[i,j,k] ||
                        b[i+1,j,k] ||
                        b[i,j+1,k] ||
                        b[i,j,k+1]
                end
            end
        end
    end

    return mb
end

function tfill!(x::AbstractArray{T, N}, v, A::Poisson2{<:AbstractFloat, N}) where {T, N}
    vT = convert(T, v)
    @inbounds @batch minbatch=1024 for I in A.RI
        x[I] = vT
    end
    return x
end


text:
include("adaptation/Adaptation.jl")

using .Adaptation

function PreConditioner(::UnitEuclideanMetric)
    return UnitPreConditioner()
end

function PreConditioner(m::DiagEuclideanMetric)
    return DiagPreConditioner(m.dim)
end

function PreConditioner(m::DenseEuclideanMetric)
    return DensePreConditioner(m.dim)
end

function update(h::Hamiltonian, prop::AbstractProposal, dpc::Adaptation.AbstractPreConditioner)
    return h(getM⁻¹(dpc)), prop
end

function update(h::Hamiltonian, prop::AbstractProposal, da::NesterovDualAveraging)
    return h, prop(prop.integrator(getϵ(da)))
end

function update(h::Hamiltonian, prop::AbstractProposal, ca::Adaptation.AbstractCompositeAdaptor)
    return h(getM⁻¹(ca.pc)), prop(prop.integrator(getϵ(ca.ssa)))
end


text:
# an abstract named tuple type
abstract Attributed

# as a tuple
Base.length(a::Attributed) = nfields(a)
Base.start(a::Attributed) = 1
Base.next(a::Attributed, s) = getfield(a, s), s+1
Base.done(a::Attributed, s) = s > nfields(a)

function Base.isequal{A<:Attributed}(a1::A, a2::A)
    for (el1, el2) in zip(a1,a2)
        isequal(el1, el2) || return false
    end
    return true
end

const attributedhash_seed = UInt === UInt64 ? 0x8d02b517474715f1 : 0xf5e5be51
@generated function Base.hash{A<:Attributed}(a::A, h::UInt)
    tuple_ex = Expr(:tuple)
    for j in 1:length(A.parameters)
        push!(tuple_ex.args, Expr(:call, :getfield, :a, j))
    end
    return quote
        $(Expr(:meta, :inline))
        hash($tuple_ex, hash(h, attributedhash_seed))
    end
end

Base.show(io::IO, a::Attributed) = Base.show_delim_array(io, a, '(', ',', ')', true)


text:
using DataFrames
using Gadfly

# Input parameters
path = "/scratch/cns7ae/ciris/sigtest/"
protonFile = "proton_sigmas.csv"
electronFile = "electron_sigmas.csv"
plotname = "cross_sections.pdf"
xmin = 1
xmax = 1e6

plottheme = Gadfly.Theme(
                        line_width=2pt,
                        minor_label_font_size=13pt,
                        minor_label_color=colorant"black",
                        major_label_font_size=15pt,
                        major_label_color=colorant"black",
                        key_label_color=colorant"black",
                        key_title_color=colorant"black",
                        key_title_font_size=15pt,
                        key_label_font_size=13pt,
                        grid_color=colorant"black")

println("Reading input file")
σ_proton = readtable(
                  path*protonFile,
                  names = [:energy, :σ, :Type, :Species],
                  eltypes = [Float64, Float64, String, String],
                  header = false
                 )

σ_electron = readtable(
                  path*electronFile,
                  names = [:energy, :σ, :Type, :Species],
                  eltypes = [Float64, Float64, String, String],
                  header = false
                 )

#println("Now filtering out Total elastic cross-section")

println("Now plotting data")

σ_proton = σ_proton[σ_proton[:σ].>10^(-20.0),:]
σ_proton = σ_proton[σ_proton[:σ].<10^(-15.0),:]
σ_electron = σ_electron[σ_electron[:σ].>10^(-20.0),:]
σ_electron = σ_electron[σ_electron[:σ].<10^(-15.0),:]

red_sig_elec = DataFrame()
red_sig_elec = vcat(red_sig_elec,
                     σ_electron[σ_electron[:Type].=="Elastic",:],
                     σ_electron[σ_electron[:Type].=="Excitation",:],
                     σ_electron[σ_electron[:Type].=="Ionization",:])

σ_electron = red_sig_elec

sort!(σ_proton,cols=[:Type])
sort!(σ_electron,cols=[:Type])

sigo = plot(
           σ_proton[σ_proton[:Species].=="O",:],
           x="energy",
           y="σ",
           color="Type",
           Geom.line,
           Scale.x_log10,
           Scale.y_log10(minvalue=10^(-20.0),maxvalue=10^(-15.0)),
           Guide.xlabel("Particle Energy (eV)"),
           Guide.ylabel("σ (cm2)"),
#           Guide.title("Proton Cross Sections vs. Particle Energy"),
           Theme(default_point_size=1.5pt)
          )

push!(sigo,plottheme)

sigo2 = plot(
           σ_proton[σ_proton[:Species].=="O2",:],
           x="energy",
           y="σ",
           color="Type",
           Geom.line,
           Scale.x_log10,
           Scale.y_log10(minvalue=10^(-20.0),maxvalue=10^(-15.0)),
           Guide.xlabel("Particle Energy (eV)"),
           Guide.ylabel("σ (cm2)"),
#           Guide.title("Proton Cross Sections vs. Particle Energy"),
           Theme(default_point_size=1.5pt)
          )

push!(sigo2,plottheme)

sigo3 = plot(
           σ_proton[σ_proton[:Species].=="O3",:],
           x="energy",
           y="σ",
           color="Type",
           Geom.line,
           Scale.x_log10,
           Scale.y_log10(minvalue=10^(-20.0),maxvalue=10^(-15.0)),
           Guide.xlabel("Particle Energy (eV)"),
           Guide.ylabel("σ (cm2)"),
#           Guide.title("Proton Cross Sections vs. Particle Energy"),
           Theme(default_point_size=1.5pt)
          )

push!(sigo3,plottheme)

esigo = plot(
           σ_electron[σ_electron[:Species].=="O",:],
           x="energy",
           y="σ",
           color="Type",
           Geom.line,
           Scale.x_log10,
           Scale.y_log10(minvalue=10^(-20.0),maxvalue=10^(-15.0)),
           Guide.xlabel("Particle Energy (eV)"),
           Guide.ylabel("σ (cm2)"),
#           Guide.title("Electron Cross Sections vs. Particle Energy"),
           Theme(default_point_size=1.5pt)
          )

push!(esigo,plottheme)

esigo2 = plot(
           σ_electron[σ_electron[:Species].=="O2",:],
           x="energy",
           y="σ",
           color="Type",
           Geom.line,
           Scale.x_log10,
           Scale.y_log10(minvalue=10^(-20.0),maxvalue=10^(-15.0)),
           Guide.xlabel("Particle Energy (eV)"),
           Guide.ylabel("σ (cm2)"),
#           Guide.title("Electron Cross Sections vs. Particle Energy"),
           Theme(default_point_size=1.5pt)
          )

push!(esigo2,plottheme)

esigo3 = plot(
           σ_electron[σ_electron[:Species].=="O3",:],
           x="energy",
           y="σ",
           color="Type",
           Geom.line,
           Scale.x_log10,
           Scale.y_log10(minvalue=10^(-20.0),maxvalue=10^(-15.0)),
           Guide.xlabel("Particle Energy (eV)"),
           Guide.ylabel("σ (cm2)"),
#           Guide.title("Electron Cross Sections vs. Particle Energy"),
           Theme(default_point_size=1.5pt)
          )

push!(esigo3,plottheme)

println("Now exporting plot")
draw(PDF(path*"f1a.pdf", 6inch, 4.5inch), sigo)
draw(PDF(path*"f1b.pdf", 6inch, 4.5inch), sigo2)
draw(PDF(path*"f1c.pdf", 6inch, 4.5inch), sigo3)
draw(PDF(path*"f2a.pdf", 6inch, 4.5inch), esigo)
draw(PDF(path*"f2b.pdf", 6inch, 4.5inch), esigo2)
draw(PDF(path*"f2c.pdf", 6inch, 4.5inch), esigo3)
#draw(PDF(path*"f2a.pdf", 6inch, 4.5inch), p1)
#draw(PDF(path*"f2b.pdf", 6inch, 4.5inch), p2)
#draw(PDF(path*plotname, 6inch, 9inch), vstack(p1,p2))
println("Now ending script")


text:
abstract type PreallocatedDataNormal{T <: Real, S} <: PreallocatedData{T, S} end

abstract type PreallocatedDataNormalKrylov{T <: Real, S} <: PreallocatedDataNormal{T, S} end

include("K1Krylov.jl")


text:
#=
The number 145 is well known for the property that the sum of the factorial of its digits is equal to 145:

1! + 4! + 5! = 1 + 24 + 120 = 145

Perhaps less well known is 169, in that it produces the longest chain of numbers that link back to 169; it turns out that there are only three such loops that exist:

169 → 363601 → 1454 → 169
871 → 45361 → 871
872 → 45362 → 872

It is not difficult to prove that EVERY starting number will eventually get stuck in a loop. For example,

69 → 363600 → 1454 → 169 → 363601 (→ 1454)
78 → 45360 → 871 → 45361 (→ 871)
540 → 145 (→ 145)

Starting with 69 produces a chain of five non-repeating terms, but the longest non-repeating chain with a starting number below one million is sixty terms.

How many chains, with a starting number below one million, contain exactly sixty non-repeating terms?
=#
function sumOfFactorialOfDigits(x)
  total = sum(map(factorial, digits(x)))
end

#TODO: Kinda slowish, could be sped up by memoizing each start number -> count of terms
function lengthOfLoop(x)
  arr = Int64[]
  last = x
  while true
    current = sumOfFactorialOfDigits(last)
    if in(current, arr)
      return length(arr) + 1
    end
    push!(arr, current)
    last = current
  end
end

function calc()
  count = 0
  for i in 69:999_999
    if lengthOfLoop(i) == 60
      count = count + 1
    end
  end
  count
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
# # Statistical plots
#
# Statistical plots are supported automatically, provided that they are implemented
# in Makie. Integration with the styling framework should also happen automatically.
# Some particular statistical plots have specific keyword options that can be applied
# in mapping, such as `dodge` (for `boxplot` and `violin`) or `side` (for `violin`).
#
# ## Examples

using RDatasets: dataset
using AlgebraOfGraphics, CairoMakie
mpg = dataset("ggplot2", "mpg");
mpg.IsAudi = mpg.Manufacturer .== "audi"
geom = visual(BoxPlot, layout_x = 1) + visual(Violin, layout_x = 2)
data(mpg) *
    mapping(:Cyl => categorical, :Hwy) *
    mapping(dodge = :IsAudi => categorical, color = :IsAudi => categorical) *
    geom |> draw
AbstractPlotting.save("boxplot.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](boxplot.svg)

data(mpg) *
    mapping(:Cyl => categorical, :Hwy) *
    mapping(side = :IsAudi => categorical, color = :IsAudi => categorical) *
    visual(Violin) |> draw
AbstractPlotting.save("violin.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](violin.svg)

data(mpg) *
    mapping(:Displ, :Hwy, layout_x = :IsAudi => categorical) *
    visual(QQPlot, qqline = :fit) |> draw
AbstractPlotting.save("qqplot.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](qqplot.svg)


text:
struct ConstantTexture <: Texture
    color::AbstractVec
end

function value(texture::ConstantTexture, u::AbstractFloat, v::AbstractFloat, point::AbstractVec)::AbstractVec
    texture.color
end

max_stars_repo_name: goooxu/JuliaRayTracing

max_issues_repo_name: goooxu/JuliaRayTracing

max_forks_repo_name: goooxu/JuliaRayTracing


text:
###################################################################
#  configuration options
#
#  config(), config_cell(),
###################################################################

########### top-level configuration ##############

config_spec = Dict(
  :viewport     => Vector,
  :background   => AbstractString,
  :timeFormat   => AbstractString,
  :numberFormat => AbstractString
)

"""
Sets top-level properties
config(;properties...)

- `viewport` (`Vector{Real}`) : The width and height of the on-screen viewport, in pixels. If necessary, clipping and scrolling will be applied. Default value: (none)
- `background` (`AbstractString`) : CSS color property to use as background of visualization. Default value: (none)
- `timeFormat` (`AbstractString`) : The default time format pattern for text and labels of axes and legends (in the form of D3 time format pattern). Default value: '%Y-%m-%d'.
- `numberFormat` (`AbstractString`) : The default number format pattern for text and labels of axes and legends (in the form of D3 number format pattern). Default value: 's'.
"""
config(;properties...) = _mkvis((:config, ), config_spec, properties)

########### cell-level configuration ##############

config_cell_spec = Dict(
  :width            =>	Real,
  :height           =>	Real,
  :fill             =>	AbstractString,
  :fillOpacity      =>	Real,
  :stroke           =>	AbstractString,
  :strokeOpacity    =>	Real,
  :strokeWidth      =>	Real,
  :strokeDash       =>	Vector,
  :strokeDashOffset =>	Vector
)

cell_help = """

- `width` (`Real`) :	The width of the single plot or each plot in a trellis plot when the visualization has continuous x-scale. (If the plot has ordinal x-scale, the width is determined by the x-scale’s bandSize and the cardinality of the x-scale. If the plot does not have a field on x, the width is derived from scale config’s bandSize for all marks except text and from scale config’s textBandWidth for text mark.) Default value: 200
- `height` (`Real`) :	The height of the single plot or each plot in a trellis plot when the visualization has continuous y-scale. (If the visualization has ordinal y-scale, the height is determined by the bandSize and the cardinality of the y-scale. If the plot does not have a field on y, the height is scale config’s bandSize.) Default value: 200
- `fill` (`AbstractString`) :	The fill color. Default value: (none)
- `fillOpacity` (`Real`) :	The fill opacity (value between [0,1]). Default value: (none)
- `stroke` (`AbstractString`) :	The stroke color. Default value: (none)
- `strokeOpacity` (`Real`) :	The stroke opacity (value between [0,1]). Default value: (none)
- `strokeWidth` (`Number`) :	The stroke width, in pixels. Default value: (none)
- `strokeDash` (`Vector{Real}`) : An array of alternating stroke, space lengths for creating dashed or dotted lines. Default value: (none)
- `strokeDashOffset` (`Vector{Real}`) : The offset (in pixels) into which to begin drawing with the stroke dash array. Default value: (none)
"""

"""
Sets cell properties
`config_cell(;properties...)`

$cell_help
"""
config_cell(;properties...) = _mkvis((:config, :cell), config_cell_spec, properties)

########### facet configuration (general, axis, scale) ##############

"""
Sets facet cell properties,
`config_facet_cell(;properties...)`

$cell_help
"""
config_facet_cell(;properties...) = _mkvis((:config, :facet, :cell), config_cell_spec, properties)

"""
Sets facet scale properties,
`config_facet_scale(;properties...)`

$scale_help
"""
config_facet_scale(;properties...) = _mkvis((:config, :facet, :scale), config_scale_spec, properties)

"""
Sets facet axis properties, see ? config_axis
`config_facet_axis(;properties...)`

$axis_help
"""
config_facet_axis(;properties...) = _mkvis((:config, :facet, :axis), config_axis_spec, properties)

config_grid_spec = Dict(
  :gridColor   => AbstractString,
  :gridOpacity => Real,
  :gridOffset  => Real
)

"""
Sets facet grid properties

config_grid(;properties...)

- `gridColor` (`AbstractString`) : Color of the grid between facets.
- `gridOpacity` (`Real`) : Opacity of the grid between facets.
- `gridOffset` (`Real`) : Offset for grid between facets.
"""
config_grid(;properties...) = _mkvis((:config, :facet, :grid), config_grid_spec, properties)


text:
function ComputeSourceFunctionOverGrid(
    LensingSourceFunction::LensingSourceFunction,
    ConvolvedDensity::AbstractConvolvedDensity,
    CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities,
    Cosmology::AbstractCosmology,
    PowerSpectrum::PowerSpectrum)
    c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
                       #physical constants involved in calculations
    for zbinidx in 1:length(ConvolvedDensity.ZBinArray)-1
        for zidx in 1:length(CosmologicalGrid.ZArray)
            LensingSourceFunction.SourceFunctionArray[zbinidx, zidx] =
            1.5 * (Cosmology.H0/c_0)^2 * Cosmology.ΩM *
            (1. + CosmologicalGrid.ZArray[zidx]) *
            LensingSourceFunction.LensingEfficiencyArray[zbinidx, zidx] *
            BackgroundQuantities.rZArray[zidx]
        end
    end
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

function nmac_clusters(clustering::DataFrame, Ds::DFSet)
    meta = getmeta(Ds)
    D = join(meta, clustering, on=[:encounter_id])
    D = D[D[:nmac], :] #nmacs only
    ids = D[:id]
    labels = D[:label]
    Dl = DFSetLabeled(Ds[ids], labels)
    Dl
end

function nonnmacs_extra_cluster(clustering::DataFrame, Ds::DFSet)
    Dl1 = nmac_clusters(clustering, Ds) #cluster labeled nmacs

    meta = getmeta(Ds)
    D = meta[!meta[:nmac], :] #non-nmacs only
    extra_label = maximum(clustering[:label]) + 1 #highest cluster label + 1
    ids = D[:id]
    labels = fill(extra_label, length(ids))
    Dl2 = DFSetLabeled(Ds[ids], labels)
    Dl = vcat(Dl1, Dl2)
    Dl
end

function nmacs_vs_nonnmacs(Ds::DFSet)
    meta = getmeta(Ds)
    labels = map(b -> b ? 1 : 2, meta[:nmac]) #remap to 1=nmac,2=nonnmac
    labels = convert(Vector{Int64}, labels)
    Dl = DFSetLabeled(Ds, labels)
    Dl
end


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# ╔═╡ 21c50854-502a-11eb-3de4-7955c27ab0de
using LightGraphs, GraphPlot

# ╔═╡ df733764-5029-11eb-0068-7f149693d12c
md"""
# Trees

An extention of Linked list, such that a node(parent) can have multiple sub-nodes(children). A child node can have only one parent node. Start from a root node (Level 1) it branches out to multiple nodes (Level 2 and further to 3,4,...) and end on leafs (without any children). Height of node is the number of edges between itself and the root. Depth is inverse of height.

A tree must be completely conneced (there's always a path (group of edges) from root to the node) and there is no cycle.

Operations:
- Tree Traversal isn't as clear as in linked list
    - Depth-first Search (DFS): Move through one entire path to a leaf before moving to another path. If children are unexplored explore them first.
        - Pre-order Traversal: Check off the node as soon as you see it before checking its children. Check the left-most nodes by convention.
        - In-order Traversal: Check off the node when we see its left child and come back to it.
        - Post-order Traversal: Check off the node only after checking all its children (starting from the right or left, generally left by convention)
    - Breadth-first Search (BFS) or Level order: Priority visit all children on the same level before moving to child nodes (By Convention move from Left to right).

![Tree Traversal](https://algorithms.tutorialhorizon.com/files/2015/11/Tree-Traversals.png)

"""

# ╔═╡ 292cacf0-502a-11eb-03fd-735f45f94dc4
md"""
## Binary Trees

Nodes can have at most two children (0,1 or 2 children).

Since there is no order here we need to go through each node to search an element. Search is O(n)

Operation
- Deletion of a middle node is bit complicated it is also O(n)
- Insertion at a node need to travel down the path to find an empty spot
"""

# ╔═╡ 323e33d6-502a-11eb-304c-994320ec1518

# ╔═╡ Cell order:
# ╠═df733764-5029-11eb-0068-7f149693d12c
# ╠═21c50854-502a-11eb-3de4-7955c27ab0de
# ╠═292cacf0-502a-11eb-03fd-735f45f94dc4
# ╠═323e33d6-502a-11eb-304c-994320ec1518


text:
include("vi_optimizers.jl")
include("analytic.jl")
include("analyticVI.jl")
include("numericalVI.jl")
include("sampling.jl")
include("optimisers.jl")

export RobbinsMonro, ALRSVI

post_process!(::AbstractGP) = nothing

# Utils to iterate over inference objects
Base.length(::AbstractInference) = 1

Base.iterate(i::AbstractInference) = (i, nothing)
Base.iterate(::AbstractInference, ::Any) = nothing

isStochastic(i::AbstractInference) = i.stoch

## Multiple accessors
conv_crit(i::AbstractInference) = i.ϵ

## Conversion from natural to standard distribution parameters ##
function global_update!(gp::AbstractLatent)
    gp.post.Σ .= -0.5 * inv(nat2(gp))
    return gp.post.μ .= cov(gp) * nat1(gp)
end

## For the online case, the size may vary and inplace updates are note valid
function global_update!(gp::OnlineVarLatent)
    gp.post.Σ = -0.5 * inv(nat2(gp))
    return gp.post.μ = cov(gp) * nat1(gp)
end

## Default function for getting a view on y
xview(inf::AbstractInference) = inf.xview
setxview!(inf::AbstractInference, xview) = inf.xview = xview

setyview!(inf::AbstractInference, yview) = inf.yview = yview
yview(inf::AbstractInference) = inf.yview

nMinibatch(inf::AbstractInference) = inf.nMinibatch
setnMinibatch!(inf::AbstractInference, n::Int) = inf.nMinibatch = n

nSamples(i::AbstractInference) = i.nSamples
setnSamples!(inf::AbstractInference, n::Int) = inf.nSamples = n

getρ(inf::AbstractInference) = inf.ρ
setρ!(inf::AbstractInference, ρ) = inf.ρ = ρ

MBIndices(inf::AbstractInference) = inf.MBIndices
function setMBIndices!(inf::AbstractInference, mbindices::AbstractVector)
    return inf.MBIndices .= mbindices
end

setHPupdated!(inf::AbstractInference, status::Bool) = inf.HyperParametersUpdated = status
isHPupdated(inf::AbstractInference) = inf.HyperParametersUpdated

nIter(inf::AbstractInference) = inf.nIter

get_opt(::AbstractInference) = nothing
get_opt(i::VariationalInference) = i.vi_opt
get_opt(i::SamplingInference) = i.opt
get_opt(i::AbstractInference, n::Int) = get_opt(i)[n]
opt_type(i::AbstractInference) = first(get_opt(i))

# Initialize the final version of the inference objects
# using the right parametrization and size
function tuple_inference(
    i::AbstractInference,
    nLatent::Int,
    nFeatures::Int,
    nSamples::Int,
    nMinibatch::Int,
    xview::AbstractVector,
    yview::AbstractVector,
)
    return tuple_inference(
        i, nLatent, fill(nFeatures, nLatent), nSamples, nMinibatch, xview, yview
    )
end


text:
# The goal of this file is to compute the k_statistics polynomials as given by
# https://discourse.julialang.org/t/help-needed-translate-10-lines-from-maple-to-julia/60650
# and
# https://iris.unito.it/retrieve/handle/2318/1561388/399488/FASTpolykaysMAPLE.pdf
# also see
# Di Nardo E., G. Guarino, D. Senato (2007), A new method for fast computing unbiased estimators of cumulants. In press Statistics and Computing. http://www.springer.com/statistics/computational/journal/11222 (download from http://www.unibas.it/utenti/dinardo/lavori.html)

# We also compute their variances, but thoughr bootstrap, although there is probably a closed form for them.

function _better_stirling(n::T,k::T) where T<:Integer
# Algorithm taken directly from Combinatorics.jl, but we compute the whole array instead of only one by one.
arr = zeros(T,(n,k))
@inbounds for ni in 1:n
    @inbounds for ki in 1:min(ni,k)
        if ki ∈ (1,ni)
            arr[ni,ki] = 1
        elseif ki == ni - 1
            arr[ni,ki] = binomial(ni, 2)
        elseif ki == 2
            arr[ni,ki] = 2^(ni-1) - 1
        else
            @views arr[ni,ki] = ki * arr[ni - 1, ki] + arr[ni - 1, ki- 1]
        end
    end
end
return arr
end
#

_fd(k,j) = Polynomials.fromroots([i+j for i in 0:(k-j)])

function _compute_k(N,S)
m = BigInt(length(S))
result = []
stir = _better_stirling(m,m)
_facts = factorial.(0:m)
alternate__facts = (-1) .^ (0:m) .* _facts
u = Array{Polynomials.Polynomial{BigInt}}(undef,(m,))

for k in 1:m
    part = Combinatorics.integer_partitions(k)
    u[k] = Polynomials.Polynomial([0,(stir[k,1:k] .* alternate__facts[1:k])...])
    u2 = alternate__facts[1:k] .* _fd.(k-1,1:k)
    num = 0
    denom = prod(N .- (0:(k-1)))

    for i in eachindex(part)
        Cst = BigInt(_facts[k+1]//prod(_facts[v+1]^n * _facts[n+1] for (v,n) in StatsBase.countmap(part[i])))
        Ns = (prod(u[part[i]]).coeffs[2:end]'u2)(N)
        Ss = prod(S[part[i]])
        num += Cst * Ns * Ss
    end

    push!(result,num/denom)
end
return result
end
# # As a simple check, we can verify that the polynomials holds versus the theoretical ones.

# # More examples could be found in Fisher's work. 1929
struct FastPols
    N::Int64
    m::Int64
    pols::Vector{StaticPolynomials.Polynomial{BigFloat, SE, Nothing} where SE}
end

function FastPols(N,m)
    DynamicPolynomials.@polyvar S[1:m]
    pols = _compute_k(N,S)
    faster = StaticPolynomials.Polynomial.(pols)
    return FastPols(N,m,faster)
end

function _check_n_m!(fs::FastPols,N,m)
    if isnothing(fs)
        fs = FastPols(n,m)
    else
        @assert fs.N ==N
        @assert fs.m >= m
    end
end

# function k_statistics(data::Vector{T},t_star,m; k_stats_funs::FastPols) where T

#     # Empirical verxion from data
#     n = length(data)
#     ET_weights = exp.(t_star .* data) ./ n
#     if isnothing(k_stats_funs)
#         k_stats_funs = FastPols(n,m)
#     else
#         _check_n_m(k_stats_funs,n,m)
#     end
#     sw = sum(ET_weights)
#     zeroth_cumulant = log(sw)
#     ET_weights ./= sw
#     ET_S = vec(sum(data .^ (1:m)' .* ET_weights,dims=1))
#     ET_k_stats = [k_stats_funs.pols[i](ET_S)::T for i in 1:m]::Vector{T}
#     ET_k_stats ./= factorial.(0:(m-1))
#     return [zeroth_cumulant, ET_k_stats...]
# end

function k_statistics(data::Vector{T},t_star,m; k_stats_funs::FastPols) where T
    # Empirical verxion from data
    n = length(data)
    wts = exp.(t_star .* data) ./ n
    _check_n_m!(k_stats_funs,n,m)

    k0 = log(sum(wts))
    wts ./= sum(wts)
    S = vec(sum(data .^ (1:m)' .* wts,dims=1))*n
    k = [k_stats_funs.pols[i](S)::T for i in 1:m]::Vector{T}
    k ./= factorial.(0:(m-1))
    return [k0, k...]
end

function resampled_k_statistics(M,D::Vector{T},t_star,m; k_stats_funs::FastPols) where T
    # first step, resample the dataset D:
    n = length(D)
    _check_n_m!(k_stats_funs,n,m)

    ks = Array{T}(undef,(M,m+1))
    facts = factorial.(big.(0:(m-1)))
    arr = exp.(t_star .* D) .* (D.^(0:m)')
    arr[:,1] ./= n
    Threads.@threads for i in 1:M
        println(i)
        ks[i,:] = sum(arr[StatsBase.sample(1:n,n,replace=true),:], dims=1)
        ks[i,2:end] ./= ks[i,1] # Added line to match the non-bootstrapped k-staistics.
        ks[i,1] = log(ks[i,1])
        # the choice of the facts noramlisation highly condition the weights matrix.
        ks[i,2:end] = map(f->f(ks[i,2:end])::T,k_stats_funs.pols) ./ facts #or .* facts, which gave better fitting results, but is wrong...
    end
    return ks
end

function ugc_k_stats_old(α::Vector{T},θ::Vector{T},t_star,m) where T
    # We should make a test out of this one.
    s = θ ./ (1 .- t_star .* θ)
    return [sum(α .* log.(1 .-s)), [sum(α.*(s.^i)) for i in 1:m]...]
end
ugc_k_stats(D::ThorinDistributions.UnivariateGammaConvolution,t_star,m) = ugc_k_stats(D.α,D.θ,t_star,m)

function ugc_k_stats(α::Vector{T},θ::Vector{T},t_star,m) where T
    # Theoretical version for univariate gamma convolutions.
    s = θ ./ (1 .- t_star .* θ)
    s_pow = deepcopy(s)
    rez = Vector{T}(undef,(m+1,))
    rez[1] = sum(α .* log.(1 .-s))
    for i in 1:(m-1)
        rez[i+1] = α's_pow
        s_pow .*= s
    end
    rez[m+1] = α's_pow
    return rez
end

function ugc_k_stats!(rez,s_pow,s,α,θ,t_star,m)
    s .=  θ ./ (1 .- t_star .* θ)
    s_pow .= s
    rez[1] = sum(α .* log.(1 .-s))
    for i in 1:(m-1)
        rez[i+1] = α's_pow
        s_pow .*= s
    end
    rez[m+1] = α's_pow
end

# weird idea, but might work...
function correct_k_stats(k)
    lk = log.(k[2:end])
    ii = findfirst(lk[2:end] .> lk[1:end-1])
    a,b = [ones(ii) (1:ii)] \ lk[1:ii]
    lk[(ii+1):end] = a .+ b .* collect((ii+1):length(lk))
    return [k[1],exp.(lk)...]
end


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

"""
    Random

Support for generating random numbers. Provides [`rand`](@ref), [`randn`](@ref),
[`AbstractRNG`](@ref), [`MersenneTwister`](@ref), and [`RandomDevice`](@ref).
"""
module Random

include("DSFMT.jl")

using .DSFMT
using Base.GMP.MPZ
using Base.GMP: Limb
import SHA

using Base: BitInteger, BitInteger_types, BitUnsigned, require_one_based_indexing

import Base: copymutable, copy, copy!, ==, hash, convert,
             rand, randn, show

export rand!, randn!,
       randexp, randexp!,
       bitrand,
       randstring,
       randsubseq, randsubseq!,
       shuffle, shuffle!,
       randperm, randperm!,
       randcycle, randcycle!,
       AbstractRNG, MersenneTwister, RandomDevice, TaskLocalRNG, Xoshiro

## general definitions

"""
    AbstractRNG

Supertype for random number generators such as [`MersenneTwister`](@ref) and [`RandomDevice`](@ref).
"""
abstract type AbstractRNG end

Base.broadcastable(x::AbstractRNG) = Ref(x)

gentype(::Type{X}) where {X} = eltype(X)
gentype(x) = gentype(typeof(x))

### integers

# we define types which encode the generation of a specific number of bits
# the "raw" version means that the unused bits are not zeroed

abstract type UniformBits{T<:BitInteger} end

struct UInt10{T}    <: UniformBits{T} end
struct UInt10Raw{T} <: UniformBits{T} end

struct UInt23{T}    <: UniformBits{T} end
struct UInt23Raw{T} <: UniformBits{T} end

struct UInt52{T}    <: UniformBits{T} end
struct UInt52Raw{T} <: UniformBits{T} end

struct UInt104{T}    <: UniformBits{T} end
struct UInt104Raw{T} <: UniformBits{T} end

struct UInt2x52{T}    <: UniformBits{T} end
struct UInt2x52Raw{T} <: UniformBits{T} end

uint_sup(::Type{<:Union{UInt10,UInt10Raw}}) = UInt16
uint_sup(::Type{<:Union{UInt23,UInt23Raw}}) = UInt32
uint_sup(::Type{<:Union{UInt52,UInt52Raw}}) = UInt64
uint_sup(::Type{<:Union{UInt104,UInt104Raw}}) = UInt128
uint_sup(::Type{<:Union{UInt2x52,UInt2x52Raw}}) = UInt128

for UI = (:UInt10, :UInt10Raw, :UInt23, :UInt23Raw, :UInt52, :UInt52Raw,
          :UInt104, :UInt104Raw, :UInt2x52, :UInt2x52Raw)
    @eval begin
        $UI(::Type{T}=uint_sup($UI)) where {T} = $UI{T}()
        # useful for defining rand generically:
        uint_default(::$UI) = $UI{uint_sup($UI)}()
    end
end

gentype(::Type{<:UniformBits{T}}) where {T} = T

### floats

abstract type FloatInterval{T<:AbstractFloat} end

struct CloseOpen01{T<:AbstractFloat} <: FloatInterval{T} end # interval [0,1)
struct CloseOpen12{T<:AbstractFloat} <: FloatInterval{T} end # interval [1,2)

const FloatInterval_64 = FloatInterval{Float64}
const CloseOpen01_64   = CloseOpen01{Float64}
const CloseOpen12_64   = CloseOpen12{Float64}

CloseOpen01(::Type{T}=Float64) where {T<:AbstractFloat} = CloseOpen01{T}()
CloseOpen12(::Type{T}=Float64) where {T<:AbstractFloat} = CloseOpen12{T}()

gentype(::Type{<:FloatInterval{T}}) where {T<:AbstractFloat} = T

const BitFloatType = Union{Type{Float16},Type{Float32},Type{Float64}}

### Sampler

abstract type Sampler{E} end

gentype(::Type{<:Sampler{E}}) where {E} = E

# temporarily for BaseBenchmarks
RangeGenerator(x) = Sampler(default_rng(), x)

# In some cases, when only 1 random value is to be generated,
# the optimal sampler can be different than if multiple values
# have to be generated. Hence a `Repetition` parameter is used
# to choose the best one depending on the need.
const Repetition = Union{Val{1},Val{Inf}}

# these default fall-back for all RNGs would be nice,
# but generate difficult-to-solve ambiguities
# Sampler(::AbstractRNG, X, ::Val{Inf}) = Sampler(X)
# Sampler(::AbstractRNG, ::Type{X}, ::Val{Inf}) where {X} = Sampler(X)

"""
    Sampler(rng, x, repetition = Val(Inf))

Return a sampler object that can be used to generate random values from `rng` for `x`.

When `sp = Sampler(rng, x, repetition)`, `rand(rng, sp)` will be used to draw random values,
and should be defined accordingly.

`repetition` can be `Val(1)` or `Val(Inf)`, and should be used as a suggestion for deciding
the amount of precomputation, if applicable.

[`Random.SamplerType`](@ref) and [`Random.SamplerTrivial`](@ref) are default fallbacks for
*types* and *values*, respectively. [`Random.SamplerSimple`](@ref) can be used to store
pre-computed values without defining extra types for only this purpose.
"""
Sampler(rng::AbstractRNG, x, r::Repetition=Val(Inf)) = Sampler(typeof_rng(rng), x, r)
Sampler(rng::AbstractRNG, ::Type{X}, r::Repetition=Val(Inf)) where {X} =
    Sampler(typeof_rng(rng), X, r)

typeof_rng(rng::AbstractRNG) = typeof(rng)

# this method is necessary to prevent rand(rng::AbstractRNG, X) from
# recursively constructing nested Sampler types.
Sampler(T::Type{<:AbstractRNG}, sp::Sampler, r::Repetition) =
    throw(MethodError(Sampler, (T, sp, r)))

# default shortcut for the general case
Sampler(::Type{RNG}, X) where {RNG<:AbstractRNG} = Sampler(RNG, X, Val(Inf))
Sampler(::Type{RNG}, ::Type{X}) where {RNG<:AbstractRNG,X} = Sampler(RNG, X, Val(Inf))

#### pre-defined useful Sampler types

"""
    SamplerType{T}()

A sampler for types, containing no other information. The default fallback for `Sampler`
when called with types.
"""
struct SamplerType{T} <: Sampler{T} end

Sampler(::Type{<:AbstractRNG}, ::Type{T}, ::Repetition) where {T} = SamplerType{T}()

Base.getindex(::SamplerType{T}) where {T} = T

struct SamplerTrivial{T,E} <: Sampler{E}
    self::T
end

"""
    SamplerTrivial(x)

Create a sampler that just wraps the given value `x`. This is the default fall-back for
values.
The `eltype` of this sampler is equal to `eltype(x)`.

The recommended use case is sampling from values without precomputed data.
"""
SamplerTrivial(x::T) where {T} = SamplerTrivial{T,gentype(T)}(x)

Sampler(::Type{<:AbstractRNG}, x, ::Repetition) = SamplerTrivial(x)

Base.getindex(sp::SamplerTrivial) = sp.self

# simple sampler carrying data (which can be anything)
struct SamplerSimple{T,S,E} <: Sampler{E}
    self::T
    data::S
end

"""
    SamplerSimple(x, data)

Create a sampler that wraps the given value `x` and the `data`.
The `eltype` of this sampler is equal to `eltype(x)`.

The recommended use case is sampling from values with precomputed data.
"""
SamplerSimple(x::T, data::S) where {T,S} = SamplerSimple{T,S,gentype(T)}(x, data)

Base.getindex(sp::SamplerSimple) = sp.self

# simple sampler carrying a (type) tag T and data
struct SamplerTag{T,S,E} <: Sampler{E}
    data::S
    SamplerTag{T}(s::S) where {T,S} = new{T,S,gentype(T)}(s)
end

#### helper samplers

# TODO: make constraining constructors to enforce that those
# types are <: Sampler{T}

##### Adapter to generate a randome value in [0, n]

struct LessThan{T<:Integer,S} <: Sampler{T}
    sup::T
    s::S    # the scalar specification/sampler to feed to rand
end

function rand(rng::AbstractRNG, sp::LessThan)
    while true
        x = rand(rng, sp.s)
        x <= sp.sup && return x
    end
end

struct Masked{T<:Integer,S} <: Sampler{T}
    mask::T
    s::S
end

rand(rng::AbstractRNG, sp::Masked) = rand(rng, sp.s) & sp.mask

##### Uniform

struct UniformT{T} <: Sampler{T} end

uniform(::Type{T}) where {T} = UniformT{T}()

rand(rng::AbstractRNG, ::UniformT{T}) where {T} = rand(rng, T)

### machinery for generation with Sampler

# This describes how to generate random scalars or arrays, by generating a Sampler
# and calling rand on it (which should be defined in "generation.jl").
# NOTE: this section could be moved into a separate file when more containers are supported.

#### scalars

rand(rng::AbstractRNG, X)                                           = rand(rng, Sampler(rng, X, Val(1)))
# this is needed to disambiguate
rand(rng::AbstractRNG, X::Dims)                                     = rand(rng, Sampler(rng, X, Val(1)))
rand(rng::AbstractRNG=default_rng(), ::Type{X}=Float64) where {X} = rand(rng, Sampler(rng, X, Val(1)))::X

rand(X)                   = rand(default_rng(), X)
rand(::Type{X}) where {X} = rand(default_rng(), X)

#### arrays

rand!(A::AbstractArray{T}, X) where {T}             = rand!(default_rng(), A, X)
rand!(A::AbstractArray{T}, ::Type{X}=T) where {T,X} = rand!(default_rng(), A, X)

rand!(rng::AbstractRNG, A::AbstractArray{T}, X) where {T}             = rand!(rng, A, Sampler(rng, X))
rand!(rng::AbstractRNG, A::AbstractArray{T}, ::Type{X}=T) where {T,X} = rand!(rng, A, Sampler(rng, X))

function rand!(rng::AbstractRNG, A::AbstractArray{T}, sp::Sampler) where T
    for i in eachindex(A)
        @inbounds A[i] = rand(rng, sp)
    end
    A
end

rand(r::AbstractRNG, dims::Integer...) = rand(r, Float64, Dims(dims))
rand(                dims::Integer...) = rand(Float64, Dims(dims))

rand(r::AbstractRNG, X, dims::Dims)  = rand!(r, Array{gentype(X)}(undef, dims), X)
rand(                X, dims::Dims)  = rand(default_rng(), X, dims)

rand(r::AbstractRNG, X, d::Integer, dims::Integer...) = rand(r, X, Dims((d, dims...)))
rand(                X, d::Integer, dims::Integer...) = rand(X, Dims((d, dims...)))
# note: the above methods would trigger an ambiguity warning if d was not separated out:
# rand(r, ()) would match both this method and rand(r, dims::Dims)
# moreover, a call like rand(r, NotImplementedType()) would be an infinite loop

rand(r::AbstractRNG, ::Type{X}, dims::Dims) where {X} = rand!(r, Array{X}(undef, dims), X)
rand(                ::Type{X}, dims::Dims) where {X} = rand(default_rng(), X, dims)

rand(r::AbstractRNG, ::Type{X}, d::Integer, dims::Integer...) where {X} = rand(r, X, Dims((d, dims...)))
rand(                ::Type{X}, d::Integer, dims::Integer...) where {X} = rand(X, Dims((d, dims...)))

# SamplerUnion(X, Y, ...}) == Union{SamplerType{X}, SamplerType{Y}, ...}
SamplerUnion(U...) = Union{Any[SamplerType{T} for T in U]...}
const SamplerBoolBitInteger = SamplerUnion(Bool, BitInteger_types...)

include("Xoshiro.jl")
include("RNGs.jl")
include("generation.jl")
include("normal.jl")
include("misc.jl")
include("XoshiroSimd.jl")

## rand & rand! & seed! docstrings

"""
    rand([rng=default_rng()], [S], [dims...])

Pick a random element or array of random elements from the set of values specified by `S`;
`S` can be

* an indexable collection (for example `1:9` or `('x', "y", :z)`),
* an `AbstractDict` or `AbstractSet` object,
* a string (considered as a collection of characters), or
* a type: the set of values to pick from is then equivalent to `typemin(S):typemax(S)` for
  integers (this is not applicable to [`BigInt`](@ref)), to ``[0, 1)`` for floating
  point numbers and to ``[0, 1)+i[0, 1)`` for complex floating point numbers;

`S` defaults to [`Float64`](@ref).
When only one argument is passed besides the optional `rng` and is a `Tuple`, it is interpreted
as a collection of values (`S`) and not as `dims`.

!!! compat "Julia 1.1"
    Support for `S` as a tuple requires at least Julia 1.1.

# Examples
```julia-repl
julia> rand(Int, 2)
2-element Array{Int64,1}:
 1339893410598768192
 1575814717733606317

julia> using Random

julia> rand(MersenneTwister(0), Dict(1=>2, 3=>4))
1=>2

julia> rand((2, 3))
3

julia> rand(Float64, (2, 3))
2×3 Array{Float64,2}:
 0.999717  0.0143835  0.540787
 0.696556  0.783855   0.938235
```

!!! note
    The complexity of `rand(rng, s::Union{AbstractDict,AbstractSet})`
    is linear in the length of `s`, unless an optimized method with
    constant complexity is available, which is the case for `Dict`,
    `Set` and dense `BitSet`s. For more than a few calls, use `rand(rng,
    collect(s))` instead, or either `rand(rng, Dict(s))` or `rand(rng,
    Set(s))` as appropriate.
"""
rand

"""
    rand!([rng=default_rng()], A, [S=eltype(A)])

Populate the array `A` with random values. If `S` is specified
(`S` can be a type or a collection, cf. [`rand`](@ref) for details),
the values are picked randomly from `S`.
This is equivalent to `copyto!(A, rand(rng, S, size(A)))`
but without allocating a new array.

# Examples
```jldoctest
julia> rng = MersenneTwister(1234);

julia> rand!(rng, zeros(5))
5-element Vector{Float64}:
 0.5908446386657102
 0.7667970365022592
 0.5662374165061859
 0.4600853424625171
 0.7940257103317943
```
"""
rand!

"""
    seed!([rng=default_rng()], seed) -> rng
    seed!([rng=default_rng()]) -> rng

Reseed the random number generator: `rng` will give a reproducible
sequence of numbers if and only if a `seed` is provided. Some RNGs
don't accept a seed, like `RandomDevice`.
After the call to `seed!`, `rng` is equivalent to a newly created
object initialized with the same seed.

If `rng` is not specified, it defaults to seeding the state of the
shared task-local generator.

# Examples
```julia-repl
julia> Random.seed!(1234);

julia> x1 = rand(2)
2-element Vector{Float64}:
 0.32597672886359486
 0.5490511363155669

julia> Random.seed!(1234);

julia> x2 = rand(2)
2-element Vector{Float64}:
 0.32597672886359486
 0.5490511363155669

julia> x1 == x2
true

julia> rng = Xoshiro(1234); rand(rng, 2) == x1
true

julia> Xoshiro(1) == Random.seed!(rng, 1)
true

julia> rand(Random.seed!(rng), Bool) # not reproducible
true

julia> rand(Random.seed!(rng), Bool) # not reproducible either
false

julia> rand(Xoshiro(), Bool) # not reproducible either
true
```
"""
seed!(rng::AbstractRNG, ::Nothing) = seed!(rng)

end # module

max_stars_repo_name: Physolia/julia

max_issues_repo_name: Eurus-Holmes/julia

max_forks_repo_name: Nexuscompute/julia


text:
using DelimitedFiles
using WriteVTK
using MeshCore: AbsShapeDesc, P1, L2, T3, Q4, T4, H8, T6
using MeshCore: VecAttrib, @_check, datavaluetype
using MeshCore: shapedesc, nshapes, IncRel
using MeshCore: ShapeColl, ir_skeleton, ir_bbyfacets, nshapes
using LinearAlgebra: norm

"""
    T3blockx(xs::Vector{T}, ys::Vector{T}, orientation::Symbol; intbytes = 8) where {T}

Generate a graded triangle mesh  of a 2D block.

The mesh is produced by splitting each logical  rectangular cell into two
triangles. Orientation may be chosen as `:a`, `:b`.

Keyword argument `intbytes` controls the size of the integer indexes.

# Return
By convention the function returns an incidence relation (`connectivity`)
between the elements (`connectivity.left`) and the vertices
(`connectivity.right`). The geometry is stored as the attribute "geom" of the
vertices.
"""
function T3blockx(xs::Vector{T}, ys::Vector{T}, orientation::Symbol; intbytes = 8) where {T}
    if (orientation==:a) || (orientation==:b)
        # nothing
    else
        error("Unknown orientation")
    end
    inttype = Int64
    if intbytes == 4
        inttype = Int32
    end
    nL = length(xs)-1;
    nW = length(ys)-1;
    nnodes = (nL+1)*(nW+1);
    ncells = 2*(nL)*(nW);
    xys = zeros(T, nnodes, 2);
    conns = zeros(inttype, ncells, 3);

    f=1;
    for j in 1:(nW+1)
        for i in 1:(nL+1)
            xys[f, 1] = xs[i]
            xys[f, 2] = ys[j]
            f=f+1;
        end
    end

    gc=1;
    for i in 1:nL
        for j in 1:nW
            f=(j-1)*(nL+1)+i;
            if     (orientation==:a)
                conns[gc,:] .= f, (f+1), f+(nL+1)
            elseif (orientation==:b)
                conns[gc,:] .= f, (f+1), f+(nL+1)+1
            end
            gc=gc+1;
            if     (orientation==:a)
                conns[gc,:] .= (f+1), f+(nL+1)+1, f+(nL+1)
            elseif (orientation==:b)
                conns[gc,:] .= f, f+(nL+1)+1, f+(nL+1)
            end
            gc=gc+1;
        end
    end

    C = inttype.(conns[1:gc-1, :]);
    N, TC = size(xys, 2), eltype(xys)
    locs =  VecAttrib([SVector{N, TC}(xys[i, :]) for i in 1:size(xys, 1)])
    vertices = ShapeColl(P1, length(locs), "vertices")
    vertices.attributes["geom"] = locs
    elements = ShapeColl(T3, size(C, 1), "elements")
    return IncRel(elements, vertices, C)
end

"""
    T3block(Length, Width, nL, nW, orientation = :a)

Generate a triangle mesh  of a 2D block.

See also: T3blockx
"""
function T3block(Length, Width, nL, nW, orientation = :a; kwargs...)
    return T3blockx(collect(linearspace(0.0, Length, nL+1)),
                    collect(linearspace(0.0, Width, nW+1)), orientation; kwargs...);
end

"""
    T3toT6(ir)

Convert three node triangles (T3) to six node triangles (T6).
"""
function T3toT6(ir)
    @_check shapedesc(ir.left) == T3
    locs = ir.right.attributes["geom"]
    sk = ir_skeleton(ir)  # skeleton consists of edges
    bf = ir_bbyfacets(ir, sk) # edges bounding triangles
    n = nshapes(ir.right) # number of vertices in the three-node triangle mesh
    med = [idx+n for idx in 1:nshapes(sk.left)]
    C = fill(zero(indextype(ir)), nshapes(ir.left), 6)
    # Create array to hold the new coordinates
    nx = fill(zero(eltype(datavaluetype(locs))), length(locs) + length(med), length(locs[1]))
    for i in 1:length(locs) # copy in the old locations
        nx[i, :] .= locs[i]
    end
    for i in 1:nshapes(ir.left)
        C[i, 1:3] .= ir[i]
        ec = abs.(bf[i]) # ignore the orientation
        for (j, w) in enumerate((5, 6, 4))
            en = med[ec[j]]
            C[i, w] = en
            ev = sk[ec[j]]
            nx[en, :] = (locs[ev[1]] + locs[ev[2]]) / 2.0
        end
    end
    locs =  VecAttrib([SVector{size(nx, 2), eltype(locs[1])}(nx[i, :]) for i in 1:size(nx, 1)])
    vertices = ShapeColl(P1, size(nx, 1), "vertices")
    vertices.attributes["geom"] = locs
    elements = ShapeColl(T6, size(C, 1), "elements")
    return IncRel(elements, vertices, C)
end

"""
    T6blockx(xs::Vector{T}, ys::Vector{T}, orientation::Symbol; intbytes = 8) where
        {T}

Generate a graded quadratic triangle mesh  of a 2D block.

The mesh is produced by splitting each logical  rectangular cell into two
triangles. Orientation may be chosen as `:a`, `:b`.

Keyword argument `intbytes` controls the size of the integer indexes.

# Return
By convention the function returns an incidence relation (`connectivity`)
between the elements (`connectivity.left`) and the vertices
(`connectivity.right`). The geometry is stored as the attribute "geom" of the
vertices.
"""
function T6blockx(xs::Vector{T}, ys::Vector{T}, orientation::Symbol; intbytes = 8) where
    {T}
    ir = T3blockx(xs, ys, orientation; intbytes = intbytes)
    return T3toT6(ir)
end

"""
    T6block(Length, Width, nL, nW, orientation = :a)

Generate a quadratic triangle mesh of a 2D block.

See also: T6blockx
"""
function T6block(Length, Width, nL, nW, orientation = :a; kwargs...)
    return T6blockx(collect(linearspace(0.0, Length, nL+1)),
                    collect(linearspace(0.0, Width, nW+1)), orientation; kwargs...);
end

"""
    T6toT3(ir)

Convert six node triangles (T6) to three node triangles (T3).

Note: Unconnected vertices will be removed from the right-hand side shape
collection.
"""
function T6toT3(ir)
    @_check shapedesc(ir.left) == T6
    # The same vertices
    vertices = ShapeColl(P1, nshapes(ir.right), "vertices")
    locs = ir.right.attributes["geom"]
    vertices.attributes["geom"] = locs
    # Elements have only the three corner nodes
    C = [SVector{3}(ir[idx][1:3]) for idx in 1:nrelations(ir)]
    elements = ShapeColl(T3, length(C), "elements")
    newir = IncRel(elements, vertices, C)
    # Remove unconnected vertices.
    ic = vconnected(newir)
    nn = vnewnumbering(newir, ic)
    return compactify(newir, nn)
end


text:
using RoguePkg
using Compat
using Compat.Test

mktempdir() do dir
    pkg_dir = joinpath(dir, "Foo")
    src_dir = joinpath(pkg_dir, "src")
    test_dir = joinpath(pkg_dir, "test")
    deps_dir = joinpath(pkg_dir, "deps")
    mkpath(src_dir)
    mkpath(test_dir)
    mkpath(deps_dir)
    write(joinpath(src_dir, "Foo.jl"), "module Foo\nend")
    # this will test Pkg.Entry.resolve
    write(joinpath(test_dir, "REQUIRE"), "Parameters")
    # running tests will create an empty file, which we can check
    test_path = tempname()
    deps_path = tempname()
    write(joinpath(test_dir, "runtests.jl"), "touch(\"$(test_path)\")")
    write(joinpath(deps_dir, "build.jl"), "touch(\"$(deps_path)\")")

    function all_tests(pkg)
        @test Pkg.dir(pkg, "src") == src_dir
        @test Pkg.dir("something", pkg) == Pkg.dir("something", "Foo") # NOTE hardcoded
        Pkg.build(pkg)
        @test isfile(deps_path)
        rm(deps_path; force = true) # may not exist if test doesn't work, that should not error
        @test (Pkg.test(pkg); true)
        @test isfile(test_path)
        rm(test_path; force = true) # may not exist if test doesn't work, that should not error
    end

    # find by path
    pkg = RoguePkg.PkgAtPath(@eval @pkg_at_str($pkg_dir))
    all_tests(pkg)
    @test repr(pkg) == "package at path $pkg_dir"
    # find by module
    push!(LOAD_PATH, dir)
    pkg = pkg_for"Foo"
    @test repr(pkg) == "package for module \"Foo\" (located dynamically)"
    all_tests(pkg)
    # find by path in environment
    @test repr(pkg_local"Foo") ==
        "package \"Foo\" in ENV[\"JULIA_LOCAL_PACKAGES\"] \e[31m(not set)\e[39m"
    withenv("JULIA_LOCAL_PACKAGES" => dir) do
        pkg = pkg_local"Foo"
        @test repr(pkg) == "package \"Foo\" in $(dir)"
        all_tests(pkg)
    end
    @test_throws ErrorException Pkg.dir(pkg_for"this_better_not_exist")
end


text:
@testset "Get All Symbols Function Tests" begin

  @test isdefined(Julz, :get_all_symbols) == true

end


text:
## consumer

mutable struct PartitionList
    rkparlist::Ptr{CKafkaTopicPartitionList}
end

function PartitionList()
    ptr = kafka_topic_partition_list_new()
    parlist = PartitionList(ptr)
    finalizer(parlist -> kafka_topic_partition_list_destroy(ptr), parlist)
    return parlist
end

mutable struct KafkaConsumer
    client::KafkaClient
    parlist::PartitionList
end

function KafkaConsumer(conf::Dict; dr_cb=nothing, err_cb=nothing)
    @assert haskey(conf, "bootstrap.servers") "`bootstrap.servers` should be specified in conf"
    @assert haskey(conf, "group.id") "`group.id` should be specified in conf"
    client = KafkaClient(KAFKA_TYPE_CONSUMER, conf; dr_cb=dr_cb, err_cb=err_cb)
    parlist = PartitionList()
    return KafkaConsumer(client, parlist)
end

function KafkaConsumer(bootstrap_servers::String, group_id::String, conf::Dict=Dict(); dr_cb=nothing, err_cb=nothing)
    conf["bootstrap.servers"] = bootstrap_servers
    conf["group.id"] = group_id
    return KafkaConsumer(conf; dr_cb=dr_cb, err_cb=err_cb)
end

function Base.show(io::IO, c::KafkaConsumer)
    group_id = c.client.conf["group.id"]
    bootstrap_servers = c.client.conf["bootstrap.servers"]
    print(io, "KafkaConsumer($group_id @ $bootstrap_servers)")
end

function subscribe(c::KafkaConsumer, topic_partitions::Vector{Tuple{String, Int}})
    rkparlist = c.parlist.rkparlist
    for (topic, partition) in topic_partitions
        kafka_topic_partition_list_add(rkparlist, topic, partition)
    end
    kafka_subscribe(c.client.rk, rkparlist)
end

function offsets_for_timestamp(c::KafkaConsumer, timestamp::Int64, timeout::Int=1000)
    rkparlist = c.parlist.rkparlist
    kafka_assignment(c.client.rk, rkparlist)
    tpl = unsafe_load(rkparlist)::CKafkaTopicPartitionList
    for i in 1:tpl.cnt
        e = unsafe_load(tpl.elems, i)
        e.offset = timestamp
        unsafe_store!(tpl.elems, e, i)
    end
    kafka_offsets_for_times(c.client.rk, rkparlist, timeout)
end

function seek(c::KafkaConsumer, timestamp::Int64, timeout::Int=1000)
    offsets_for_timestamp(c, timestamp, timeout)
    rkparlist = c.parlist.rkparlist
    kafka_assign(c.client.rk, rkparlist)
    tpl = unsafe_load(rkparlist)::CKafkaTopicPartitionList
    for i in 1:tpl.cnt
        e = unsafe_load(tpl.elems, i)
        topic = unsafe_string(e.topic)
        kt = KafkaTopic(c.client, topic, Dict())
        kafka_seek(kt.rkt, e.partition, e.offset, timeout)
    end
end

function poll(::Type{K}, ::Type{P}, c::KafkaConsumer, timeout::Int=1000) where {K,P}
    c_msg_ptr = kafka_consumer_poll(c.client.rk, timeout)
    if c_msg_ptr != nothing
        c_msg = unsafe_load(c_msg_ptr)
        msg = Message{K,P}(c_msg)
        kafka_message_destroy(c_msg_ptr)
        return msg
    else
        return nothing
    end
end

poll(c::KafkaConsumer, timeout::Int=1000) = poll(Vector{UInt8}, Vector{UInt8}, c, timeout)


text:
function error_correction_solution_non_iterative()

end

function error_correction_solution_iterative()

end

"""
To Do:
- fix bug in Autumn compiler where default behavior is not applied to list elements when only some elements of the list have on-clauses applied
- mapping/parsing fixes:
  - DONE bias later frame decompositions with object types from earlier, to avoid issues like the Ice-parsing problem
  - DONE switch between singlecell and non-singlecell parsing versions when observation trace has overlaps
  - DONE support object types having a color field, at least in the multi-cellular case
  - POSTPONE supporting rotations for object types, i.e. not defining new object types for shapes that have just been rotated
- dynamics matrix parsing
  - generating event predicates via genBool
  - augmenting dynamics matrix with abstractions of atomic update rules (i.e. abstracting concrete function parameters via gen[Position, Int, Bool])
  - implementing non-iterative and iterative parsing algorithms:
    - non-iterative approach
      - eliminating elements of dynamics matrix cells that are redundant
      - finding the simplest object to handle inside each type (i.e. fewest possibilities, most single-update rule cells)
      - minimizing the number of update function changes (dynamic programming?)
    - iterative approach
      -
"""