text:
#-------------------------------------------------------------------------------
# Flags hold auxilary information about tokens/nonterminals which the Kind
# doesn't capture in a nice way.
const RawFlags = UInt32
const EMPTY_FLAGS = RawFlags(0)
const TRIVIA_FLAG = RawFlags(1<<0)
# Some of the following flags are head-specific and could probably be allowed
# to cover the same bits...
const INFIX_FLAG  = RawFlags(1<<1)
# Record whether syntactic operators were dotted
const DOTOP_FLAG = RawFlags(1<<2)
# Set when kind == K"String" was triple-delimited as with """ or ```
const TRIPLE_STRING_FLAG = RawFlags(1<<3)
# Set when a string or identifier needs "raw string" unescaping
const RAW_STRING_FLAG = RawFlags(1<<4)
# try-finally-catch
const TRY_CATCH_AFTER_FINALLY_FLAG = RawFlags(1<<5)
# Flags holding the dimension of an nrow or other UInt8 not held in the source
const NUMERIC_FLAGS = RawFlags(RawFlags(0xff)<<8)
# Todo ERROR_FLAG = 0x80000000 ?

function set_numeric_flags(n::Integer)
    f = RawFlags((n << 8) & NUMERIC_FLAGS)
    if numeric_flags(f) != n
        error("Numeric flags unable to hold large integer $n")
    end
    f
end

function numeric_flags(f::RawFlags)
    Int((f >> 8) % UInt8)
end

# Return true if any of `test_flags` are set
has_flags(flags::RawFlags, test_flags) = (flags & test_flags) != 0

#-------------------------------------------------------------------------------
struct SyntaxHead
    kind::Kind
    flags::RawFlags
end

kind(head::SyntaxHead) = head.kind
flags(head::SyntaxHead) = head.flags
has_flags(head::SyntaxHead, test_flags) = has_flags(flags(head), test_flags)

is_trivia(head::SyntaxHead) = has_flags(head, TRIVIA_FLAG)
is_infix(head::SyntaxHead)  = has_flags(head, INFIX_FLAG)
is_dotted(head::SyntaxHead) = has_flags(head, DOTOP_FLAG)
numeric_flags(head::SyntaxHead) = numeric_flags(flags(head))
is_error(head::SyntaxHead)  = kind(head) == K"error"

function Base.summary(head::SyntaxHead)
    untokenize(head, unique=false, include_flag_suff=false)
end

function untokenize(head::SyntaxHead; unique=true, include_flag_suff=true)
    str = untokenize(kind(head); unique=unique)
    if is_dotted(head)
        str = "."*str
    end
    if include_flag_suff && flags(head) ∉ (EMPTY_FLAGS, DOTOP_FLAG)
        str = str*"-"
        is_trivia(head)  && (str = str*"t")
        is_infix(head)   && (str = str*"i")
        has_flags(head, TRIPLE_STRING_FLAG) && (str = str*"s")
        has_flags(head, RAW_STRING_FLAG) && (str = str*"r")
        has_flags(head, TRY_CATCH_AFTER_FINALLY_FLAG) && (str = str*"f")
        n = numeric_flags(head)
        n != 0 && (str = str*string(n))
    end
    str
end

#-------------------------------------------------------------------------------
"""
`SyntaxToken` is a token covering a contiguous byte range in the input text.
Information about preceding whitespace is added for use by the parser.
"""
struct SyntaxToken
    kind::Kind
    first_byte::UInt32
    last_byte::UInt32
    # Flags for leading whitespace
    is_dotted::Bool
    is_suffixed::Bool
    had_whitespace::Bool
    had_newline::Bool
end

function SyntaxToken(raw::RawToken, had_whitespace, had_newline)
    SyntaxToken(raw.kind, raw.startbyte + 1, raw.endbyte + 1, raw.dotop, raw.suffix,
                had_whitespace, had_newline)
end

function Base.show(io::IO, tok::SyntaxToken)
    range = string(lpad(first_byte(tok), 3), ":", rpad(last_byte(tok), 3))
    print(io, rpad(range, 17, " "), rpad(kind(tok), 15, " "))
end

kind(tok::SyntaxToken) = tok.kind
flags(tok::SyntaxToken) = tok.is_dotted ? DOTOP_FLAG : EMPTY_FLAGS
first_byte(tok::SyntaxToken) = tok.first_byte
last_byte(tok::SyntaxToken) = tok.last_byte
span(tok::SyntaxToken) = last_byte(tok) - first_byte(tok) + 1

is_dotted(tok::SyntaxToken)    = tok.is_dotted
is_suffixed(tok::SyntaxToken)  = tok.is_suffixed
is_decorated(tok::SyntaxToken) = is_dotted(tok) || is_suffixed(tok)

Base.:(==)(tok::SyntaxToken, k::Kind) = (kind(tok) == k && !is_decorated(tok))

#-------------------------------------------------------------------------------

"""
Range in the source text which will become a node in the tree. Can be either a
token (leaf node of the tree) or an interior node, depending on how the
start_mark compares to previous nodes.

TODO: Optimize this data structure?  It's very large at the moment.
"""
struct TaggedRange
    head::SyntaxHead # Kind,flags
    orig_kind::Kind  # Kind of the original token for leaf tokens, or K"None"
    first_byte::UInt32  # First byte in the input text
    last_byte::UInt32   # Last byte in the input text
    start_mark::UInt32  # Index of first emitted range which this range covers
end

head(range::TaggedRange)       = range.head
kind(range::TaggedRange)       = kind(range.head)
flags(range::TaggedRange)      = flags(range.head)
first_byte(range::TaggedRange) = Int(range.first_byte)
last_byte(range::TaggedRange)  = Int(range.last_byte)
span(range::TaggedRange)       = 1 + last_byte(range) - first_byte(range)

#-------------------------------------------------------------------------------
struct ParseStreamPosition
    input_byte::Int    # Index of next byte in input
    output_index::Int  # Index of last span in output
end

const NO_POSITION = ParseStreamPosition(0,0)

#-------------------------------------------------------------------------------
"""
ParseStream provides an IO interface for the parser. It
- Wraps the lexer with a lookahead buffer
- Removes insignificant whitespace and comment tokens, shifting them into the
  output implicitly (newlines may be significant depending on `skip_newlines`)
"""
mutable struct ParseStream
    # `textbuf` is a buffer of UTF-8 encoded text of the source code. This is a
    # natural representation as we desire random access and zero-copy parsing
    # of UTF-8 text from various containers, and unsafe_wrap(Vector{UInt8},
    # ...) allows us to use a Vector here.
    #
    # We want `ParseStream` to be concrete so that all `parse_*` functions only
    # need to be compiled once. Thus `textbuf` must not be parameterized here.
    textbuf::Vector{UInt8}
    # GC root for the object which owns the memory in `textbuf`. `nothing` if
    # the `textbuf` owner was unknown (eg, ptr,length was passed)
    text_root::Any
    # Lexer, transforming the input bytes into a token stream
    lexer::Tokenize.Lexers.Lexer{IOBuffer,RawToken}
    # Lookahead buffer for already lexed tokens
    lookahead::Vector{SyntaxToken}
    # Parser output as an ordered sequence of ranges, parent nodes after children.
    ranges::Vector{TaggedRange}
    # Parsing diagnostics (errors/warnings etc)
    diagnostics::Vector{Diagnostic}
    # First byte of next token
    next_byte::Int
    # Counter for number of peek()s we've done without making progress via a bump()
    peek_count::Int
    # (major,minor) version of Julia we're parsing this code for.
    # May be different from VERSION!
    version::Tuple{Int,Int}

    function ParseStream(text_buf::Vector{UInt8}, text_root, next_byte::Integer,
                         version::VersionNumber)
        io = IOBuffer(text_buf)
        seek(io, next_byte-1)
        lexer = Tokenize.Lexers.Lexer(io, RawToken)
        # To avoid keeping track of the exact Julia development version where new
        # features were added or comparing prerelease strings, we treat prereleases
        # or dev versons as the release version using only major and minor version
        # numbers. This means we're inexact for old dev versions but that seems
        # like an acceptable tradeoff.
        ver = (version.major, version.minor)
        new(text_buf, text_root, lexer,
            Vector{SyntaxToken}(),
            Vector{TaggedRange}(),
            Vector{Diagnostic}(),
            next_byte,
            0,
            ver)
    end
end

function ParseStream(text::Vector{UInt8}, index::Integer=1; version=VERSION)
    ParseStream(text, text, index, version)
end

# Buffer with unknown owner. Not exactly recommended, but good for C interop
function ParseStream(ptr::Ptr{UInt8}, len::Integer, index::Integer=1; version=VERSION)
    ParseStream(unsafe_wrap(Vector{UInt8}, ptr, len), nothing, index, version)
end

# Buffers originating from strings
function ParseStream(text::String, index::Integer=1; version=VERSION)
    ParseStream(unsafe_wrap(Vector{UInt8}, text),
                text, index, version)
end
function ParseStream(text::SubString, index::Integer=1; version=VERSION)
    # See also IOBuffer(SubString("x"))
    ParseStream(unsafe_wrap(Vector{UInt8}, pointer(text), sizeof(text)),
                text, index, version)
end
function ParseStream(text::AbstractString, index::Integer=1; version=VERSION)
    ParseStream(String(text), index; version=version)
end

# IO-based cases
function ParseStream(io::IOBuffer; version=VERSION)
    ParseStream(io.data, io, position(io)+1, version)
end
function ParseStream(io::Base.GenericIOBuffer; version=VERSION)
    textbuf = unsafe_wrap(Vector{UInt8}, pointer(io.data), length(io.data))
    ParseStream(textbuf, io, position(io)+1, version)
end
function ParseStream(io::IOStream; version=VERSION)
    textbuf = Mmap.mmap(io)
    ParseStream(textbuf, io, position(io)+1, version)
end
function ParseStream(io::IO; version=VERSION)
    textbuf = read(io)
    ParseStream(textbuf, textbuf, 1, version)
end

function Base.show(io::IO, mime::MIME"text/plain", stream::ParseStream)
    println(io, "ParseStream at position $(stream.next_byte)")
end

function show_diagnostics(io::IO, stream::ParseStream, code)
    show_diagnostics(io, stream.diagnostics, code)
end

#-------------------------------------------------------------------------------
# Stream input interface - the peek_* family of functions

# Buffer up until the next non-whitespace token.
# This can buffer more than strictly necessary when newlines are significant,
# but this is not a big problem.
function _buffer_lookahead_tokens(stream::ParseStream)
    had_whitespace = false
    had_newline    = false
    while true
        raw = Tokenize.Lexers.next_token(stream.lexer)
        k = TzTokens.exactkind(raw)
        was_whitespace = k in (K"Whitespace", K"Comment", K"NewlineWs")
        was_newline    = k == K"NewlineWs"
        had_whitespace |= was_whitespace
        had_newline    |= was_newline
        push!(stream.lookahead, SyntaxToken(raw, had_whitespace, had_newline))
        if !was_whitespace
            break
        end
    end
end

# Find the index of the first nontrivia token in the lookahead buffer.
#
# TODO: Store this as part of _buffer_lookahead_tokens to avoid redoing this
# work all the time!
function _lookahead_index(stream::ParseStream, n::Integer, skip_newlines::Bool)
    i = 1
    while true
        if i > length(stream.lookahead)
            _buffer_lookahead_tokens(stream)
        end
        k = kind(stream.lookahead[i])
        is_skipped =  k ∈ (K"Whitespace", K"Comment") ||
                     (k == K"NewlineWs" && skip_newlines)
        if !is_skipped
            if n == 1
                return i
            end
            n -= 1
        end
        i += 1
    end
end

"""
    peek(stream [, n=1])

Look ahead in the stream `n` tokens, returning the token kind. Comments and
non-newline whitespace are skipped automatically. Whitespace containing a
single newline is returned as kind `K"NewlineWs"` unless `skip_newlines` is
true.
"""
function Base.peek(stream::ParseStream, n::Integer=1; skip_newlines::Bool=false)
    kind(peek_token(stream, n; skip_newlines=skip_newlines))
end

"""
    peek_token(stream [, n=1])

Like `peek`, but return the full token information rather than just the kind.
"""
function peek_token(stream::ParseStream, n::Integer=1; skip_newlines=false)
    stream.peek_count += 1
    if stream.peek_count > 100_000
        error("The parser seems stuck at byte $(stream.next_byte)")
    end
    stream.lookahead[_lookahead_index(stream, n, skip_newlines)]
end

function _peek_behind_fields(ranges, i)
    r = ranges[i]
    return (kind=kind(r),
            flags=flags(r),
            orig_kind=r.orig_kind,
            is_leaf=r.start_mark == i)
end

"""
    peek_behind(ps; skip_trivia=true)
    peek_behind(ps, pos::ParseStreamPosition)

Return information about a span which was previously inserted into the output,
defaulting to the most previous nontrivia node when `skip_trivia` is true, or
at the provided position `pos`.

Retroactively inspecting or modifying the parser's output can be confusing, so
using this function should be avoided where possible.
"""
function peek_behind(stream::ParseStream; skip_trivia::Bool=true)
    if skip_trivia
        for i = length(stream.ranges):-1:1
            r = stream.ranges[i]
            if !is_trivia(head(r)) && kind(r) != K"TOMBSTONE"
                return _peek_behind_fields(stream.ranges, i)
            end
        end
    elseif !isempty(stream.ranges)
        return _peek_behind_fields(stream.ranges, lastindex(stream.ranges))
    end
    internal_error("Can't peek behind at start of stream")
end

function peek_behind(stream::ParseStream, pos::ParseStreamPosition)
    return _peek_behind_fields(stream.ranges, pos.output_index)
end

#-------------------------------------------------------------------------------
# Stream output interface - the `bump_*` and `emit_*` family of functions
#
# Though note bump() really does both input and output

# Bump the next `n` tokens
# flags and remap_kind are applied to any non-trivia tokens
function _bump_n(stream::ParseStream, n::Integer, flags, remap_kind=K"None")
    if n <= 0
        return
    end
    for i = 1:n
        tok = stream.lookahead[i]
        k = kind(tok)
        if k == K"EndMarker"
            break
        end
        is_trivia = k ∈ (K"Whitespace", K"Comment", K"NewlineWs")
        f = is_trivia ? TRIVIA_FLAG : flags
        is_dotted(tok) && (f |= DOTOP_FLAG)
        outk = (is_trivia || remap_kind == K"None") ? k : remap_kind
        range = TaggedRange(SyntaxHead(outk, f), k, first_byte(tok),
                            last_byte(tok), lastindex(stream.ranges)+1)
        push!(stream.ranges, range)
    end
    Base._deletebeg!(stream.lookahead, n)
    stream.next_byte = last_byte(last(stream.ranges)) + 1
    # Defuse the time bomb
    stream.peek_count = 0
end

"""
    bump(stream [, flags=EMPTY_FLAGS];
         skip_newlines=false, error, remap_kind)

Shift the current token from the input to the output, adding the given flags.
"""
function bump(stream::ParseStream, flags=EMPTY_FLAGS; skip_newlines=false,
              error=nothing, remap_kind::Kind=K"None")
    emark = position(stream)
    _bump_n(stream, _lookahead_index(stream, 1, skip_newlines), flags, remap_kind)
    if !isnothing(error)
        emit(stream, emark, K"error", flags, error=error)
    end
    # Return last token location in output if needed for reset_node!
    return position(stream)
end

"""
Bump comments and whitespace tokens preceding the next token

**Skips newlines** by default.  Set skip_newlines=false to avoid that.
"""
function bump_trivia(stream::ParseStream, flags=EMPTY_FLAGS;
                     skip_newlines=true, error=nothing)
    emark = position(stream)
    _bump_n(stream, _lookahead_index(stream, 1, skip_newlines) - 1, EMPTY_FLAGS)
    if !isnothing(error)
        emit(stream, emark, K"error", flags, error=error)
    end
    return position(stream)
end

"""
Bump an invisible zero-width token into the output

This is useful when surrounding syntax implies the presence of a token.  For
example, `2x` means `2*x` via the juxtoposition rules.
"""
function bump_invisible(stream::ParseStream, kind, flags=EMPTY_FLAGS;
                        error=nothing)
    emit(stream, position(stream), kind, flags, error=error)
    return position(stream)
end

"""
Bump several tokens, gluing them together into a single token

This is for use in special circumstances where the parser needs to resolve
lexing ambiguities. There's no special whitespace handling — bump any
whitespace if necessary with bump_trivia.
"""
function bump_glue(stream::ParseStream, kind, flags, num_tokens)
    span = TaggedRange(SyntaxHead(kind, flags), K"None",
                       first_byte(stream.lookahead[1]),
                       last_byte(stream.lookahead[num_tokens]),
                       lastindex(stream.ranges) + 1)
    Base._deletebeg!(stream.lookahead, num_tokens)
    push!(stream.ranges, span)
    stream.next_byte = last_byte(last(stream.ranges)) + 1
    stream.peek_count = 0
    return position(stream)
end

"""
    bump_split(stream, token_spec1, [token_spec2 ...])

Bump the next token, splitting it into several pieces

Tokens are defined by a number of `token_spec` of shape `(nbyte, kind, flags)`.
The number of input bytes of the last spec is taken from the remaining bytes of
the input token, with the associated `nbyte` ignored.

This is a hack which helps resolves the occasional lexing ambiguity. For
example
* Whether .+ should be a single token or the composite (. +) which is used for
  standalone operators.
* Whether ... is splatting (most of the time) or three . tokens in import paths

TODO: Are these the only cases?  Can we replace this general utility with a
simpler one which only splits preceding dots?
"""
function bump_split(stream::ParseStream, split_spec...)
    tok = popfirst!(stream.lookahead)
    fbyte = first_byte(tok)
    for (i, (nbyte, k, f)) in enumerate(split_spec)
        lbyte = (i == length(split_spec)) ? last_byte(tok) : fbyte + nbyte - 1
        push!(stream.ranges, TaggedRange(SyntaxHead(k, f), kind(tok),
                                         fbyte, lbyte,
                                         lastindex(stream.ranges) + 1))
        fbyte += nbyte
    end
    stream.next_byte = last_byte(last(stream.ranges)) + 1
    stream.peek_count = 0
    # Returning position(stream) like the other bump* methods would be
    # ambiguous here; return nothing instead.
    nothing
end

"""
Reset kind or flags of an existing node in the output stream

This is a hack, but in some limited occasions the trailing syntax may change
the kind or flags of a token in a way which would require unbounded lookahead
in a recursive descent parser. Modifying the output with reset_node! is useful
in those cases.
"""
function reset_node!(stream::ParseStream, pos::ParseStreamPosition;
                     kind=nothing, flags=nothing)
    range = stream.ranges[pos.output_index]
    k = isnothing(kind)  ? (@__MODULE__).kind(range)  : kind
    f = isnothing(flags) ? (@__MODULE__).flags(range) : flags
    stream.ranges[pos.output_index] =
        TaggedRange(SyntaxHead(k, f), range.orig_kind,
                    first_byte(range), last_byte(range), range.start_mark)
end

"""
Move `numbytes` from the range at output position `pos+1` to the output
position `pos`. If the donor range becomes empty, mark it dead with
K"TOMBSTONE" and return `true`, otherwise return `false`.

Hack alert! This is used only for managing the complicated rules related to
dedenting triple quoted strings.
"""
function steal_node_bytes!(stream::ParseStream, pos::ParseStreamPosition, numbytes)
    i = pos.output_index
    r1 = stream.ranges[i]
    r2 = stream.ranges[i+1]
    @assert span(r1) == 0
    @assert numbytes <= span(r2)
    fb2 = r2.first_byte + numbytes
    rhs_empty = fb2 > last_byte(r2)
    head2 = rhs_empty ? SyntaxHead(K"TOMBSTONE", EMPTY_FLAGS) : r2.head
    stream.ranges[i]   = TaggedRange(r1.head, r1.orig_kind,
                                     r2.first_byte, fb2 - 1,
                                     r1.start_mark)
    stream.ranges[i+1] = TaggedRange(head2, r2.orig_kind,
                                     fb2, r2.last_byte,
                                     r2.start_mark)
    return rhs_empty
end

function Base.position(stream::ParseStream)
    ParseStreamPosition(stream.next_byte, lastindex(stream.ranges))
end

"""
    emit(stream, mark, kind, flags = EMPTY_FLAGS; error=nothing)

Emit a new text span into the output which covers source bytes from `mark` to
the end of the most recent token which was `bump()`'ed. The starting `mark`
should be a previous return value of `position()`.
"""
function emit(stream::ParseStream, mark::ParseStreamPosition, kind::Kind,
              flags::RawFlags = EMPTY_FLAGS; error=nothing)
    range = TaggedRange(SyntaxHead(kind, flags), K"None", mark.input_byte,
                        stream.next_byte-1, mark.output_index+1)
    if !isnothing(error)
        _emit_diagnostic(stream, first_byte(range), last_byte(range), error=error)
    end
    push!(stream.ranges, range)
    return position(stream)
end

function _emit_diagnostic(stream::ParseStream, fbyte, lbyte; kws...)
    push!(stream.diagnostics, Diagnostic(fbyte, lbyte; kws...))
    return nothing
end

"""
Emit a diagnostic at the position of the next token

If `whitespace` is true, the diagnostic is positioned on the whitespace before
the next token. Otherwise it's positioned at the next token as returned by `peek()`.

FIXME: Rename? This doesn't emit normal tokens into the output event list!
"""
function emit_diagnostic(stream::ParseStream; whitespace=false, kws...)
    i = _lookahead_index(stream, 1, true)
    begin_tok_i = i
    end_tok_i = i
    if whitespace
        # It's the whitespace which is the error. Find the range of the current
        # whitespace.
        begin_tok_i = 1
        end_tok_i = is_whitespace(stream.lookahead[i]) ? i : max(1, i-1)
    end
    fbyte = first_byte(stream.lookahead[begin_tok_i])
    lbyte = last_byte(stream.lookahead[end_tok_i])
    _emit_diagnostic(stream, fbyte, lbyte; kws...)
    return nothing
end

function emit_diagnostic(stream::ParseStream, mark::ParseStreamPosition; kws...)
    _emit_diagnostic(stream, mark.input_byte, stream.next_byte-1; kws...)
end

function emit_diagnostic(stream::ParseStream, r::NTuple{2,ParseStreamPosition}; kws...)
    emit_diagnostic(stream, first(r), last(r); kws...)
end

function emit_diagnostic(stream::ParseStream, mark::ParseStreamPosition,
                         end_mark::ParseStreamPosition; kws...)
    _emit_diagnostic(stream, mark.input_byte, end_mark.input_byte-1; kws...)
end

#-------------------------------------------------------------------------------
# Tree construction from the list of text ranges held by ParseStream

# API for extracting results from ParseStream

"""
    build_tree(::Type{NodeType}, stream::ParseStream;
               wrap_toplevel_as_kind=nothing, kws...)

Construct a tree with `NodeType` nodes from a ParseStream using depth-first
traversal. `NodeType` must have the constructors

    NodeType(head::SyntaxHead, span::Integer)
    NodeType(head::SyntaxHead, span::Integer, children::Vector{NodeType})

A single node which covers the input is expected, but if the ParseStream has
multiple nodes at the top level, `wrap_toplevel_as_kind` may be used to wrap
them in a single node.

The tree here is constructed depth-first, but it would also be possible to use
a bottom-up tree builder interface similar to rust-analyzer. (In that case we'd
traverse the list of ranges backward rather than forward.)
"""
function build_tree(::Type{NodeType}, stream::ParseStream;
                    wrap_toplevel_as_kind=nothing, kws...) where NodeType
    stack = Vector{NamedTuple{(:range,:node),Tuple{TaggedRange,NodeType}}}()
    for (span_index, range) in enumerate(stream.ranges)
        if kind(range) == K"TOMBSTONE"
            # Ignore invisible tokens which were created but never finalized.
            # See bump_invisible()
            continue
        end

        if isempty(stack) || range.start_mark > stack[end].range.start_mark
            # A leaf node (span covering a single token):
            # [a][b][stack[end]]
            #                   [range]
            node = NodeType(head(range), span(range))
            push!(stack, (range=range, node=node))
            continue
        end
        # An interior node, span covering multiple tokens:
        #
        # [a][b][stack[end]]
        #    [        range]
        #
        # We use start_mark rather than first_byte to determine node overlap.
        # This solve the following ambiguity between invisible nodes 1 and 2:
        #
        # [a][b]|[...]
        #       |--- invisible node 1
        #       `--- invisible node 2
        #
        # Does node 2 contain node 1? Using start_mark, we can distinguish the
        # cases:
        #
        # [a][b][2][1]  [a][b][2]...
        #                     [   1]
        j = length(stack)
        while j > 1 && range.start_mark < stack[j].range.start_mark
            j -= 1
        end
        children = [stack[k].node for k = j:length(stack)]
        resize!(stack, j-1)
        node = NodeType(head(range), span(range), children)
        push!(stack, (range=range, node=node))
    end
    # show(stdout, MIME"text/plain"(), stack[1].node)
    if length(stack) == 1
        return only(stack).node
    elseif !isnothing(wrap_toplevel_as_kind)
        # Mostly for debugging
        children = [x.node for x in stack]
        return NodeType(SyntaxHead(wrap_toplevel_as_kind, EMPTY_FLAGS), children...)
    else
        error("Found multiple nodes at top level")
    end
end

"""
    sourcetext(stream::ParseStream; steal_textbuf=true)

Return the source text being parsed by this `ParseStream` as a UTF-8 encoded
string.

If `steal_textbuf==true`, this is permitted to steal the content of the
stream's text buffer. Note that this leaves the `ParseStream` in an invalid
state for further parsing.
"""
function sourcetext(stream::ParseStream; steal_textbuf=false)
    if stream.text_root isa AbstractString && codeunit(stream.text_root) == UInt8
        return stream.text_root
    elseif steal_textbuf
        return String(stream.textbuf)
    else
        # Safe default for other cases is to copy the buffer. Technically this
        # could possibly be avoided in some situations, but might have side
        # effects such as mutating stream.text_root or stealing the storage of
        # stream.textbuf
        return String(copy(stream.textbuf))
    end
end

"""
    textbuf(stream)

Return the `Vector{UInt8}` text buffer being parsed by this `ParseStream`.
"""
textbuf(stream) = stream.textbuf

first_byte(stream::ParseStream) = first_byte(first(stream.ranges))
last_byte(stream::ParseStream) = last_byte(last(stream.ranges))
any_error(stream::ParseStream) = any_error(stream.diagnostics)


text:
@testset "PLS1 Pediction Tests (in sample)" begin

	@testset "Single Column Prediction Test" begin

		X        = MLJ.table([1; 2; 3.0][:,:])
		Y        = [1; 2; 3.0]

		pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=1)  target = MLJ.UnivariateStandardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJ.fit!(pls_machine, rows=train,force=true)
        ŷ     = MLJ.predict(pls_machine, rows=train);

		@test isequal(round.(ŷ),[1; 2; 3.0])

	end

	@testset "Constant Values Prediction Tests (Ax + b) | A=0, b=1 " begin

		X        = MLJ.table([1 3;2 1;3.0 2.0])
		Y        = [1; 1; 1.0]
		try

    		pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=2)  target = MLJ.UnivariateStandardizer()

			pls_machine    = MLJ.machine(pls_pipe, X, Y)

			train = range(1,stop=length(X))
			MLJ.fit!(pls_machine, rows=train,force=true)
		catch
			@test true
		end

	end

	@testset "Linear Prediction Tests " begin

		X        = MLJ.table([1 2; 2 4; 4.0 6.0])
		Y        = [2; 4; 6.0]

    	pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=1)  target = MLJ.UnivariateStandardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJ.fit!(pls_machine, rows=train,force=true)
        ŷ = MLJ.predict(pls_machine, rows=train);

		@test isequal(round.(ŷ),[2; 4; 6.0])

		X           = MLJ.table([1 -2; 2 -4; 4.0 -6])
		Y           = [-2; -4; -6.0]

    	pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=1)  target = MLJ.UnivariateStandardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJ.fit!(pls_machine, rows=train,force=true)
        ŷ = MLJ.predict(pls_machine, rows=train);

		@test isequal(round.(ŷ),[-2; -4; -6.0])

	end

	@testset "Linear Prediction Tests (Ax + b)" begin

		X = MLJ.table([1 2; 2 4; 4.0 6;6 8; 8 10; 10.0 12.0])
		Y = [2; 4; 6.0; 8; 10; 12.0]

		train = range(1,stop=3)
		test  = range(4,stop=6)

		pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=2)  target = MLJ.UnivariateStandardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        MLJ.fit!(pls_machine, rows=train,force=true)
        pred = MLJ.predict(pls_machine, rows=test);

		@test isequal(round.(pred),[8; 10; 12.0])

		X        = MLJ.table([1 2; 2 4.0; 4.0 6; 6 8; 1 2; 2.0 4.0])
		Y        = [2; 4; 6.0; 8; 2; 4.0]

		train = range(1,stop=4)
		test  = range(5,stop=6)

    	pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=2)  target = MLJ.UnivariateStandardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        MLJ.fit!(pls_machine, rows=train,force=true)
        pred = MLJ.predict(pls_machine, rows=test);

		@test isequal(round.(pred),[2; 4])

	end

end;

@testset "PLS1 Pediction Tests (out of sample)" begin

	@testset "Linear Prediction Tests (Ax + b) | A>0" begin

		X        = MLJ.table([1 2; 2 4; 4.0 6;6 8; 8 10; 10.0 12.0])
		Y        = [2; 4; 6.0; 8.0; 10; 12]

		train = range(1,stop=3)
		test  = range(4,stop=6)

		pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=2)  target = MLJ.UnivariateStandardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        MLJ.fit!(pls_machine, rows=train,force=true)
        pred = MLJ.predict(pls_machine, rows=test);

		@test isequal(round.(pred),[8; 10; 12.0])

		X        = MLJ.table([1 2; 2 4; 4.0 6; 6 8; 8 10; 10.0 12.0])
		Y        = [4; 6; 8.0; 10; 12; 14.0]

		train = range(1,stop=3)
		test  = range(4,stop=6)

		pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=2)  target = MLJ.UnivariateStandardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        MLJ.fit!(pls_machine, rows=train,force=true)
        pred = MLJ.predict(pls_machine, rows=test);

		@test isequal(round.(pred),[10; 12; 14.0])

	end

	@testset "Linear Prediction Tests (Ax + b) | A<0" begin

		X        = MLJ.table([1 -2; 2 -4; 4.0 -6; 6 -8; 8 -10; 10.0 -12])
		Y        = [-2; -4; -6.0; -8; -10; -12.0]

		train = range(1,stop=3)
		test  = range(4,stop=6)

		pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=2)  target = MLJ.UnivariateStandardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        MLJ.fit!(pls_machine, rows=train,force=true)
        pred = MLJ.predict(pls_machine, rows=test);

		@test isequal(round.(pred),[-8; -10; -12.0])

		X        = MLJ.table([1 -2; 2 -4; 4.0 -6; 6 -8; 8 -10; 10.0 -12.0])
		Y        = [-4; -6; -8.0; -10; -12; -14.0]

		train = range(1,stop=3)
		test  = range(4,stop=6)

		pls_pipe       = MLJ.@pipeline prediction_type=:deterministic MLJ.Standardizer() PLSRegressor.PLS(n_factors=2)  target = MLJ.UnivariateStandardizer()

		pls_machine    = MLJ.machine(pls_pipe, X, Y)

        MLJ.fit!(pls_machine, rows=train,force=true)
        pred = MLJ.predict(pls_machine, rows=test);

		@test isequal(round.(pred),[-10; -12; -14.0])

	end

end;

ext: jl


text:
export GraphPlot, graphplot, graphplot!

"""
    graphplot(graph::AbstractGraph)
    graphplot!(ax, graph::AbstractGraph)

Creates a plot of the network `graph`. Consists of multiple steps:
- Layout the nodes: see `layout` attribute. The node position is accesible from outside
  the plot object `p` as an observable using `p[:node_positions]`.
- plot edges as `linesegments`-plot
- if `arrows_show` plot arrowheads as `scatter`-plot
- plot nodes as `scatter`-plot
- if `nlabels!=nothing` plot node labels as `text`-plot
- if `elabels!=nothing` plot edge labels as `text`-plot

The main attributes for the subplots are exposed as attributes for `graphplot`.
Additional attributes for the `scatter`, `linesegments` and `text` plots can be provided
as a named tuples to `node_attr`, `edge_attr`, `nlabels_attr` and `elabels_attr`.

Most of the arguments can be either given as a vector of length of the
edges/nodes or as a single value. One might run into errors when changing the
underlying graph and therefore changing the number of Edges/Nodes.

## Attributes
### Main attributes
- `layout=Spring()`: function `AbstractGraph->Vector{Point}` determines the base layout
- `node_color=scatter_theme.color`
- `node_size=scatter_theme.markersize`
- `node_marker=scatter_theme.marker`k
- `node_attr=(;)`: List of kw arguments which gets passed to the `scatter` command
- `edge_color=lineseg_theme.color`: Pass a vector with 2 colors per edge to get
  color gradients.
- `edge_width=lineseg_theme.linewidth`: Pass a vector with 2 width per edge to
  get pointy edges.
- `edge_attr=(;)`: List of kw arguments which gets passed to the `linesegments` command
- `arrow_show=Makie.automatic`: `Bool`, indicate edge directions with arrowheads?
  Defaults to `true` for `SimpleDiGraph` and `false` otherwise.
- `arrow_size=scatter_theme.markersize`: Size of arrowheads.
- `arrow_shift=0.5`: Shift arrow position from source (0) to dest (1) node.
- `arrow_attr=(;)`: List of kw arguments which gets passed to the `scatter` command

### Node labels
The position of each label is determined by the node position plus an offset in
data space.

- `nlabels=nothing`: `Vector{String}` with label for each node
- `nlabels_align=(:left, :bottom)`: Anchor of text field.
- `nlabels_distance=0.0`: Pixel distance from node in direction of align.
- `nlabels_color=labels_theme.color`
- `nlabels_offset=nothing`: `Point` or `Vector{Point}` (in data space)
- `nlabels_textsize=labels_theme.textsize`
- `nlabels_attr=(;)`: List of kw arguments which gets passed to the `text` command

### Edge labels
The base position of each label is determinded by `src + shift*(dst-src)`. The
additional `distance` parameter is given in pixels and shifts the text away from
the edge.

- `elabels=nothing`: `Vector{String}` with label for each edge
- `elabels_align=(:center, :bottom)`: Anchor of text field.
- `elabels_distance=0.0`: Pixel distance of anchor to edge.
- `elabels_shift=0.5`: Position between src and dst of edge.
- `elabels_opposite=Int[]`: List of edge indices, for which the label should be
  displayed on the opposite side
- `elabels_rotation=nothing`: Angle of text per label. If `nothing` this will be
  determined by the edge angle!
- `elabels_offset=nothing`: Additional offset in data space
- `elabels_color=labels_theme.color`
- `elabels_textsize=labels_theme.textsize`
- `elabels_attr=(;)`: List of kw arguments which gets passed to the `text` command

"""
@recipe(GraphPlot, graph) do scene
    # TODO: figure out this whole theme business
    scatter_theme = default_theme(scene, Scatter)
    lineseg_theme = default_theme(scene, LineSegments)
    labels_theme = default_theme(scene, Makie.Text)
    Attributes(
        layout = Spring(),
        # node attributes (Scatter)
        node_color = scatter_theme.color,
        node_size = scatter_theme.markersize,
        node_marker = scatter_theme.marker,
        node_attr = (;),
        # edge attributes (LineSegements)
        edge_color = lineseg_theme.color,
        edge_width = lineseg_theme.linewidth,
        edge_attr = (;),
        # arrow attributes (Scatter)
        arrow_show = automatic,
        arrow_size = scatter_theme.markersize,
        arrow_shift = 0.5,
        arrow_attr = (;),
        # node label attributes (Text)
        nlabels = nothing,
        nlabels_align = (:left, :bottom),
        nlabels_distance = 0.0,
        nlabels_color = labels_theme.color,
        nlabels_offset = nothing,
        nlabels_textsize = labels_theme.textsize,
        nlabels_attr = (;),
        # edge label attributes (Text)
        elabels = nothing,
        elabels_align = (:center, :bottom),
        elabels_distance = 0.0,
        elabels_shift = 0.5,
        elabels_opposite = Int[],
        elabels_rotation = nothing,
        elabels_offset = nothing,
        elabels_color = labels_theme.color,
        elabels_textsize = labels_theme.textsize,
        elabels_attr = (;),
    )
end

function Makie.plot!(gp::GraphPlot)
    graph = gp[:graph]

    # create initial vertex positions, will be updated on changes to graph or layout
    # make node_position-Observable available as named attribute from the outside
    gp[:node_positions] = @lift [Point(p) for p in ($(gp.layout))($graph)]

    node_pos = gp[:node_positions]

    # create two arrays for src pos and dst pos triggered by node_pos changes
    # in case of a graph change the node_position will change anyway
    edge_pos = @lift ([$node_pos[e.src] for e in edges(graph[])],
                      [$node_pos[e.dst] for e in edges(graph[])])

    # in case the edge_with is same as the number of edges
    # create a new observable which doubles the values for compat with line segments
    # https://github.com/JuliaPlots/Makie.jl/pull/992
    if length(gp.edge_width[]) == ne(graph[])
        lineseg_width = @lift repeat($(gp.edge_width), inner=2)
    else
        lineseg_width = gp.edge_width
    end

    # calculate the vectors for each edge in pixel space
    sc = Makie.parent_scene(gp)
    edge_vec_px = lift(edge_pos, sc.px_area, sc.camera.projectionview) do epos, pxa, pv
        # project should transform to 2d point in px space
        map((src, dst)->project(sc, dst) - project(sc, src), epos[1], epos[2])
    end

    # get rotation in px space (for arrow markers and edge_labels)
    edge_rotations_px = @lift map(v -> atan(v.data[2], v.data[1]), $edge_vec_px)

    # plot edges
    edge_segments = @lift vec(permutedims(hcat($edge_pos[1], $edge_pos[2])))
    edge_plot = linesegments!(gp, edge_segments;
                              color=gp.edge_color,
                              linewidth=lineseg_width,
                              gp.edge_attr...)

    # plott arrow heads
    arrow_pos = @lift $edge_pos[1] .+ $(gp.arrow_shift) .* ($edge_pos[2] .- $edge_pos[1])
    arrow_show = @lift $(gp.arrow_show) === automatic ? $graph isa SimpleDiGraph : $(gp.arrow_show)
    arrow_heads = scatter!(gp,
                           arrow_pos,
                           marker = '➤',
                           markersize = gp.arrow_size,
                           color = gp.edge_color,
                           rotations = @lift(Billboard($edge_rotations_px)),
                           strokewidth = 0.0,
                           markerspace = Pixel,
                           visible = arrow_show,
                           gp.arrow_attr...)

    # plot vertices
    vertex_plot = scatter!(gp, node_pos;
                           color=gp.node_color,
                           marker=gp.node_marker,
                           markersize=gp.node_size,
                           gp.node_attr...)

    # plot node labels
    if gp.nlabels[] !== nothing
        positions = @lift begin
            if $(gp.nlabels_offset) != nothing
                $node_pos .+ $(gp.nlabels_offset)
            else
                copy($node_pos)
            end
        end

        offset = @lift if $(gp.nlabels_align) isa Vector
            $(gp.nlabels_distance) .* align_to_dir.($(gp.nlabels_align))
        else
            $(gp.nlabels_distance) .* align_to_dir($(gp.nlabels_align))
        end

        nlabels_plot = text!(gp, gp.nlabels;
                             position=positions,
                             align=gp.nlabels_align,
                             color=gp.nlabels_color,
                             offset=offset,
                             textsize=gp.nlabels_textsize,
                             gp.nlabels_attr...)
    end

    # plot edge labels
    if gp.elabels[] !== nothing
        # rotations based on the edge_vec_px and opposite argument
        rotation = @lift begin
            if $(gp.elabels_rotation) isa Real
                # fix rotation to a single angle
                rot = $(gp.elabels_rotation)
            else
                # determine rotation by edge vector
                rot = copy($edge_rotations_px)
                for i in $(gp.elabels_opposite)
                    rot[i] += π
                end
                # if there are user provided rotation for some labels, use those
                if $(gp.elabels_rotation) isa Vector
                    for (i, α) in enumerate($(gp.elabels_rotation))
                        α !== nothing && (rot[i] = α)
                    end
                end
            end
            return rot
        end

        # positions: center point between nodes + offset + distance*normal + shift*edge direction
        positions = @lift begin
            pos = $edge_pos[1] .+ $(gp.elabels_shift) .* ($edge_pos[2] .- $edge_pos[1])

            if $(gp.elabels_offset) !== nothing
                pos .= pos .+ $(gp.elabels_offset)
            end
            pos
        end

        # calculate the offset in pixels in normal direction to the edge
        offsets = @lift begin
            offsets = map(p -> Point(-p.data[2], p.data[1])/norm(p), $edge_vec_px)
            offsets .= $(gp.elabels_distance) .* offsets
            for i in $(gp.elabels_opposite)
                offsets[i] = -1.0 * offsets[i] # flip offset if in opposite
            end
            offsets
        end

        elabels_plot = text!(gp, gp.elabels;
                             position=positions,
                             rotation=rotation,
                             offset=offsets,
                             align=gp.elabels_align,
                             color=gp.elabels_color,
                             textsize=gp.elabels_textsize,
                             gp.elabels_attr...)
    end

    return gp
end

function align_to_dir(align)
    halign, valign = align

    x = 0.0
    if halign === :left
        x = 1.0
    elseif halign === :right
        x = -1.0
    end

    y = 0.0
    if valign === :top
        y = -1.0
    elseif valign === :bottom
        y = 1.0
    end
    norm = x==y==0.0 ? 1 : sqrt(x^2 + y^2)
    return Point(x/norm, y/norm)
end


text:
export ExactSampler

"""
    ExactSampler(n_samples; seed=rand)

Constructs an exact sampler, which builds the full pdf of the quantum state and
samples it exactly.

This sampler can only be used on `indexable` spaces, and should be used only for
somewhat small systems (Closed N<14, Open N<7), as the computational cost increases
exponentially with the number of sites.

Initial seed can be set bu specifying `seed`.
"""
mutable struct ExactSampler{S} <: MonteCarloSampler
    samples_length::Int
    seed::S
end

ExactSampler(n_samples; seed=rand(UInt)) =
    ExactSampler(n_samples, seed)

## Cache
mutable struct ExactSamplerCache{H,A<:AbstractRNG,T} <: SamplerCache{ExactSampler}
    hilb::H
    loc_samples_length::Int

    rng::A
    steps_done::Int
    pdf::Vector{T}
end

function _sampler_cache(s::ExactSampler, v, hilb, net, par_cache)
    if net isa CachedNet
        if net.cache isa NNBatchedCache
            return _sampler_batched_cache(s, v, hilb, net, par_cache)
        end
    end

    loc_chain_length = Int(ceil(s.samples_length/num_workers(par_cache)))

    loc_seed         = worker_local_seed(s.seed, par_cache)
    rng              = build_rng_generator_T(prob, loc_seed)

    ExactSamplerCache(hilb,
                      loc_chain_length,
                      rng,
                      0,
                      zeros(real(eltype(v)), spacedimension(hilb)))
end
##
function init_sampler!(sampler::ExactSampler, net::Union{MatrixNet,KetNet}, σ, c::ExactSamplerCache)
    # Compute the distribution
    for i=1:spacedimension(c.hilb)
        set!(σ, c.hilb, i)
        c.pdf[i] = exp(log_prob_ψ(net, σ))
    end
    c.pdf ./= sum(c.pdf)
    cumsum!(c.pdf, c.pdf)
    c.steps_done = 0

    samplenext!(σ, sampler, net, c)
    return c
end

chain_length(s::ExactSampler, c::ExactSamplerCache) = s.samples_length

done(s::ExactSampler, σ, c) = c.steps_done >= chain_length(s, c)

function samplenext!(σ, s::ExactSampler, net::Union{MatrixNet,KetNet}, c)
    # Check termination condition, and return if verified
    done(s, σ, c) && return false

    # sample and set
    r = rand(c.rng)
    i = searchsortedfirst(c.pdf, r)
    set!(σ, i-1)
    c.steps_done += 1

    return true
end

# return findfirst(probs .> number)
function sample_from_distr(number::Real, probs::Vector)
    @inbounds for (i,p)=enumerate(probs)
        if p > number
            return i
        end
    end
end

## Batched

mutable struct ExactSamplerBatchedCache{H,A<:AbstractRNG,I,Ia,T,V,Pc} <: SamplerCache{ExactSampler}
    hilb::H
    loc_samples_length::Int

    loc_hilb_numbers::I
    all_hilb_numbers::Ia

    rng::A
    steps_done::Int
    pdf::Vector{T}

    logψ::V

    parallel_cache::Pc
end

function _sampler_batched_cache(s::ExactSampler, v, hilb, net, par_cache)
    loc_chain_length = Int(ceil(s.samples_length/num_workers(par_cache)))

    loc_seed         = worker_local_seed(s.seed, par_cache)
    rng              = build_rng_generator_T(out_similar(net), loc_seed)

    # all blocks of iterations per worker
    all_hilb_numbers = [_iterator_blocks(1:spacedimension(hilb),
                        i,
                        num_workers(par_cache)) for i=1:num_workers(par_cache)]
    loc_hilb_iter    = my_block(all_hilb_numbers, par_cache)
    all_hilb_numbers = Int32.(length.(all_hilb_numbers))

    ExactSamplerBatchedCache(hilb,
                             loc_chain_length,
                             loc_hilb_iter,
                             all_hilb_numbers,
                             rng,
                             0,
                             zeros(real(out_type(net)), spacedimension(hilb)),
                             out_similar(net),
                             par_cache)
end

##
function init_sampler!(sampler::ExactSampler,
                       net::Union{MatrixNet,KetNet},
                       σ, c::ExactSamplerBatchedCache)
    batch_sz = batch_size(net)
    # Compute the distribution
    for batch_i=Iterators.partition(c.loc_hilb_numbers, batch_sz)
        for (i, h_i)=enumerate(batch_i)
            @inbounds σᵢ = state_uview(σ, i)
            set!(σᵢ, c.hilb, h_i)
        end
        log_prob_ψ!(c.logψ, c.logψ, net, σ)

        for (i, h_i)=enumerate(batch_i)
            c.pdf[h_i] = exp.(c.logψ[i])
        end
    end
    worker_allgatherv!(c.pdf, c.all_hilb_numbers, c.parallel_cache)

    #c.pdf ./= sum(c.pdf)
    tot = sum(c.pdf)
    cumsum!(c.pdf, c.pdf)
    c.pdf ./= tot

    c.steps_done = 0

    samplenext!(σ, σ, sampler, net, c)
    return c
end

chain_length(s::ExactSampler, c::ExactSamplerBatchedCache) = c.loc_samples_length

function samplenext!(σ_out::T, σ_in::T,
                     s::ExactSampler,
                     net::Union{MatrixNet,KetNet},
                     c::ExactSamplerBatchedCache) where {T<:Union{AStateBatch, ADoubleStateBatch}}
    # Check termination condition, and return if verified
    done(s, σ, c) && return false

    for σ=states(σ_out)
        r = rand(c.rng)
        hi = searchsortedfirst(c.pdf, r)
        set!(σ, c.hilb, hi)
    end

    c.steps_done += 1
    return true
end


text:
"""
    inputfile

Abstract type for all kinds of input files
"""
abstract type inputfile end

"""
    Inputconstants = new(lx, ly, maxruntime, dumping, gravity, γ, δ, kbt)

Struct containing input parameters.

Contains `.lx` lattice points in x-direction, `.ly` lattice points in y-direction.
Other fields are `.maxruntime` for the maximal number of time steps and `.dumping` to limit the number of output files.
On top of these there are physical quantities such as `.gravity`, `.γ`, `.δ` and `.kbt`.
for the values of gravitational acceleration, fluids surface tension and the slip length.
The example relates to an quadratic lattice 20 times 20 lattice units in area.
Run for 100 lattice Boltzmann time steps only printing output every 10 time steps.
Having no gravity and a surface tension of 0.01 and a slip length of 1.

# Example
```jldoctest firsttest
julia> using JuSwalbe

julia> new_input = JuSwalbe.Inputconstants()
JuSwalbe.Inputconstants
  lx: Int64 512
  ly: Int64 512
  maxruntime: Int64 100000
  dumping: Int64 1000
  τ: Float64 1.0
  gravity: Float64 0.0
  γ: Float64 0.01
  δ: Float64 1.0
  μ: Float64 0.16666666666666666
  kbt: Float64 0.0

julia> new_input.γ
0.01
```

# References
See also: [`readinput`](@ref), [`findargument`](@ref), [`computeslip`](@ref)
"""
@with_kw struct Inputconstants <: inputfile
    lx = 512
    ly = 512
    maxruntime = 100000
    dumping = 1000
    τ = 1.0
    gravity = 0.0
    γ = 0.01
    δ = 1.0
    μ = 1 / 3 * (2 - τ) / 2 * τ
    kbt = 0.0
end

"""
    readinput(file)

Reads input parameters from a `file`.

The expected amount of parameters can be addressed with [`Inputconstants`](@ref).
For now it expects seven values for different runtime constants.

# Example
```jldoctest secondtest
julia> using JuSwalbe, DelimitedFiles

julia> args = ["Lattice_points_x" 10; "Lattice_points_y" 5; "Max_run_time" 1000; "Output_dump" 100; "Relaxation_rate" 1.0; "gravity" 0.0; "surface_tension" 0.01; "slippage" 1.0] # Generate a text file with input
8×2 Array{Any,2}:
 "Lattice_points_x"    10
 "Lattice_points_y"     5
 "Max_run_time"      1000
 "Output_dump"        100
 "Relaxation_rate"      1.0
 "gravity"              0.0
 "surface_tension"      0.01
 "slippage"             1.0

julia> writedlm("test.txt", args)

julia> test = readinput("test.txt")
JuSwalbe.Inputconstants
  lx: Int64 10
  ly: Int64 5
  maxruntime: Int64 1000
  dumping: Int64 100
  τ: Float64 1.0
  gravity: Float64 0.0
  γ: Float64 0.01
  δ: Float64 1.0
  μ: Float64 0.16666666666666666
  kbt: Float64 0.0

julia> test.lx
10

julia> test.γ
0.01

julia> test.γ + test.δ
1.01

julia> isa(test.lx + test.gravity, Int32)
false

julia> rm("test.txt")
```
"""
function readinput(file)
    # Length of the input file
    num = countlines(file)
    # Actually reading of the file and saving to an array
    input = readdlm(file)
    # Extracting of numerical values
    values = []
    arguments = ["Lattice_points_x",
                "Lattice_points_y",
                "Max_run_time",
                "Output_dump",
                "Relaxation_rate",
                "gravity",
                "surface_tension",
                "slippage",
                "viscosity",
                "k_BT"]

    for i in arguments
        push!(values, findfirst(x -> x == i, vec(input)))
    end

    lx = findargument(input, "Lattice_points_x")
    ly = findargument(input, "Lattice_points_y")
    maxruntime = findargument(input, "Max_run_time")
    dumping = findargument(input, "Output_dump")
    τ = findargument(input, "Relaxation_rate")
    gravity = findargument(input, "gravity")
    γ = findargument(input, "surface_tension")
    δ = findargument(input, "slippage")
    kbt = findargument(input, "k_BT")
    if kbt == nothing
        kbt = 0.0
    end

    runtimeconstants = Inputconstants(lx=lx, ly=ly, maxruntime=maxruntime, dumping=dumping, τ=τ, gravity=gravity, γ=γ, δ=δ, kbt=kbt)

    return runtimeconstants
end

"""
    findargument(arr, str)

Searches for a numerical value based on a str input and returns the value.

# Example
```jldoctest
julia> using JuSwalbe

julia> arr = ["hmm" 1; "yeah" 0.01; "world" 1090]
3×2 Array{Any,2}:
 "hmm"       1
 "yeah"      0.01
 "world"  1090

julia> world = findargument(arr, "world")
1090

```

# References
See also: [`readinput`](@ref)
"""
function findargument(arr::Array{Any,2}, argument::String)
    for (i,arg) in enumerate(arr[:,1])
        if arr[i,1] == argument
            value = arr[i,2]
            return value
        end
    end
end

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
using NNSparseCoding ; const N=NNSparseCoding
using Test
using Statistics

@testset "Initialization" begin
    n = 200
    b = n
    k = 80
    A,S = N.initialize_mats(n,b,k)
    @test all(A .>=  0.)
    @test all(S .>=  0.)
    @test N.nnsc_objective(A*S,A,S,0.0) ≈ 0.
    @test N.nnsc_objective(rand(n,k),A,S,1.0) > 0.
    # projection
    M = randn(n,k)
    MM = copy(M)
    N.project!(M)
    @test all(M[MM .< 0.] .== 0.0)
    mynorms = mapslices(N.norm,M;dims=1)
    @test all(isapprox.(mynorms,1.0))
end

@testset "Convergence of S alone (basis coef)" begin
    n = 80
    b = n
    k = 100
    A,S = N.initialize_mats(n,b,k)
    # let's have a sparse S
    map!(s-> rand() < 0.8 ? 0.0 : 3s ,S,S)
    X = A*S
    _,S̃ = N.initialize_mats(n,b,k)
    lambda=0.1
    nsteps=5000
    for _ in 1:nsteps
      N.updateS!(S̃,A,X,lambda)
    end
    @test all(isapprox.(S,S̃;atol=0.5))
end

@testset "Convergence of S, new interface" begin
    n = 80
    b = n
    k = 100
    lambda=0.1
    A,S = N.initialize_mats(n,b,k)
    # let's have a sparse S
    map!(s-> rand() < 0.8 ? 0.0 : 3s ,S,S)
    X = A*S
    obj = N.NNSC_optimizer(X,b,lambda)
    copy!(obj.A,A)
    N.updateAS!(obj)
    nsteps=5000
    for _ in 1:nsteps
      N.updateS!(obj)
    end
    @test all(isapprox.(S,obj.S;atol=0.5))
end

@testset "Convergence of A alone" begin
    n = 80
    b = n
    k = 100
    A,S = N.initialize_mats(n,b,k)
    # let's have a sparse S
    map!(s-> rand() < 0.85 ? 0.0 : 3s ,S,S)
    X = A*S
    Ã,_ = N.initialize_mats(n,b,k)
    lambda=0.1
    nsteps=3_000
    stepsize=0.1
    for _ in 1:nsteps
      temp = N.updateA!(Ã,S,X,lambda,stepsize)
      stepsize=temp
    end
    @test all(isapprox.(A,Ã;atol=0.3))
end
#
# @testset "Convergence of A, other interface" begin
#     n = 80
#     b = n
#     k = 100
#     lambda=0.1
#     A,S = N.initialize_mats(n,b,k)
#     # let's have a sparse S
#     map!(s-> rand() < 0.85 ? 0.0 : 3s ,S,S)
#     X = A*S
#     obj = N.NNSC_optimizer(X,b,lambda)
#     copy!(obj.S,S)
#     N.updateAS!(obj)
#     nsteps=5_000
#     stepsize=0.02
#     objval = N.nnsc_objective(obj)
#     for _ in 1:nsteps
#       stepsize,objval,_= N.updateA!(obj,stepsize,objval)
#     end
#     @test all(isapprox.(A,obj.A;atol=0.3))
# end
#
# @testset "Optimization of A and S" begin
#     n = 80
#     b = n
#     k = 300
#     A,S = N.initialize_mats(n,b,k)
#     # let's have a sparse S
#     map!(s-> rand() < 0.95 ? 0.0 : 3s ,S,S)
#     X = A*S
#     lambda=0.1
#     o = N.NNSC_optimizer(X,b,lambda)
#     N.nnsc_objective(o)
#     nsteps=10_000
#     stepsize=0.02
#     obj_temp = N.nnsc_objective(o)
#     for _ = 1:nsteps
#       (stepsize,obj_temp,_) = N.updateA!(o, stepsize,obj_temp)
#       N.updateS!(o)
#     end
#     @test all(isapprox.(X,o.AS;atol=1.0))
# end
#
@testset "Transformations" begin
    n = 80
    b = n
    k = 150
    A,S = N.initialize_mats(n,b,k)
    map!(s-> rand() < 0.85 ? 0.0 : 3s ,S,S)
    X = A*S
    broadcast!(-,X,X,mean(X;dims=2))
    whit,Mp,Mm = N.preprocess_plus_minus(X,N.PreprocessWhiten)
    @test (!any(Mp .< 0)) && (!any(Mm .< 0))
    X̃ = N.revert_plus_minus(Mp,Mm,whit)
    @test all(isapprox.(X,X̃))
    pcadims = 40
    whit, Mp,Mm = N.preprocess_plus_minus(X,N.PreprocessPCA(pcadims))
    @test (!any(Mp .< 0)) && (!any(Mm .< 0))
    X̃ = N.revert_plus_minus(Mp,Mm,whit)
    @test all(isapprox.(X,X̃ ; atol=2))
    whit,Mp,Mm = N.preprocess_plus_minus(X,N.PreprocessNot)
    @test (!any(Mp .< 0)) && (!any(Mm .< 0))
    X̃ = N.revert_plus_minus(Mp,Mm,whit)
    @test all(isapprox.(X,X̃))
end


text:
immutable HessianFunction{T <: Function}
    f::T
    name::String
end

hess_funcs = (
    HessianFunction(
        FiniteDiff.hessian,
        "FiniteDiff.hessian(f, x)",
    ),
    HessianFunction(
        Calculus.finite_difference_hessian,
        "Calculus.finite_difference_hessian(f, x)",
    ),
)


text:
function test_data()
	x = simulate_random_snparray(1000, 1000, undef)
	z = ones(1000, 1)
	y = rand(1000)
    v = IHTVariables(x, z, y, 1, 10, Int[], Float64[]) #J = 1, k = 10

	return (x, z, y, v)
end

function test_correlated_data()
    x = simulate_correlated_snparray(1000, 1000, undef)
    z = ones(1000, 1)
    y = rand(1000)
    v = IHTVariables(x, z, y, 1, 10, Int[], Float64[]) #J = 1, k = 10

    return (x, z, y, v)
end

@testset "loglikelihood" begin
	Random.seed!(2019)

    d = Normal()
    μ = rand(1000000)
    y = [rand(Normal(μi)) for μi in μ]	#simulate deviance ϕ ≈ 1

    @test isapprox(MendelIHT.loglikelihood(d, y, μ), sum(logpdf.(Normal.(μ), y)), atol=1e-1)

    d = Bernoulli()
    p = rand(10000)
    y = rand(0.0:1.0, 10000)

    @test isapprox(MendelIHT.loglikelihood(d, y, p), sum(logpdf.(Bernoulli.(p), y)), atol=1e-8)

    d = Poisson()
    λ = rand(10000)
    y = Float64.(rand(1:100, 10000))

    @test isapprox(MendelIHT.loglikelihood(d, y, λ), sum(logpdf.(Poisson.(λ), y)), atol=1e-8)

    d = NegativeBinomial()
    p = rand(1000000)
    r = ones(1000000)
    y = Float64.([rand(NegativeBinomial(1.0, p_i)) for p_i in p])

    @test isapprox(MendelIHT.loglikelihood(d, y, p), sum(logpdf.(NegativeBinomial.(r, r ./ (p .+ r)), y)), atol=1e-6)
end

@testset "deviance" begin
	# Need more test other than normal distribution!!
	d = Normal
	y = randn(1000)
	μ = randn(1000)
    @test sum(abs2, y .- μ) ≈ MendelIHT.deviance(d(), y, μ)
end

@testset "update_μ!" begin
	Random.seed!(2019)

    d = Normal
    μ = zeros(1000)
    xb = rand(1000)
    MendelIHT.update_μ!(μ, xb, canonicallink(d()))

    @test all(isapprox(μ, xb, atol=1e-8))

    d = Bernoulli
    p = zeros(1000)
    xb = rand(1000)
    MendelIHT.update_μ!(p, xb, canonicallink(d()))

    @test all(isapprox(p, 1 ./ (1 .+ exp.(-1.0 .* xb)), atol=1e-8))

    d = Poisson
    λ = zeros(1000)
    xb = rand(1000)
    MendelIHT.update_μ!(λ, xb, canonicallink(d()))

    @test all(isapprox(λ, (exp.(xb)), atol=1e-8))

    d = NegativeBinomial
    p = zeros(1000)
    xb = rand(1000)
    MendelIHT.update_μ!(p, xb, LogLink()) #use loglink for NegativeBinomial

    @test all(isapprox(p, exp.(xb), atol=1e-8))
end

@testset "update_xb! and check_covariate_supp!" begin
   	Random.seed!(1111)
	x, z, y, v = test_data()
    v.idx[1:10] .= trues(10)
    v.b .= rand(1000)
	tmp_x = convert(Matrix{Float64}, @view(x[:, v.idx]), center=true, scale=true)

    @test size(v.xk, 2) == 9 #IHTVariable is initialized to have k - 1 columns

	MendelIHT.check_covariate_supp!(v)

	@test size(v.xk, 2) == 10

    MendelIHT.update_xb!(v, x, z)

    @test all(v.xb .≈ tmp_x * v.b[1:10])
    @test all(v.zc .== 0.0)
    @test all(abs.(v.b) .<= 30)
end

@testset "score!" begin
    # Not sure how to "really" test this function, so I'm throwing in a lot of input/outputs
    # as they are calculated right now, because the code seems to work well

   	Random.seed!(1993)
    d = Normal()
    l = IdentityLink()
   	(x, z, y, v) = test_data()
	xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
	score!(d, l, v, xbm, z, y)

	@test all(v.df[1:3] .≈ [8.057330896198419; -2.229495636684423; 9.948528223937274])
	@test v.df2[1] ≈ 500.8665816573597

   	(x, z, y, v) = test_data()
	xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
	score!(d, l, v, xbm, z, y)

	@test all(v.df[1:3] .≈ [-12.569757729532215; 8.237144382138629; 9.625154193038197])
	@test v.df2[1] ≈ 503.8433996263735

   	(x, z, y, v) = test_data()
	xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
	score!(d, l, v, xbm, z, y)

	@test all(v.df[1:3] .≈ [-2.2634046444623226; -5.266454260596382; -1.5449038662162529])
	@test v.df2[1] ≈ 507.55935638764254
end

@testset "_iht_gradstep" begin
   	Random.seed!(1111)

	x, z, y, v = test_data()
	xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);

	v.b[1:3] .= [1; 2; 3]
	v.df[1:3] .= [1; 2; 3]
	b = copy(v.b)
	df = copy(v.df)
	J = 1
	k = 2
	η = 0.9

	MendelIHT._iht_gradstep(v, η, J, k, [v.df; v.df2]) # this should keep 1 * 2 = 2 elements

	@test v.b[1] ≈ 0.0 # because first entry is smallest, it should be set to 0
	@test v.b[2] ≈ (b + η*df)[2]
	@test v.b[3] ≈ (b + η*df)[3]
	@test all(v.b[4:end] .≈ 0.0)
	@test all(v.c .≈ 0.0)
	@test all(v.df .≈ df)
end

@testset "_choose!" begin
	J = 1
	k = 10

   	Random.seed!(1111)
	x, z, y, v = test_data()
	v.idx[1:10] .= trues(10)
	MendelIHT._choose!(v, J, k)

	@test all(v.idx[1:10] .== trues(10))
	@test v.idc[1] == false

	v.idx[1:11] .= trues(11)
	MendelIHT._choose!(v, J, k)

	@test sum(v.idx[1:11]) == 10

	x, z, y, v = test_data()
	v.idc .= true
	MendelIHT._choose!(v, J, 1)

	@test v.idc[1] == true
	@test sum(v.idc) == 1
	@test sum(v.idx) == 0

	x, z, y, v = test_data()
	v.idc .= true
	v.idx[1:10] .= trues(10)
	MendelIHT._choose!(v, J, k)

	@test sum(v.idx) + sum(v.idc) == 10
end

@testset "_iht_backtrack" begin
	@test MendelIHT._iht_backtrack_(-150.0, -100.0, 1, 10) == true
	@test MendelIHT._iht_backtrack_(-50.0, -100.0, 1, 10) == false
	@test MendelIHT._iht_backtrack_(150.0, 100.0, 1, 10) == false
	@test MendelIHT._iht_backtrack_(50.0, 100.0, 1, 10) == true
	@test MendelIHT._iht_backtrack_(50.0, 100.0, 0, 3) == true
	@test MendelIHT._iht_backtrack_(50.0, 100.0, 11, 10) == false
	@test MendelIHT._iht_backtrack_(50.0, 100.0, 11, 11) == false
end

@testset "std_reciprocal" begin
	Random.seed!(2019)

	# first compute the correct answer by converting each column to floats and call std() directly
	n, p = 10000, 10000
	x = simulate_random_snparray(n, p, undef)
	storage = zeros(n)
	answer  = zeros(p)
    for i in 1:p
    	copyto!(storage, @view(x[:, i]))
        answer[i] = std(storage)
    end
    answer .= 1.0 ./ answer

    # now compute the std of ecah snp
	xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
    μ = 2.0maf(x)
    std_vector = std_reciprocal(xbm, μ)

    @test all(isapprox.(std_vector, answer, atol=1e-3))
end

@testset "standardize!" begin
	Random.seed!(1017)
    z = rand(1000, 1000)
    standardize!(z)
    @test isapprox(0.0, mean(z[:, 1]), atol=1e-8)
    @test isapprox(0.0, mean(z[:, 10]), atol=1e-8)
    @test isapprox(0.0, mean(z[:, 100]), atol=1e-8)
    @test isapprox(1.0, std(z[:, 1]), atol=1e-8)
    @test isapprox(1.0, std(z[:, 10]), atol=1e-8)
    @test isapprox(1.0, std(z[:, 100]), atol=1e-8)

    z = rand(1000, 1000)
    col1_μ = mean(z[:, 1])
    col1_σ = std(z[:, 1])
    standardize!(@view(z[:, 2:end]))
    @test isapprox(col1_μ, mean(z[:, 1]), atol=1e-8)
    @test isapprox(0.0, mean(z[:, 10]), atol=1e-8)
    @test isapprox(0.0, mean(z[:, 100]), atol=1e-8)
    @test isapprox(col1_σ, std(z[:, 1]), atol=1e-8)
    @test isapprox(1.0, std(z[:, 10]), atol=1e-8)
    @test isapprox(1.0, std(z[:, 100]), atol=1e-8)
end

@testset "project_k!" begin
    x = rand(100000)
    k = 100
    p = sortperm(x, rev = true)
    top_k_index = p[1:k]
	last_k_index = p[k+1:end]
	project_k!(x, k)

	@test all(x[top_k_index] .!= 0.0)
	@test all(x[last_k_index] .== 0.0)
end

@testset "project_group_sparse!" begin
	Random.seed!(1914)

    m, n, k = 2, 3, 10 #2 active groups, 3 active predictors per group, 10 total predictors
	y = randn(k);
	group = rand(1:5, k);
	x = copy(y)
	project_group_sparse!(x, group, m, n)

	non_zero_position = findall(!iszero, x)
	non_zero_entries = x[non_zero_position]
	@test all(non_zero_entries .== y[non_zero_position])
	@test all(non_zero_position .== [2; 6; 8; 10])
	@test all(non_zero_entries .≈ [0.44267476307372516; -1.4199877945915547;
 								  -0.8857806660404711; -0.06177816577797742;])

	# test project_group_sparse! is equivalent to project_k! when max group = 1
	m, n, k = 1, 10, 100000
	group = ones(Int, k) #everybody is in the same group
	y = randn(k)
	x = copy(y)
	project_k!(x, n)
	project_group_sparse!(y, group, m, n)
	@test all(x .== y)
end

@testset "maf_weights" begin
    Random.seed!(33)
    x = simulate_random_snparray(1000, 10000, undef)
    m = maf(x)
    p = maf_weights(x)

	@test eltype(p) <: AbstractFloat
	@test all(p .>= 1.0)
	@test p[1] ≈ 1 / (2.0sqrt(m[1] * (1 - m[1])))
	@test p[2] ≈ 1 / (2.0sqrt(m[2] * (1 - m[2])))

	p = maf_weights(x, max_weight=2.0)
	@test all(1.0 .<= p .<= 2)
	@test p[1] ≈ 1 / (2.0sqrt(m[1] * (1 - m[1])))
	@test p[15] == 2.0
end

@testset "save_prev!" begin
   	Random.seed!(1111)
	x, z, y, v = test_data()

	v.b .= rand(1000)
	v.xb .= rand(1000)
	v.idx .= bitrand(1000)
	v.idc .= true
	v.c .= rand()
	v.zc .= rand(1000)

	save_prev!(v)

	@test all(v.b .== v.b0)
	@test all(v.xb .== v.xb0)
	@test all(v.idx .== v.idx0)
	@test all(v.idc .== v.idc0)
	@test all(v.c .== v.c0)
	@test all(v.zc .== v.zc0)
end

@testset "iht_stepsize" begin
    # Not sure how to "really" test this function, so I'm throwing in a lot of input/outputs
    # as they are calculated right now, because the code seems to work well

    Random.seed!(1234)
    d = Normal()
    l = canonicallink(d)
    x, z, y, v = test_data()
    v.df .= rand(1000)
    v.xk .= rand(1000, 9)
    v.μ .= rand(1000)
    v.idx[1:9] .= trues(9)

    @test MendelIHT.iht_stepsize(v, z, d, l) ≈ 0.0006414336637728858

    d = Poisson()
    l = canonicallink(d)
    x, z, y, v = test_data()
    v.df .= rand(1000)
    v.xk .= rand(1000, 9)
    v.μ .= rand(1000)
    v.idx[1:9] .= trues(9)

    @test MendelIHT.iht_stepsize(v, z, d, l) ≈ 6.440657641991925e-5

    d = NegativeBinomial()
    l = LogLink()
    x, z, y, v = test_data()
    v.df .= rand(1000)
    v.xk .= rand(1000, 9)
    v.μ .= rand(1000)
    v.idx[1:9] .= trues(9)

    @test MendelIHT.iht_stepsize(v, z, d, l) ≈ 6.606936547001659e-5

    d = Bernoulli()
    l = canonicallink(d)
    x, z, y, v = test_data()
    v.df .= rand(1000)
    v.xk .= rand(1000, 9)
    v.μ .= rand(1000)
    v.idx[1:9] .= trues(9)

    @test MendelIHT.iht_stepsize(v, z, d, l) ≈ 0.0004864093320174198
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Main function to print the numbers.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export pretty_number

"""
    pretty_number([io::IO, | String, ]number::Number; kwargs...)

Print the `number`.

If the first argument is an `io`, then the number is printed to it. If it is
a `String`, then a string is returned with the printed number. It it is omitted,
then it defaults to `stdout`.

# Keywords

The following keywords are available to modify the number printing. Those
arguments depends on the type of the number.

## Rational

If `number` is `Rational`, then the following keywords are available:

- `compact::Bool`: If `true`, then the rational number will be printed
    compactly, in one line like `³/₄`. Otherwise, the rational number is printed
    using multiple lines, like (**Default** = `true`):

        123
        ————
        4567

## Numbers

Otherwise, the `number` is printed using the scientific notation in the base 10.
In this case, the following keywords are available:

- `always_print_base::Bool`: If `true`, then the base is always printed even if
    the base exponent is 0. (**Default** = `false`)
- `multiplication_sign::Char`: The multiplication sign that will be used between
    the significand and the decimal base, common options are `'⋅'` and `'×'`.
    (**Default** = `'×'`)
- `significand_format::String`: The format that will be used to print the
    signifcand, as described by the function [`Printf.@printf`](@ref).
    (**Default** = `"%g"`)
- `show_base::Bool`: If `true`, then the base will be printed. Otherwise, it
    will be omitted. (**Default** = `true`)
- `show_significand::Bool`: If `true`, then the significand will be printed.
    Otherwise, it will be omitted. (**Default** = `true`)
- `new_decimal_base::Union{Nothing, Number}`: If it is a number, then the
    decimal base of the number will be converted to this number. If it is
    `nothing`, then the base is not changed. (**Default** = `nothing`)

# Examples

```julia
julia> pretty_number(19//86)
¹⁹/₈₆

julia> pretty_number(19//86; compact = false)
19
——
86

julia> pretty_number(1906.1896)
1.90619 · 10³

julia> pretty_number(1906.1896, significand_format = "%.10f")
1.9061896000 · 10³

julia> pretty_number(1906.1896; new_decimal_base = 4)
0.190619 · 10⁴
```
"""
function pretty_number(number::Number; kwargs...)
    return pretty_number(stdout, number; kwargs...)
end

function pretty_number(io::IO, number::Number; kwargs...)
    _pn_text(io, number; kwargs...)
    return nothing
end

function pretty_number(::Type{String}, number::Number; kwargs...)
    io = IOBuffer()
    pretty_number(io, number; kwargs...)
    str = String(take!(io))
    return str
end


text:
@testset "test_discrete" begin

C_111 = ss([-5], [2], [3], [0])
C_212 = ss([-5 -3; 2 -9], [1; 2], [1 0; 0 1], [0; 0])
C_221 = ss([-5 -3; 2 -9], [1 0; 0 2], [1 0], [0 0])
C_222_d = ss([-5 -3; 2 -9], [1 0; 0 2], [1 0; 0 1], [1 0; 0 1])

@test c2d(ss(4*[1 0; 0 1]), 0.5, :zoh) == (ss(4*[1 0; 0 1], 0.5), zeros(0, 2))
@test c2d(ss(4*[1 0; 0 1]), 0.5, :foh) == (ss(4*[1 0; 0 1], 0.5), zeros(0, 2))
@test_c2d(c2d(C_111, 0.01, :zoh),
    ss([0.951229424500714], [0.019508230199714396], [3], [0], 0.01), [1 0], true)
@test_c2d(c2d(C_111, 0.01, :foh),
    ss([0.951229424500714], [0.01902855227625244], [3], [0.029506188017136226], 0.01),
    [1 -0.009835396005712075], true)
@test_c2d(c2d(C_212, 0.01, :zoh),
    ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
    [0.009466805409666932; 0.019219966830212765], [1 0; 0 1], [0; 0], 0.01),
    [1 0 0; 0 1 0], true)
@test_c2d(c2d(C_212, 0.01, :foh),
    ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482],
    [0.008957940478201584; 0.018468989584974498], [1 0; 0 1], [0.004820885889482196;
    0.009738343195298675], 0.01), [1 0 -0.004820885889482196; 0 1 -0.009738343195298675], true)
@test_c2d(c2d(C_221, 0.01, :zoh),
    ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
    [0.009753161420545834 -0.0002863560108789034; 9.54520036263011e-5 0.019124514826586465],
    [1.0 0.0], [0.0 0.0], 0.01), [1 0 0 0; 0 1 0 0], true)
@test_c2d(c2d(C_221, 0.01, :foh),
    ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482],
    [0.009511049106772921 -0.0005531086285713394; 0.00018436954285711309 0.018284620042117387],
    [1 0], [0.004917457305816479 -9.657141633428213e-5], 0.01),
    [1 0 -0.004917457305816479 9.657141633428213e-5;
    0 1 -3.219047211142736e-5 -0.009706152723187249], true)
@test_c2d(c2d(C_222_d, 0.01, :zoh),
    ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
    [0.009753161420545834 -0.0002863560108789034; 9.54520036263011e-5 0.019124514826586465],
    [1 0; 0 1], [1 0; 0 1], 0.01), [1 0 0 0; 0 1 0 0], true)
@test_c2d(c2d(C_222_d, 0.01, :foh),
    ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482],
    [0.009511049106772921 -0.0005531086285713394; 0.00018436954285711309 0.018284620042117387],
    [1 0; 0 1], [1.0049174573058164 -9.657141633428213e-5; 3.219047211142736e-5 1.0097061527231872], 0.01),
    [1 0 -0.004917457305816479 9.657141633428213e-5; 0 1 -3.219047211142736e-5 -0.009706152723187249], true)

# Test some false
@test_c2d(c2d(2C_111, 0.01, :zoh),  # Factor 2
    ss([0.951229424500714], [0.019508230199714396], [3], [0], 0.01), [1 0], false)
@test_c2d(c2d(C_111, 0.01, :foh),   #Wrong C
    ss([0.951229424500714], [0.01902855227625244], [3*3], [0.029506188017136226], 0.01),
    [1 -0.009835396005712075], false)
@test_c2d(c2d(C_212, 0.1, :zoh),    # Wrong Ts
    ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
    [0.009466805409666932; 0.019219966830212765], [1 0; 0 1], [0; 0], 0.01),
    [1 0 0; 0 1 0], false)

# Test c2d for transfer functions
G = tf([1, 1], [1, 3, 1])
Gd = c2d(G, 0.2)
@test Gd ≈ tf([0, 0.165883310712090, -0.135903621603238], [1.0, -1.518831946985175, 0.548811636094027], 0.2) rtol=1e-14

# Issue #391
@test c2d(tf(C_111), 0.01, :zoh) ≈ tf(c2d(C_111, 0.01, :zoh)[1])
@test c2d(tf(C_111), 0.01, :foh) ≈ tf(c2d(C_111, 0.01, :foh)[1])

# c2d on a zpk model should arguably return a zpk model
@test_broken typeof(c2d(zpk(G), 1)) <: TransferFunction{<:ControlSystems.SisoZpk}

# ERRORS
@test_throws ErrorException c2d(ss([1], [2], [3], [4], 0.01), 0.01)   # Already discrete
@test_throws ErrorException c2d(ss([1], [2], [3], [4], -1), 0.01)     # Already discrete

# d2c
@static if VERSION > v"1.4" # log(matrix) is buggy on previous versions, should be fixed in 1.4 and back-ported to 1.0.6
    @test d2c(c2d(C_111, 0.01)[1]) ≈ C_111
    @test d2c(c2d(C_212, 0.01)[1]) ≈ C_212
    @test d2c(c2d(C_221, 0.01)[1]) ≈ C_221
    @test d2c(c2d(C_222_d, 0.01)[1]) ≈ C_222_d
    @test d2c(Gd) ≈ G

    sys = ss([0 1; 0 0], [0;1], [1 0], 0)
    sysd = c2d(sys, 1)[1]
    @test d2c(sysd) ≈ sys
end

end


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using ContactImplicitMPC
using LinearAlgebra

# ## Simulation
s = get_simulation("pushbot", "flat_2D_lc", "flat");
model = s.model
env = s.env

# ## Reference Trajectory
h = 0.04
H = 100
ref_traj = contact_trajectory(model, env, H, h)
ref_traj.h
qref = [0.0; 0.0]
ur = zeros(model.nu)
γr = zeros(model.nc)
br = zeros(model.nc * friction_dim(env))
ψr = zeros(model.nc)
ηr = zeros(model.nc * friction_dim(env))
wr = zeros(model.nw)

# ## Set Reference
for t = 1:H
	ref_traj.z[t] = pack_z(model, env, qref, γr, br, ψr, ηr)
	ref_traj.θ[t] = pack_θ(model, qref, qref, ur, wr, model.μ_world, ref_traj.h)
end

# ## Initial conditions
# q0 = ContactImplicitMPC.SVector{2}([0.0 * π, 0.0])
q1 = [0.0 * π, 0.0]
v1 = [0.0; 0.0]

# ## Simulator
sim = simulator(s, H, h=h)

# ## Simulate
status = simulate!(sim, q1, v1)

# ## MPC setup
N_sample = 2
H_mpc = 40
h_sim = h / N_sample
H_sim = 1000
κ_mpc = 1.0e-4

# ## Slow Recovery
obj = TrackingVelocityObjective(model, env, H_mpc,
	q = [Diagonal([12*(t/H_mpc)^2; 2.0*(t/H_mpc)^4]) for t = 1:H_mpc-0],
	v = [Diagonal([1; 0.01] ./ (h^2.0)) for t = 1:H_mpc-0],
	u = [Diagonal([100; 1]) for t = 1:H_mpc-0],
	γ = [Diagonal(1.0e-100 * ones(model.nc)) for t = 1:H_mpc-0],
	b = [Diagonal(1.0e-100 * ones(model.nc * friction_dim(env))) for t = 1:H_mpc]);

# ## Fast Recovery
obj = TrackingVelocityObjective(model, env, H_mpc,
    q = [Diagonal([12*(t/H_mpc)^2; 12*(t/H_mpc)^2]) for t = 1:H_mpc-0],
	v = [Diagonal([1; 0.01] ./ (h^2.0)) for t = 1:H_mpc-0],
    u = [Diagonal([100; 1]) for t = 1:H_mpc-0],
    γ = [Diagonal(1.0e-100 * ones(model.nc)) for t = 1:H_mpc-0],
    b = [Diagonal(1.0e-100 * ones(model.nc * friction_dim(env))) for t = 1:H_mpc]);

# ## Policy
p = ci_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
    n_opts = NewtonOptions(
		r_tol = 3e-4,
		max_iter = 10,
		max_time = ref_traj.h/2, # HARD REAL TIME
		),
    mpc_opts = CIMPCOptions());

# ## Disturbances
idx_d1 = 20
idx_d2 = idx_d1 + 200
idx_d3 = idx_d2 + 80
idx_d4 = idx_d3 + 200
idx_d5 = idx_d4 + 30
idx = [idx_d1, idx_d2, idx_d3, idx_d4, idx_d5]
impulses = [[-5.5; 0.0], [+5.5; 0.0], [+5.5; 0.0], [-1.5; 0.0], [-6.5; 0.0]]
d = impulse_disturbances(impulses, idx);

# ## Initial Conditions
q1_sim = [0.0, 0.0]
v1_sim = [0.0; 0.0]

# ## Simulator
sim = simulator(s, H_sim, h=h_sim, policy=p, dist=d)

# ## Simulate
status = simulate!(sim, verbose = true)

# ## Visualizer
vis = ContactImplicitMPC.Visualizer()
ContactImplicitMPC.render(vis)

# ## Visualize
anim = visualize_robot!(vis, model, sim.traj, sample = 1)
pθ_right = generate_pusher_traj(d, sim.traj, side=:right)
pθ_left  = generate_pusher_traj(d, sim.traj, side=:left)
visualize_disturbance!(vis, model, pθ_right, anim=anim, sample=1, offset=0.05, name=:PusherRight);
visualize_disturbance!(vis, model, pθ_left,  anim=anim, sample=1, offset=0.05, name=:PusherLeft);

# ## Timing result
# Julia is [JIT-ed](https://en.wikipedia.org/wiki/Just-in-time_compilation) so re-run the MPC setup through Simulate for correct timing results.
process!(sim.stats, N_sample) # Time budget
H_sim * h_sim / sum(sim.stats.policy_time) # Speed ratio


text:
using transD_GP, Distributed
## make options for the multichannel lengthscale GP
fileprefix = "OrstUhn_"
Klog10λ = transD_GP.GP.OrstUhn()
nminlog10λ, nmaxlog10λ = 2, 30
pnorm = 2.
log10bounds = [-0.2 0.5]
λlog10λ = [0.05abs(diff([extrema(znall)...])[1])]
δlog10λ = 0.06
demean = false
sdev_poslog10λ = [0.07abs(diff([extrema(znall)...])[1])]
sdev_proplog10λ = 0.2*diff(log10bounds, dims=2)[:]
xall = permutedims(collect(znall))
xbounds = permutedims([extrema(znall)...])
## Initialize a lengthscale model using these options
Random.seed!(12)
optlog10λ = transD_GP.OptionsStat(nmin = nminlog10λ,
                        nmax = nmaxlog10λ,
                        xbounds = xbounds,
                        fbounds = log10bounds,
                        fdataname = fileprefix,
                        xall = xall,
                        λ = λlog10λ,
                        δ = δlog10λ,
                        demean = demean,
                        sdev_prop = sdev_proplog10λ,
                        sdev_pos = sdev_poslog10λ,
                        pnorm = pnorm,
                        quasimultid = false,
                        needλ²fromlog = true,
                        updatenonstat = true,
                        K = Klog10λ,
                        timesλ = 3.7,
                        peskycholesky = true
                        )
## make options for the nonstationary actual properties GP
nmin, nmax = 2, 20
fbounds = [-0.5 2.5]
δ = 0.1
sdev_prop = 0.04*diff(fbounds, dims=2)[:]
sdev_pos = [0.02abs(diff([extrema(znall)...])[1])]
demean_ns = true
K = transD_GP.GP.OrstUhn()
## Initialize model for the nonstationary properties GP
Random.seed!(13)
opt = transD_GP.OptionsNonstat(optlog10λ,
                        nmin = nmin,
                        nmax = nmax,
                        fbounds = fbounds,
                        δ = δ,
                        demean = demean_ns,
                        sdev_prop = sdev_prop,
                        sdev_pos = sdev_pos,
                        pnorm = pnorm,
                        K = K
                        )

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
using Plots
using LaTeXStrings
using FrictionModels

# Develop a function describing the relative velocity between the two bodies.
velocity(t) = 1.5 * sin(10.0 * pi * t)

# Develop a function describing the normal force between the two bodies.
normal(t) = 100.0

# A function describing the relative position between two bodies.  Notice, this
# is not needed by the Lu-Gre model, and as such, any function can be defined
# that only takes a single argument.
position(t) = 0.0

# Define the friction model
mdl = ElastoPlasticModel(
    0.25,
    0.15,
    0.25,
    1e6,
    5e2,
    1e-6,
    7e-6,
    0.05
)

# Compute the friction model over a segment of time
tspan = [0.0, 1.0]
rsp = friction(mdl, tspan, normal, position, velocity, [0.0])

# Plot the results
plot(
    velocity.(rsp.t),
    rsp.f ./ (mdl.static_coefficient .* normal.(rsp.t)),
    xlabel = L"v(t)",
    ylabel = L"\frac{F}{\mu_s N}",
    label = false
)


text:
function print_modifier(io, mod, substring)
    if mod == "^"
        to_superscript(io, substring)
    elseif mod == "_"
        to_subscript(io, substring)
    elseif mod == "bb"
        to_blackboardbold(io, substring)
    elseif mod == "bf"
        to_boldface(io, substring)
    elseif mod == "it"
        to_italic(io, substring)
    elseif mod == "cal"
        to_caligraphic(io, substring)
    elseif mod == "frak"
        to_frakture(io, substring)
    elseif mod == "mono"
        to_mono(io, substring) # leave unmodified for now
    else
        error("Modifier $mod not supported")
    end
end
"""
Base findnext doesn't handle utf8 strings correctly
"""
function utf8_findnext(A::AbstractString, v::Char, idx::Integer)
    while true
        lastidx = idx
        elem_idx = iterate(A, idx)
        elem_idx === nothing && break
        elem, idx = elem_idx
        elem == v && return lastidx
    end
    0
end

function to_latex(text)
    io = IOBuffer()
    charidx = iterate(text)
    charidx === nothing && return ""
    char, idx = charidx
    started = true
    while true
        started || (charidx = iterate(text, idx))
        started = false
        charidx === nothing && break
        char, idx = charidx
        if char in ('^', '_', '\\')
            mod = string(char)
            if mod == "\\"
                ss = SubString(text, idx, lastindex(text))
                for mod_candidate in ("bb", "bf", "it", "cal", "frak", "mono")
                    if startswith(ss, mod_candidate)
                        mod = mod_candidate
                        break
                    end
                end
                if mod == "\\" # no match was found
                    # is this a latex symbol?
                    for (k, v) in latex_symbol_map
                        if startswith(ss, k)
                            print(io, v) # replace
                            for i in 1:length(k) # move forward
                                idx = nextind(text, idx)
                            end
                            break
                        end
                    end
                    continue # ignore '\' mod
                else
                    for i in 1:length(mod) # move forward
                        idx = nextind(text, idx)
                    end
                end
            end
            char, idx = iterate(text, idx)
            if char == '{'
                i = utf8_findnext(text, '}', idx)
                if i == 0
                    error("Invalid latex. Couldn't find matching } in $(text[idx:end])")
                end
                print_modifier(io, mod, SubString(text, idx, prevind(text, i)))
                char, idx = iterate(text, i)
            else
                print_modifier(io, mod, char)
            end
        else
            print(io, char)
        end
    end
    return String(take!(io))
end


text:
include("ow_tools.jl")
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    OW2LevelEstimator(data::Array{Float64,2}, hdfm::HDFMStruct)

Description:
Estimate a two-level HDFM using the Otrok-Whiteman approach.
Both the latent factors and hyperparameters are estimated using the Bayesian approach outlined in Otrok and Whiteman (1998).

Inputs:
- data = Matrix with each column being a data series.
- hdfm = Model structure specification.

Outputs:
- results = HDMF Bayesian estimator-generated MCMC posterior distribution samples and their means for latent factors and hyperparameters.
"""
function OW2LevelEstimator(data::Array{Float64,2}, hdfm::HDFMStruct)

    # Unpack two-level HDFM parameters
    @unpack nlevels, nfactors, factorassign, factorlags, errorlags, ndraws, burnin = hdfm

    fassign = factorassign
    flags = factorlags
    varlags = errorlags

    # Save data & its size
    y = data                    # save data in new matrix
    capt, nvars = size(y)       # nvar = # of variables; capt = # of time periods in complete sample
    nvar = nvars

    # Store factor and parameter counts
    nfact = sum(nfactors)       # # of factors
    arlag = max(flags...)       # autoregressive lags in the dynamic factors
    arterms = max(varlags...)   # number of AR lags to include in each observable equation
    nreg = 1 + nlevels          # # of regressors in each obs. eq. (intercept + factors)

    # Count number of variables each 2nd-level factor loads on
    # using a vector called `fnvars`
    # where the i-th entry represents the i-th 2nd-level factors
    fnvars = zeros(Int, nfact - 1)
    for i in 1:(nfact-1)
        for j in 1:nvars
            if fassign[j, 2] == i
                fnvars[i] = fnvars[i] + 1
            end
        end
    end

    # Track variables assigned to each 2nd-level factor
    # using a list of vectors called `varassign`
    # where the i-th entry is a vector
    # that contains the indices of the corresponding variables
    # in the dataset
    varassign = Any[]
    for i in 1:(nfact-1) # iterate over level-2 factors
        push!(varassign, Any[])
        for j in 1:nvar # iterate over variables
            if fassign[j, 2] == i
                push!(varassign[i], j)
            end
        end
    end

    # Load data
    ytemp = y

    # Express each variable in the dataset
    # in deviation-from-mean form
    # (de-mean each variable in the dataset)
    y = ytemp - repeat(mean(ytemp, dims = 1), capt, 1)

    # Set up some matrices for storage
    Xtsave = zeros(capt, nfact, (ndraws + burnin))          # Keep draw of factor, not all states (others are trivial)
    bsave = zeros(ndraws + burnin, nreg * nvar)             # Observable equation regression coefficients
    ssave = zeros(ndraws + burnin, nvar)                    # Innovation variances
    psave = zeros(ndraws + burnin, nfact * arlag)           # Factor autoregressive polynomials
    psave2 = zeros(ndraws + burnin, nvar * arterms)         # Idiosyncratic disturbance autoregressive polynomials

    ##############################################
    ##############################################
    ##### include a necessary procedure #####
    ##### set priors for observable equations in factor model

    # Specify prior mean of factor loading
    # and precision of factor loading
    b0_ = ones(nreg, 1)                 # prior mean of factor loading
    B0__ = 0.01 * ident(nreg)           # prior precision of factor loading
    B0__[1, 1] = 1.0

    # Specify prior mean of idiosyncratic
    # disturbance AR lag coefficients
    # and precision of AR lag coefficients
    r0_ = zeros(arterms, 1)             # prior mean of phi (idiosyncratic AR polynomial)
    phipri = 0.25
    R0__ = phipri * ident(arterms)      # prior precision of phi

    # Specify prior parameters of
    # innovation variances
    v0_ = trunc(Int, ceil(capt * 0.05)) # inverted gamma parameters of for innovation variances
    d0_ = (0.25)^2                      # v0 is an integer

    # Specify prior for factor AR lag coefficients
    # and precision of factor AR lag coefficients
    prem = 1
    for i = 1:(arlag-1)
        prem = [prem; (1 / 2)^i]
    end
    phiprif = 1 ./ prem
    r0f_ = zeros(arlag)                        # prior mean of phi
    R0f__ = diagrv(ident(arlag), phiprif)
    #R0f__ = phipri * ident(arlag)

    # Normalize innovation variance for the factor vector
    # since diagonal variance is set to 1
    sigU = [1; zeros(arlag - 1)]
    for i = 2:nfact
        sigU = [sigU; 1; zeros(arlag - 1)]
    end

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

    # Initialize hyperparameters
    SigE = ones(nvar, 1) * 0.0001           # Idiosyncratic error precision vector
    phi = zeros(arlag, nfact)               # Factor AR companion matrix
    bold = zeros(nvar, nreg)                # Obs. eq. regression starting coefficient matrix
    phimat0 = zeros(arterms, nvar)          # Idiosyncratic error AR companion matrix

    # Initialize factor series
    facts = rand(capt, nfact)           # Random starting factor series matrix
    facts[:, 1] = mean(y, dims = 2)     # Starting global factor = crosssectional mean of obs. series
    for i in 1:(nfact-1)              # Set level-2 factors equal to their respective group means
        facts[:, 1+i] = mean(y[:, varassign[i]], dims = 2)
    end

    # Begin Monte Carlo Loop
    # (start iteratively drawing hyperparameters and factors)
    for dr = 1:(ndraws+burnin)

        println(dr)

        ############################################
        ## Draw observation equation hyperparameters
        ############################################
        for i = 1:nvar # Iterate over all obs. variables

            # Save the index of the factor assigned
            # to observable variable i
            nf = fassign[i, 2]

            # Create matrix containing all regressors
            # corresponding to variable i including:
            # (1) an intercept, (2) global factor, (3) level-2 factor
            xft = [ones(capt, 1) facts[:, 1] facts[:, 1+nf]]

            # Update variable i observation equation
            # hyperparameters, and update corresponding
            # factor orientation (if appropriate)
            b1, s21, phi1, facts = ar_LJ(y[:, i], xft, arterms, b0_, B0__, r0_, R0__, v0_, d0_, transp_dbl(bold[i, :]), SigE[i], phimat0[:, i], i, nf, facts, capt, nreg, fnvars[fassign[i, 2]], varassign)
            #b1, s21, phi1, facts = ar2(y[:, i], xft, arterms, b0_, B0__, r0_, R0__, v0_, d0_, vec(bold[i, :]), phimat0[:, i], SigE[i], i, nf, facts, capt, varassign)
            bold[i, :] = b1                                     # Update obs. regression coefficients
            phimat0[:, i] = phi1                                # Update idiosyncratic error AR coefficients
            SigE[i] = s21                                       # Update innovation variance parameter
            bsave[dr, (((i-1)*nreg)+1):(i*nreg)] = b1           # Save current obs. regression coefficient draw
            ssave[dr, i] = s21                                  # Save current innovation variance parameter draw
            psave2[dr, (((i-1)*arterms)+1):(i*arterms)] = phi1  # Save current idiosyncratic error AR coefficient draw

        end

        ############################################
        ## Draw factor autoregression coefficients
        ############################################
        for i = 1:nfact # Iterate over all factors
            j = 1 + (i - 1) * arlag
            phi[:, i] = arfac(facts[:, i], arlag, r0f_, R0f__, phi[:, i], sigU[j, 1], capt)
            psave[dr, ((i-1)*arlag+1):((i-1)*arlag+arlag)] = transp_dbl(phi[:, i])
        end

        ############################################
        ## Draw global (level-1) factor
        ############################################

        # Initialize all important objects
        sinvf1 = sigbig(phi[:, 1], arlag, capt)             # (T×T) S^{-1} quasi-differencing matrix for global factor
        f = zeros(capt, 1)                                  # Empty vector for global factor to fill out
        H = ((1 / sigU[1]) * (transp_dbl(sinvf1) * sinvf1)) # First term of (T×T) H matrix, implying b_0 = 0 (H^{-1} is factor covariance matrix)

        # Fill out important objects
        for i = 1:nvar  # Iterate over all observable variable

            # Save level-2 factor index assigned to obs. variable i
            nfC = fassign[i, 2]

            # Partial out variation in variable i due to intercept + level-2 factor
            yW = y[:, i] - ones(capt, 1) * bold[i, 1] - facts[:, 1+nfC] * bold[i, 3]'

            # S_i^{-1} for i > 2
            sinv1 = sigbig(phimat0[:, i], arterms, capt)

            # Add next term in equation for H (pg. 1004, Otrok-Whiteman 1998)
            H = H + ((bold[i, 2]^2 / (SigE[i])) * (transp_dbl(sinv1) * sinv1))

            # Add next term of within-parenthesis sum in equation for f (pg. 1004, Otrok-Whiteman 1998)
            f = f + (bold[i, 2] / SigE[i]) * (transp_dbl(sinv1) * sinv1) * yW

        end
        Hinv = invpd(H)     # Invert H to save H^{-1}
        f = Hinv * f        # Obtain mean of f by pre-multiplying existing sum by H^{-1}

        # Obtain new draw of the global factor
        fact1 = sim_MvNormal(vec(f), Hinv)

        # Save draw in output object
        Xtsave[:, 1, dr] = fact1

        # Update factor data matrix to contain
        # new global factor draw
        facts[:, 1] = fact1

        ############################################
        ## Draw level-2 factors
        ############################################
        for c = 1:(nfact-1) # Iterate over level-2 factors

            j = 1 + c * arlag

            # Store number of obs. variables
            # to which level-2 factor number c
            # gets assigned
            Size = fnvars[c]

            # Partial out variation in variables assigned to level-2 factor c
            # corresponding to variation in intercept + level-1 factor
            yC = y[:, varassign[c]] - ones(capt, 1) * transp_dbl(bold[varassign[c], 1]) - facts[:, 1] * transp_dbl(bold[varassign[c], 2])

            # Initialize all important objects
            phiC = phi[:, 1+c]                                      # Store level-2 factor c AR lag coefficients
            sinvf1 = sigbig(phiC, arlag, capt)                      # (T×T) S^{-1} quasi-differencing matrix for level-2 factor c
            f = zeros(capt, 1)                                      # Empty vector for level-2 factor c to fill out
            H = ((1 / sigU[j]) * (transp_dbl(sinvf1) * sinvf1))     # First term of (T×T) H matrix, implying b_0 = 0 (H^{-1} is factor covariance matrix)

            for i = 1:Size # Iterate over all obs. variables assigned to level-2 factor c

                # S_i^{-1} for i > 2
                sinv1 = sigbig(phimat0[:, varassign[c][i]], arterms, capt)

                # Add next term in equation for H (pg. 1004, Otrok-Whiteman 1998)
                H = H + ((bold[varassign[c][i], 3]^2 / (SigE[varassign[c][i]])) * (transp_dbl(sinv1) * sinv1))

                # Add next term of within-parenthesis sum in equation for f (pg. 1004, Otrok-Whiteman 1998)
                f = f + (bold[varassign[c][i], 3] / SigE[varassign[c][i]]) * (transp_dbl(sinv1) * sinv1) * (yC[:, i])

            end
            Hinv = invpd(H)     # Invert H to save H^{-1}
            f = Hinv * f        # Obtain mean of f by pre-multiplying existing sum by H^{-1}

            # Obtain new draw of level-2 factor c
            fact2 = sim_MvNormal(vec(f), Hinv)

            # Save draw in output object
            Xtsave[:, 1+c, dr] = fact2

            # Update factor data matrix to contain
            # new level-2 factor c draw
            facts[:, 1+c] = fact2

        end

        println(dr)
    end

    ##############################################
    ##############################################
    ## Save Monte Carlo samples
    Xtsave = Xtsave[:, :, (burnin+1):(burnin+ndraws)]
    bsave = bsave[(burnin+1):(burnin+ndraws), :]
    ssave = ssave[(burnin+1):(burnin+ndraws), :]
    psave = psave[(burnin+1):(burnin+ndraws), :]
    psave2 = psave2[(burnin+1):(burnin+ndraws), :]

    ##############################################
    ##############################################
    ## Save Monte Carlo means
    F = mean(Xtsave, dims = 3)
    B = mean(bsave, dims = 1)
    S = mean(ssave, dims = 1)
    P = mean(psave, dims = 1)
    P2 = mean(psave2, dims = 1)

    ##############################################
    ##############################################
    means = DFMMeans(F, B, S, P, P2)
    results = DFMResults(Xtsave, bsave, ssave, psave, psave2, means)

    return results
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb03vd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 20
    PMAX = 20
    LDA1 = NMAX
    LDA2 = NMAX
    LDQ1 = NMAX
    LDQ2 = NMAX
    LDTAU = NMAX-1
    LDWORK = NMAX
    ZERO = 0.0e0
    ONE = 1.0e0
    A = Array{Float64,3}(undef, LDA1,LDA2,PMAX)
    TAU = Array{Float64,2}(undef, LDTAU,PMAX)
    Q = Array{Float64,3}(undef, LDQ1,LDQ2,PMAX)
    AS = Array{Float64,3}(undef, LDA1,LDA2,PMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    QTA = Array{Float64,2}(undef, LDQ1,NMAX)
    SSQ = 0.0
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    P = parse(BlasInt, vs[2])
    ILO = parse(BlasInt, vs[3])
    IHI = parse(BlasInt, vs[4])
    if ( N<0 || N>min( LDA1, LDA2 ) )
        @error "illegal N=$N"
    end
    if ( P<=0 || P>PMAX )
        @error "illegal P=$P"
    end

    _isz,_jsz,_ksz = (N,N,P)
    for k in 1:_ksz
      vs = String[]
      while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
      end
      for i in 1:_isz
        _i0 = (i-1)*_jsz
        A[i,1:_jsz,k] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
      end
    end
# interp call 1

    #FOREIGN.dlacpy!( 'F', N, N, A(1,1,K), LDA1, AS(1,1,K), LDA1 )
    AS = copy(A)

# interp call 2

    INFO = SLICOT.mb03vd!(N, P, ILO, IHI, A, TAU)
    @test INFO == 0
    INFO == 0 || return

    if ( INFO!=0 )
        @warn "mb03vd returns info=$INFO"
    end
# interp call 3

  for K in 1:P
    # FOREIGN.dlacpy!( 'L', N, N, A(1,1,K), LDA1, Q(1,1,K), LDQ1 )
        Q[1:N,1:N,K] .= tril(A[1:N,1:N,K])

      if ( N>1 )
          if ( N>2 && K==1 )
# interp call 4

              # FOREIGN.dlaset!( 'L', N-2, N-2, ZERO, ZERO, A(3,1,K), LDA1 )
              triu!(view(A,1:N,1:N,K),-1)

          elseif ( K>1 )
# interp call 5

              # FOREIGN.dlaset!( 'L', N-1, N-1, ZERO, ZERO, A(2,1,K), LDA1 )
              triu!(view(A,1:N,1:N,K))
          end # if
      end # if
  end # for K

    println(io, "A:")
    _nc = N
    _nr = N
    _nk = P
    show(io, "text/plain", A[1:_nr,1:_nc,1:_nk])
    println(io)

# interp call 6

    INFO = SLICOT.mb03vy!(N, P, ILO, IHI, Q, TAU)
    @test INFO == 0
    INFO == 0 || return
    println(io, "Q:")
    _nc = N
    _nr = N
    _nk = P
    show(io, "text/plain", Q[1:_nr,1:_nc,1:_nk])
    println(io)

    if ( INFO!=0 )
        @warn "mb03vy returns info=$INFO"
    else
    SSQ = ZERO
    for K in 1:P
        KP1 = mod(K,P)+1
# interp call 7

    #FOREIGN.dgemm!( 'T', 'N', N, N, N, ONE, Q(1,1,K), LDQ1, AS(1,1,K), LDA1, ZERO, QTA, LDQ1 )
        qta = Q[1:N,1:N,K]' * AS[1:N,1:N,K]
# interp call 8

    #FOREIGN.dgemm!( 'N', 'N', N, N, N, ONE, QTA, LDQ1, Q(1,1,KP1), LDQ1, -ONE, A(1,1,K), LDA1 )
        mul!(view(A,1:N,1:N,K),qta,view(Q,1:N,1:N,KP1),1,-1)
        #SSQ = DLAPY2( SSQ, DLANGE( 'Frobenius', N, N, A(1,1,K), LDA1, DWORK ) )
        SSQ = hypot(SSQ, norm(view(A,1:N,1:N,K)))
    end
    println(io, "norm(Q'*A*Q - Aout) = ", SSQ)
    end
    close(f)
end # run_mb03vd()


text:
function save_data!(A::Array{T,1}, dset::HDF5Dataset, src_ind::AbstractRange{Int}, dest_ind::AbstractRange{Int}) where T
  dsel_id = HDF5.hyperslab(dset, src_ind)
  V = view(A, dest_ind)
  memtype = HDF5.datatype(A)
  memspace = HDF5.dataspace(V)
  HDF5.h5d_write(dset.id, memtype.id, memspace.id, dsel_id, dset.xfer, V)
  HDF5.close(memtype)
  HDF5.close(memspace)
  HDF5.h5s_close(dsel_id)
  return nothing
end