text:
# Nest
#
# If the line exceeds the print width it will be nested.
#
# This is done by replacing `PLACEHOLDER` nodes with `NEWLINE`
# nodes and updating the FST's indent.
#
# `extra_margin` provides additional width to consider from
# the top-level node.
#
# Example:
#
#     LHS op RHS
#
# the length of " op" will be considered when nesting LHS

function nest!(
    ds::DefaultStyle,
    nodes::Vector{FST},
    s::State,
    indent::Int;
    extra_margin = 0,
)
    style = getstyle(ds)
    for (i, n) in enumerate(nodes)
        if n.typ === NEWLINE && nodes[i+1].typ === Block
            s.line_offset = nodes[i+1].indent
        elseif n.typ === NOTCODE && nodes[i+1].typ === Block
            s.line_offset = nodes[i+1].indent
        elseif n.typ === NEWLINE
            s.line_offset = indent
        else
            n.extra_margin = extra_margin
            nest!(style, n, s)
        end
    end
end
nest!(
    style::S,
    nodes::Vector{FST},
    s::State,
    indent::Int;
    extra_margin = 0,
) where {S<:AbstractStyle} =
    nest!(DefaultStyle(style), nodes, s, indent, extra_margin = extra_margin)

function nest!(ds::DefaultStyle, fst::FST, s::State)
    style = getstyle(ds)
    if is_leaf(fst)
        s.line_offset += length(fst)
        return
    end

    if cant_nest(fst)
        walk(increment_line_offset!, fst, s)
        return
    end

    if fst.typ === Import
        n_import!(style, fst, s)
    elseif fst.typ === Export
        n_export!(style, fst, s)
    elseif fst.typ === Using
        n_using!(style, fst, s)
    elseif fst.typ === Where
        n_whereopcall!(style, fst, s)
    elseif fst.typ === Conditional
        line_margin = s.line_offset + length(fst) + fst.extra_margin
        if s.opts.conditional_to_if && line_margin > s.opts.margin
            conditional_to_if_block!(fst, s)
            nest!(style, fst, s)
        else
            n_conditionalopcall!(style, fst, s)
        end
    elseif fst.typ === Binary
        line_margin = s.line_offset + length(fst) + fst.extra_margin
        if s.opts.short_to_long_function_def &&
           line_margin > s.opts.margin &&
           fst.ref !== nothing &&
           CSTParser.defines_function(fst.ref[])
            short_to_long_function_def!(fst, s)
        end
        if fst.typ === Binary
            n_binaryopcall!(style, fst, s)
        else
            nest!(style, fst, s)
        end
    elseif fst.typ === Curly
        n_curly!(style, fst, s)
    elseif fst.typ === Call
        n_call!(style, fst, s)
    elseif fst.typ === MacroCall
        n_macrocall!(style, fst, s)
    elseif fst.typ === RefN
        n_ref!(style, fst, s)
    elseif fst.typ === TypedVcat
        n_typedvcat!(style, fst, s)
    elseif fst.typ === TupleN
        n_tuple!(style, fst, s)
    elseif fst.typ === Vect
        n_vect!(style, fst, s)
    elseif fst.typ === Vcat
        n_vcat!(style, fst, s)
    elseif fst.typ === Braces
        n_braces!(style, fst, s)
    elseif fst.typ === BracesCat
        n_bracescat!(style, fst, s)
    elseif fst.typ === Brackets
        n_invisbrackets!(style, fst, s)
    elseif fst.typ === Comprehension
        n_comprehension!(style, fst, s)
    elseif fst.typ === TypedComprehension
        n_typedcomprehension!(style, fst, s)
    elseif fst.typ === Do
        n_do!(style, fst, s)
    elseif fst.typ === Generator
        n_generator!(style, fst, s)
    elseif fst.typ === Filter
        n_filter!(style, fst, s)
    elseif fst.typ === Flatten
        n_flatten!(style, fst, s)
    elseif fst.typ === Block
        n_block!(style, fst, s)
    elseif fst.typ === Chain
        n_chainopcall!(style, fst, s)
    elseif fst.typ === Comparison
        n_comparison!(style, fst, s)
    elseif fst.typ === For
        n_for!(style, fst, s)
    elseif fst.typ === Let
        n_let!(style, fst, s)
    elseif fst.typ === Unary && length(fst.nodes) > 1 && fst[2].typ === OPERATOR
        n_unaryopcall!(style, fst, s)
    elseif fst.typ === StringN
        n_string!(style, fst, s)
    else
        nest!(style, fst.nodes, s, fst.indent, extra_margin = fst.extra_margin)
    end
end
nest!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    nest!(DefaultStyle(style), fst, s)

function n_string!(ds::DefaultStyle, fst::FST, s::State)
    style = getstyle(ds)
    # difference in positioning of the string
    # from the source document to the formatted document
    diff = s.line_offset - fst.indent

    for (i, n) in enumerate(fst.nodes)
        if n.typ === NEWLINE
            s.line_offset = fst[i+1].indent + diff
        else
            nest!(style, n, s)
        end
    end
end
n_string!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_string!(DefaultStyle(style), fst, s)

function n_unaryopcall!(ds::DefaultStyle, fst::FST, s::State)
    style = getstyle(ds)
    fst[1].extra_margin = fst.extra_margin + length(fst[2])
    nest!(style, fst[1], s)
    nest!(style, fst[2], s)
end
n_unaryopcall!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_unaryopcall!(DefaultStyle(style), fst, s)

function n_do!(ds::DefaultStyle, fst::FST, s::State)
    style = getstyle(ds)
    extra_margin = sum(length.(fst[2:3]))
    # make sure there are nodes after "do"
    if fst[4].typ === WHITESPACE
        extra_margin += length(fst[4])
        extra_margin += length(fst[5])
    end
    fst[1].extra_margin = fst.extra_margin + extra_margin
    nest!(style, fst[1], s)
    nest!(style, fst[2:end], s, fst.indent, extra_margin = fst.extra_margin)
end
n_do!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_do!(DefaultStyle(style), fst, s)

# Import,Using,Export
function n_using!(ds::DefaultStyle, fst::FST, s::State)
    style = getstyle(ds)
    line_margin = s.line_offset + length(fst) + fst.extra_margin
    idx = findfirst(n -> n.typ === PLACEHOLDER, fst.nodes)
    if idx !== nothing && (line_margin > s.opts.margin || must_nest(fst))
        fst.indent += s.opts.indent

        if can_nest(fst)
            if fst.indent + sum(length.(fst[idx+1:end])) <= s.opts.margin
                fst[idx] = Newline(length = fst[idx].len)
                walk(increment_line_offset!, fst, s)
                return
            end
        end

        for (i, n) in enumerate(fst.nodes)
            if n.typ === NEWLINE
                s.line_offset = fst.indent
            elseif n.typ === PLACEHOLDER
                fst[i] = Newline(length = n.len)
                s.line_offset = fst.indent
            else
                diff = fst.indent - fst[i].indent
                add_indent!(n, s, diff)
                i < length(fst.nodes) && (n.extra_margin = 1)
                nest!(style, n, s)
            end
        end
    else
        nest!(style, fst.nodes, s, fst.indent, extra_margin = fst.extra_margin)
    end
end
n_using!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_using!(DefaultStyle(style), fst, s)

@inline n_export!(ds::DefaultStyle, fst::FST, s::State) = n_using!(ds, fst, s)
n_export!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_export!(DefaultStyle(style), fst, s)

@inline n_import!(ds::DefaultStyle, fst::FST, s::State) = n_using!(ds, fst, s)
n_import!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_import!(DefaultStyle(style), fst, s)

function n_tuple!(ds::DefaultStyle, fst::FST, s::State)
    style = getstyle(ds)
    line_margin = s.line_offset + length(fst) + fst.extra_margin
    idx = findlast(n -> n.typ === PLACEHOLDER, fst.nodes)
    opener = findfirst(is_opener, fst.nodes) !== nothing

    if idx !== nothing && (line_margin > s.opts.margin || must_nest(fst))
        line_offset = s.line_offset
        if opener
            fst[end].indent = fst.indent
        end
        if fst.typ !== TupleN || opener
            fst.indent += s.opts.indent
        end

        for (i, n) in enumerate(fst.nodes)
            if n.typ === NEWLINE
                s.line_offset = fst.indent
            elseif n.typ === PLACEHOLDER
                fst[i] = Newline(length = n.len)
                s.line_offset = fst.indent
            elseif n.typ === TRAILINGCOMMA
                n.val = ","
                n.len = 1
                nest!(style, n, s)
            elseif n.typ === TRAILINGSEMICOLON
                n.val = ";"
                n.len = 1
                nest!(style, n, s)
            elseif n.typ === INVERSETRAILINGSEMICOLON
                n.val = ""
                n.len = 0
                nest!(style, n, s)
            elseif opener && (i == 1 || i == length(fst.nodes))
                nest!(style, n, s)
            else
                diff = fst.indent - fst[i].indent
                add_indent!(n, s, diff)
                n.extra_margin = 1
                nest!(style, n, s)
            end
        end

        if opener
            s.line_offset = fst[end].indent + 1
        end
    else
        extra_margin = fst.extra_margin
        opener && (extra_margin += 1)
        nest!(style, fst.nodes, s, fst.indent, extra_margin = extra_margin)
    end
end
n_tuple!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_tuple!(DefaultStyle(style), fst, s)

@inline n_vect!(ds::DefaultStyle, fst::FST, s::State) = n_tuple!(ds, fst, s)
n_vect!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_vect!(DefaultStyle(style), fst, s)

@inline n_vcat!(ds::DefaultStyle, fst::FST, s::State) = n_tuple!(ds, fst, s)
n_vcat!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_vcat!(DefaultStyle(style), fst, s)

@inline n_braces!(ds::DefaultStyle, fst::FST, s::State) = n_tuple!(ds, fst, s)
n_braces!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_braces!(DefaultStyle(style), fst, s)

@inline n_bracescat!(ds::DefaultStyle, fst::FST, s::State) = n_tuple!(ds, fst, s)
n_bracescat!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_bracescat!(DefaultStyle(style), fst, s)

@inline n_parameters!(ds::DefaultStyle, fst::FST, s::State) = n_tuple!(ds, fst, s)
n_parameters!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_parameters!(DefaultStyle(style), fst, s)

@inline n_invisbrackets!(ds::DefaultStyle, fst::FST, s::State) = n_tuple!(ds, fst, s)
n_invisbrackets!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_invisbrackets!(DefaultStyle(style), fst, s)

@inline n_call!(ds::DefaultStyle, fst::FST, s::State) = n_tuple!(ds, fst, s)
n_call!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_call!(DefaultStyle(style), fst, s)

@inline n_curly!(ds::DefaultStyle, fst::FST, s::State) = n_tuple!(ds, fst, s)
n_curly!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_curly!(DefaultStyle(style), fst, s)

@inline n_macrocall!(ds::DefaultStyle, fst::FST, s::State) = n_tuple!(ds, fst, s)
n_macrocall!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_macrocall!(DefaultStyle(style), fst, s)

@inline n_ref!(ds::DefaultStyle, fst::FST, s::State) = n_tuple!(ds, fst, s)
n_ref!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_ref!(DefaultStyle(style), fst, s)

@inline n_typedvcat!(ds::DefaultStyle, fst::FST, s::State) = n_tuple!(ds, fst, s)
n_typedvcat!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_typedvcat!(DefaultStyle(style), fst, s)

function n_comprehension!(ds::DefaultStyle, fst::FST, s::State; indent = -1)
    style = getstyle(ds)
    line_offset = s.line_offset
    line_margin = s.line_offset + length(fst) + fst.extra_margin

    closer = is_closer(fst[end])
    if closer && (line_margin > s.opts.margin || must_nest(fst))
        idx = findfirst(n -> n.typ === PLACEHOLDER, fst.nodes)
        if idx !== nothing
            fst[idx] = Newline(length = fst[idx].len)
        end
        idx = findlast(n -> n.typ === PLACEHOLDER, fst.nodes)
        if idx !== nothing
            fst[idx] = Newline(length = fst[idx].len)
        end

        add_indent!(fst, s, s.opts.indent)
        fst[end].indent = fst.indent - s.opts.indent
    end

    for (i, n) in enumerate(fst.nodes)
        if n.typ === NEWLINE
            s.line_offset = fst.indent
        elseif n.typ === PLACEHOLDER
            if must_nest(fst)
                fst[i] = Newline(length = n.len)
                s.line_offset = fst.indent
            else
                nest_if_over_margin!(style, fst, s, i)
            end
        elseif i == length(fst.nodes) && !closer
            nest!(style, n, s)
        else
            n.extra_margin = 1
            nest!(style, n, s)
        end
    end

    if closer
        s.line_offset = fst[end].indent + 1
    end
end

n_comprehension!(style::S, fst::FST, s::State; indent = -1) where {S<:AbstractStyle} =
    n_comprehension!(DefaultStyle(style), fst, s, indent = indent)

@inline n_typedcomprehension!(ds::DefaultStyle, fst::FST, s::State) =
    n_comprehension!(ds, fst, s)
n_typedcomprehension!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_typedcomprehension!(DefaultStyle(style), fst, s)

function n_generator!(ds::DefaultStyle, fst::FST, s::State; indent = -1)
    style = getstyle(ds)
    line_margin = s.line_offset + length(fst) + fst.extra_margin

    if line_margin > s.opts.margin || must_nest(fst)
        line_offset = s.line_offset
        fst.indent = s.line_offset
        phs = reverse(findall(n -> n.typ === PLACEHOLDER, fst.nodes))
        for (i, idx) in enumerate(phs)
            if i == 1
                fst[idx] = Newline(length = fst[idx].len)
            else
                nidx = phs[i-1]
                l1 = sum(length.(fst[1:idx-1]))
                l2 = sum(length.(fst[idx:nidx-1]))
                width = line_offset + l1 + l2
                if must_nest(fst) || width > s.opts.margin
                    fst[idx] = Newline(length = fst[idx].len)
                end
            end
        end

        for (i, n) in enumerate(fst.nodes)
            if n.typ === NEWLINE
                s.line_offset = fst.indent
            elseif i == length(fst.nodes)
                nest!(style, n, s)
            else
                n.extra_margin = 1
                nest!(style, n, s)
            end
        end

        s.line_offset = line_offset
        for (i, n) in enumerate(fst.nodes)
            if n.typ === NEWLINE && !is_comment(fst[i+1]) && !is_comment(fst[i-1])
                # +1 for newline to whitespace conversion
                width = s.line_offset + 1
                w, _ = length_to(fst, (NEWLINE,), start = i + 1)
                width += w
                # @info "" s.line_offset w width fst[i-1].typ fst[i-1].val width <= s.opts.margin
                if width <= s.opts.margin
                    fst[i] = Whitespace(1)
                    s.line_offset += length(fst[i])
                else
                    s.line_offset = fst.indent
                end
            else
                walk(increment_line_offset!, fst[i], s)
            end
        end

        s.line_offset = line_offset
        walk(increment_line_offset!, fst, s)
    else
        nest!(style, fst.nodes, s, fst.indent, extra_margin = fst.extra_margin)
    end
end

n_generator!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_generator!(DefaultStyle(style), fst, s)

@inline n_filter!(ds::DefaultStyle, fst::FST, s::State) =
    n_generator!(ds, fst, s, indent = fst.indent)
n_filter!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_filter!(DefaultStyle(style), fst, s)

@inline n_flatten!(ds::DefaultStyle, fst::FST, s::State) =
    n_generator!(ds, fst, s, indent = fst.indent)
n_flatten!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_flatten!(DefaultStyle(style), fst, s)

function n_whereopcall!(ds::DefaultStyle, fst::FST, s::State)
    style = getstyle(ds)
    line_margin = s.line_offset + length(fst) + fst.extra_margin
    if line_margin > s.opts.margin || must_nest(fst)
        line_offset = s.line_offset
        Blen = sum(length.(fst[2:end]))

        fst[1].extra_margin = Blen + fst.extra_margin
        nest!(style, fst[1], s)

        # "B"
        has_braces = is_closer(fst[end])
        if has_braces
            fst[end].indent = fst.indent
        end

        over = (s.line_offset + Blen + fst.extra_margin > s.opts.margin) || must_nest(fst)
        fst.indent += s.opts.indent
        for (i, n) in enumerate(fst[2:end])
            if n.typ === NEWLINE
                s.line_offset = fst.indent
            elseif is_opener(n)
                if fst.indent - s.line_offset > 1
                    fst.indent = s.line_offset + 1
                    fst[end].indent = s.line_offset
                end
                nest!(style, n, s)
            elseif n.typ === PLACEHOLDER && over
                fst[i+1] = Newline(length = length(n))
                s.line_offset = fst.indent
            elseif n.typ === TRAILINGCOMMA && over
                n.val = ","
                n.len = 1
                nest!(style, n, s)
            elseif has_braces
                n.extra_margin = 1 + fst.extra_margin
                nest!(style, n, s)
            else
                n.extra_margin = fst.extra_margin
                nest!(style, n, s)
            end
        end

        s.line_offset = line_offset
        walk(increment_line_offset!, fst, s)
    else
        nest!(style, fst.nodes, s, fst.indent, extra_margin = fst.extra_margin)
    end
end
n_whereopcall!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_whereopcall!(DefaultStyle(style), fst, s)

function n_conditionalopcall!(ds::DefaultStyle, fst::FST, s::State)
    style = getstyle(ds)
    line_margin = s.line_offset + length(fst) + fst.extra_margin
    if line_margin > s.opts.margin || must_nest(fst)
        line_offset = s.line_offset
        fst.indent = s.line_offset
        phs = reverse(findall(n -> n.typ === PLACEHOLDER, fst.nodes))
        for (i, idx) in enumerate(phs)
            if i == 1
                fst[idx] = Newline(length = fst[idx].len)
            else
                nidx = phs[i-1]
                l1 = sum(length.(fst[1:idx-1]))
                l2 = sum(length.(fst[idx:nidx-1]))
                width = line_offset + l1 + l2
                if must_nest(fst) || width > s.opts.margin
                    fst[idx] = Newline(length = fst[idx].len)
                end
            end
        end

        for (i, n) in enumerate(fst.nodes)
            if i == length(fst.nodes)
                n.extra_margin = fst.extra_margin
                nest!(style, n, s)
            elseif fst[i+1].typ === WHITESPACE
                n.extra_margin = length(fst[i+1]) + length(fst[i+2])
                nest!(style, n, s)
            elseif n.typ === NEWLINE
                s.line_offset = fst.indent
            else
                nest!(style, n, s)
            end
        end

        s.line_offset = line_offset
        for (i, n) in enumerate(fst.nodes)
            if n.typ === NEWLINE && !is_comment(fst[i+1]) && !is_comment(fst[i-1])
                # +1 for newline to whitespace conversion
                width = s.line_offset + 1
                if i == length(fst.nodes) - 1
                    width += sum(length.(fst[i+1:end])) + fst.extra_margin
                else
                    width += sum(length.(fst[i+1:i+3]))
                end
                if width <= s.opts.margin
                    fst[i] = Whitespace(1)
                    s.line_offset += length(fst[i])
                else
                    s.line_offset = fst.indent
                end
            else
                walk(increment_line_offset!, fst[i], s)
            end
        end

        s.line_offset = line_offset
        walk(increment_line_offset!, fst, s)
    else
        nest!(style, fst.nodes, s, fst.indent, extra_margin = fst.extra_margin)
    end
end
n_conditionalopcall!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_conditionalopcall!(DefaultStyle(style), fst, s)

function no_unnest(fst::FST)
    if fst.typ === Binary ||
       fst.typ === Conditional ||
       fst.typ === Chain ||
       fst.typ === Comparison
        return contains_comment(fst)
    end
    return false
end

function n_binaryopcall!(ds::DefaultStyle, fst::FST, s::State)
    style = getstyle(ds)
    line_margin = s.line_offset + length(fst) + fst.extra_margin
    # If there's no placeholder the binary call is not nestable
    idxs = findall(n -> n.typ === PLACEHOLDER, fst.nodes)
    rhs = fst[end]
    rhs.typ === Block && (rhs = rhs[1])

    if length(idxs) == 2 &&
       (line_margin > s.opts.margin || must_nest(fst) || must_nest(rhs))
        line_offset = s.line_offset
        i1 = idxs[1]
        i2 = idxs[2]
        fst[i1] = Newline(length = fst[i1].len)
        cst = fst.ref[]

        indent_nest =
            CSTParser.defines_function(cst) ||
            is_assignment(cst) ||
            op_kind(fst) === Tokens.PAIR_ARROW ||
            op_kind(fst) === Tokens.ANON_FUNC ||
            is_standalone_shortcircuit(cst)

        if indent_nest
            s.line_offset = fst.indent + s.opts.indent
            fst[i2] = Whitespace(s.opts.indent)

            # reset the indent of the RHS
            if fst[end].indent > fst.indent
                fst[end].indent = fst.indent
            end
            add_indent!(fst[end], s, s.opts.indent)
        else
            fst.indent = s.line_offset
        end

        # rhs
        fst[end].extra_margin = fst.extra_margin
        nest!(style, fst[end], s)

        # "lhs op" rhs
        s.line_offset = line_offset

        # extra margin for " op"
        fst[1].extra_margin = length(fst[2]) + length(fst[3])
        nest!(style, fst[1], s)
        for n in fst[2:i1]
            nest!(style, n, s)
        end

        # Undo nest if possible
        if can_nest(fst) && !no_unnest(rhs)
            line_margin = s.line_offset

            # replace IN with all of precedence level 6
            if (
                   rhs.typ === Binary &&
                   !(op_kind(rhs) === Tokens.IN || op_kind(rhs) === Tokens.DECLARATION)
               ) ||
               rhs.typ === Unary ||
               rhs.typ === Chain ||
               rhs.typ === Comparison ||
               rhs.typ === Conditional
                line_margin += length(fst[end])
            elseif rhs.typ === Do && is_iterable(rhs[1])
                rw, _ = length_to(fst, (NEWLINE,), start = i2 + 1)
                line_margin += rw
            elseif is_block(rhs)
                idx = findfirst(n -> n.typ === NEWLINE, rhs.nodes)
                if idx === nothing
                    line_margin += length(fst[end])
                else
                    line_margin += sum(length.(rhs[1:idx-1]))
                end
            else
                rw, _ = length_to(fst, (NEWLINE,), start = i2 + 1)
                line_margin += rw
            end

            if line_margin + fst.extra_margin <= s.opts.margin
                fst[i1] = Whitespace(1)
                if indent_nest
                    fst[i2] = Whitespace(0)
                    walk(unnest!, rhs, s)
                end
            end
        end

        s.line_offset = line_offset
        walk(increment_line_offset!, fst, s)
    else
        # length of op and surrounding whitespace
        oplen = sum(length.(fst[2:end]))

        for (i, n) in enumerate(fst.nodes)
            if n.typ === NEWLINE
                s.line_offset = fst.indent
            elseif i == 1
                n.extra_margin = oplen + fst.extra_margin
                nest!(style, n, s)
            elseif i == length(fst.nodes)
                n.extra_margin = fst.extra_margin
                nest!(style, n, s)
            else
                nest!(style, n, s)
            end
        end
    end
end
n_binaryopcall!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_binaryopcall!(DefaultStyle(style), fst, s)

function n_for!(ds::DefaultStyle, fst::FST, s::State)
    style = getstyle(ds)
    block_idx = findfirst(n -> !is_leaf(n), fst.nodes)
    if block_idx === nothing || length(fst[block_idx]) == 0
        nest!(style, fst.nodes, s, fst.indent, extra_margin = fst.extra_margin)
        return
    end
    nest!(style, fst.nodes, s, fst.indent, extra_margin = fst.extra_margin)
end
n_for!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_for!(DefaultStyle(style), fst, s)

@inline n_let!(ds::DefaultStyle, fst::FST, s::State) = n_for!(ds, fst, s)
n_let!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_let!(DefaultStyle(style), fst, s)

function n_block!(ds::DefaultStyle, fst::FST, s::State; indent = -1)
    style = getstyle(ds)
    line_margin = s.line_offset + length(fst) + fst.extra_margin
    idx = findfirst(n -> n.typ === PLACEHOLDER, fst.nodes)
    if idx !== nothing && (line_margin > s.opts.margin || must_nest(fst))
        line_offset = s.line_offset
        indent >= 0 && (fst.indent = indent)

        if fst.typ === Chain && is_standalone_shortcircuit(fst.ref[])
            fst.indent += s.opts.indent
        end

        for (i, n) in enumerate(fst.nodes)
            if n.typ === NEWLINE
                s.line_offset = fst.indent
            elseif n.typ === PLACEHOLDER
                fst[i] = Newline(length = n.len)
                s.line_offset = fst.indent
            elseif n.typ === TRAILINGCOMMA
                n.val = ","
                n.len = 1
                nest!(style, n, s)
            elseif i < length(fst.nodes) - 1 && fst[i+2].typ === OPERATOR
                # chainopcall / comparison
                diff = fst.indent - fst[i].indent
                add_indent!(n, s, diff)
                n.extra_margin = 1 + length(fst[i+2])
                nest!(style, n, s)
            else
                diff = fst.indent - fst[i].indent
                add_indent!(n, s, diff)
                n.extra_margin = 1
                nest!(style, n, s)
            end
        end
    else
        nest!(style, fst.nodes, s, fst.indent, extra_margin = fst.extra_margin)
    end
end
n_block!(style::S, fst::FST, s::State; indent = -1) where {S<:AbstractStyle} =
    n_block!(DefaultStyle(style), fst, s, indent = indent)

@inline n_comparison!(ds::DefaultStyle, fst::FST, s::State) =
    n_block!(ds, fst, s, indent = s.line_offset)
n_comparison!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_comparison!(DefaultStyle(style), fst, s)

@inline n_chainopcall!(ds::DefaultStyle, fst::FST, s::State) =
    n_block!(ds, fst, s, indent = s.line_offset)
n_chainopcall!(style::S, fst::FST, s::State) where {S<:AbstractStyle} =
    n_chainopcall!(DefaultStyle(style), fst, s)


text:
Base.transpose(a::AbstractJaxVector) = numpy.reshape(a, (-1,1))
Base.transpose(a::AbstractJaxMatrix) = numpy.transpose(a)

Base.adjoint(a::AbstractJaxArray{T}) where T<:Real = transpose(a)
Base.adjoint(a::AbstractJaxArray) = conj(transpose(a))

Base.PermutedDimsArray(a::AbstractJaxArray, perm) = permutedims(a, perm)
Base.permutedims(a::AbstractJaxArray) = transpose(a)
Base.permutedims(a::AbstractJaxArray{T,N}, perm::Dims{N}) where {T,N} =
    numpy.transpose(a, reverse(N.-perm))

Base.reshape(a::AbstractJaxArray, dims::Dims) = numpy.reshape(a, reverse(dims))

##
Base.similar(a::AbstractJaxArray, ::Type{T}, dims::Dims{N}) where {T,N} = JaxArray(T,dims)

##
Base.vcat(arrays::AbstractJaxVector...) = numpy.hstack(arrays)
Base.hcat(arrays::AbstractJaxVector...) = numpy.vstack(arrays)

ext: jl


text:
# A simple grid world MDP
# All cells with reward are also terminal

using CommonRLInterface
using StaticArrays
using Compose
import ColorSchemes

const RL = CommonRLInterface

mutable struct GridWorld <: AbstractEnv
    size::SVector{2, Int}
    rewards::Dict{SVector{2, Int}, Float64}
    state::SVector{2, Int}
end

function GridWorld()
    rewards = Dict(SA[9,3]=> 10.0,
                   SA[8,8]=>  3.0,
                   SA[4,3]=>-10.0,
                   SA[4,6]=> -5.0)
    return GridWorld(SA[10, 10], rewards, SA[rand(1:10), rand(1:10)])
end

RL.reset!(env::GridWorld) = (env.state = SA[rand(1:env.size[1]), rand(1:env.size[2])])
RL.actions(env::GridWorld) = (SA[1,0], SA[-1,0], SA[0,1], SA[0,-1])
RL.observe(env::GridWorld) = env.state
RL.terminated(env::GridWorld) = haskey(env.rewards, env.state)

function RL.act!(env::GridWorld, a)
    if rand() < 0.4 # 40% chance of going in a random direction (=30% chance of going in a wrong direction)
        a = rand(actions(env))
    end

    env.state = clamp.(env.state + a, SA[1,1], env.size)

    return get(env.rewards, env.state, 0.0)
end

# optional functions
@provide RL.observations(env::GridWorld) = [SA[x, y] for x in 1:env.size[1], y in 1:env.size[2]]
@provide RL.clone(env::GridWorld) = GridWorld(env.size, copy(env.rewards), env.state)
@provide RL.state(env::GridWorld) = env.state
@provide RL.setstate!(env::GridWorld, s) = (env.state = s)

@provide function RL.render(env::GridWorld)
    nx, ny = env.size
    cells = []
    for s in observations(env)
        r = get(env.rewards, s, 0.0)
        clr = get(ColorSchemes.redgreensplit, (r+10.0)/20.0)
        cell = context((s[1]-1)/nx, (ny-s[2])/ny, 1/nx, 1/ny)
        compose!(cell, rectangle(), fill(clr), stroke("gray"))
        push!(cells, cell)
    end
    grid = compose(context(), linewidth(0.5mm), cells...)
    outline = compose(context(), linewidth(1mm), rectangle(), stroke("gray"))

    s = env.state
    agent_ctx = context((s[1]-1)/nx, (ny-s[2])/ny, 1/nx, 1/ny)
    agent = compose(agent_ctx, circle(0.5, 0.5, 0.4), fill("orange"))

    sz = min(w,h)
    return compose(context((w-sz)/2, (h-sz)/2, sz, sz), agent, grid, outline)
end


text:
@sumProductRule(:node_type     => GaussianMeanVariance,
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Nothing, Message{PointMass}, Message{PointMass}),
                :name          => SPGaussianMeanVarianceOutNPP)

@sumProductRule(:node_type     => GaussianMeanVariance,
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Message{PointMass}, Nothing, Message{PointMass}),
                :name          => SPGaussianMeanVarianceMPNP)

@sumProductRule(:node_type     => GaussianMeanVariance,
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Nothing, Message{Gaussian}, Message{PointMass}),
                :name          => SPGaussianMeanVarianceOutNGP)

@sumProductRule(:node_type     => GaussianMeanVariance,
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Message{Gaussian}, Nothing, Message{PointMass}),
                :name          => SPGaussianMeanVarianceMGNP)

@sumProductRule(:node_type     => GaussianMeanVariance,
                :outbound_type => Message{Function},
                :inbound_types => (Message{Gaussian}, Message{Gaussian}, Nothing),
                :name          => SPGaussianMeanVarianceVGGN)

@sumProductRule(:node_type     => GaussianMeanVariance,
                :outbound_type => Message{Function},
                :inbound_types => (Message{PointMass}, Message{Gaussian}, Nothing),
                :name          => SPGaussianMeanVarianceVPGN)

@sumProductRule(:node_type     => GaussianMeanVariance,
                :outbound_type => Message{SampleList},
                :inbound_types => (Nothing, Message{SampleList}, Message{PointMass}),
                :name          => SPGaussianMeanVarianceOutNSP)

@sumProductRule(:node_type     => GaussianMeanVariance,
                :outbound_type => Message{SampleList},
                :inbound_types => (Message{SampleList},Nothing,Message{PointMass}),
                :name          => SPGaussianMeanVarianceMSNP)

@sumProductRule(:node_type     => GaussianMeanVariance,
                :outbound_type => Message{SampleList},
                :inbound_types => (Nothing, Message{Gaussian},Message{SampleList}),
                :name          => SPGaussianMeanVarianceOutNGS)

@sumProductRule(:node_type     => GaussianMeanVariance,
                :outbound_type => Message{SampleList},
                :inbound_types => (Message{Gaussian}, Nothing,Message{SampleList}),
                :name          => SPGaussianMeanVarianceMGNS)

@naiveVariationalRule(:node_type     => GaussianMeanVariance,
                      :outbound_type => Message{GaussianMeanVariance},
                      :inbound_types => (ProbabilityDistribution, Nothing, ProbabilityDistribution),
                      :name          => VBGaussianMeanVarianceM)

@naiveVariationalRule(:node_type     => GaussianMeanVariance,
                      :outbound_type => Message{GaussianMeanVariance},
                      :inbound_types => (Nothing, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VBGaussianMeanVarianceOut)


text:
#################
# Abstract Types
#################

export AbstractLattice, BCLattice
export isperiodic

"""
    AbstractLattice{N}

Abstract type for lattices. `N` indicated the dimension of
this lattice type.
"""
abstract type AbstractLattice{N} end

"""
    BCLattice{B, N} <: AbstractLattice{N}

Lattice with boundary condition.
"""
abstract type BCLattice{B, N} <: AbstractLattice{N} end

"""
    isperiodic(lattice) -> Bool

whether this lattice has periodic boundary.
"""
function isperiodic end

isperiodic(ltc::BCLattice{Periodic}) = true
isperiodic(ltc::BCLattice{Fixed}) = false

export LatticeIterator, SiteIterator, BondIterator
export lattice
abstract type LatticeIterator{L} end

"""
    lattice(itr)

get the lattice of `itr`.
"""
function lattice(itr::LatticeIterator) end

abstract type SiteIterator{L} <: LatticeIterator{L} end
abstract type BondIterator{L} <: LatticeIterator{L} end

####################
# Lattice Interface
####################

export shape, sites, bonds, getid

"""
    shape(lattice) -> Tuple

get the shape of a lattice
"""
function shape(lattice::AbstractLattice) end

import Base: length

"""
    length(lattice) -> Int

get the length (or the product of
size in each dimension) of the lattice
"""
function length(lattice::AbstractLattice) end

"""
    sites(lattice)

get the site iterator of the lattice.
"""
function sites(lattice::AbstractLattice) end

"""
    bonds(lattice, k)

get the ``k``th bond's iterator of the lattice
"""
function bonds(lattice::AbstractLattice) end

"""
    getid(lattice, pos)

get the given position's id on the lattice. QMTK follows
column major order for square lattice.
"""
function getid(lattice::AbstractLattice, pos) end

getid(lattice::AbstractLattice, pos::Tuple) = getid(lattice, pos...)
getid(lattice::AbstractLattice, xpos, ypos) =
    getid(lattice, xpos), getid(lattice, ypos)

getid(itr::LatticeIterator, args...) = getid(lattice(itr), args...)

################
# Fused Region
################
import IterTools
export Fusion, fuse

"""
    Fusion{T<:Tuple, L <: AbstractLattice} <: LatticeIterator

Fused Latice Iterator, like `IterTools.Chain`, but is
subtype of LatticeIterator and it only takes iterators
from the same lattice type.
"""
struct Fusion{T<:Tuple, L <: AbstractLattice} <: LatticeIterator{L}
    xss::T
end

fuse(xss::LatticeIterator{L}...) where L = Fusion{typeof(xss), L}(xss)
lattice(itr::Fusion) = lattice(first(itr.xss))

import Base: length, eltype, start, next, done

length(itr::Fusion{Tuple{}}) = 0
length(itr::Fusion) = sum(length, itr.xss)
eltype(::Type{Fusion{T}}) where T = IterTools.mapreduce_tt(eltype, typejoin, Union{}, T)

function start(itr::Fusion)
    i = 1
    xs_state = nothing
    while i <= length(itr.xss)
        xs_state = start(itr.xss[i])
        if !done(itr.xss[i], xs_state)
            break
        end
        i += 1
    end
    return i, xs_state
end

function next(itr::Fusion, state)
    i, xs_state = state
    v, xs_state = next(itr.xss[i], xs_state)
    while done(itr.xss[i], xs_state)
        i += 1
        if i > length(itr.xss)
            break
        end
        xs_state = start(itr.xss[i])
    end
    return v, (i, xs_state)
end

done(itr::Fusion, state) = state[1] > length(itr.xss)

# Implementations
include("Chain.jl")
include("Square.jl")


text:
using BSON: @load
using Flux
using GeometricFlux
using Serialization

MODELFILE="model1.bson"
DATAFILE="data.blob"

@load MODELFILE model
dat = open(DATAFILE) do io
    deserialize(io)
end

loss(gf, y) = Flux.Losses.mse(model(gf), y)
for i in 1:length(dat)
    @show loss(dat[i]["molecule"], parse(Float64, dat[i]["R2"]))
end


text:
using PrettyTables
using StatsBase

function zscore_transform(data)
    μ = mean(data)
    σ = std(data)
    z(d) = (d .- μ) ./ σ
    unz(d) = d .* σ .+ μ
    return z, unz
end

function meanlowerupper(data, PI = (0.055, 0.945))
    m = mean.(eachrow(data))
    lower = quantile.(eachrow(data), PI[1])
    upper = quantile.(eachrow(data), PI[2])
    return (mean = m,
            lower = lower,
            upper = upper,
            raw = data)
end

function estimparam(data, PI = (0.055, 0.945))
    m = mean.(eachcol(data))
    lower = quantile.(eachcol(data), PI[1])
    upper = quantile.(eachcol(data), PI[2])
    return m, lower, upper
end

function lin(a, b, c, x...)
    result = @. a + b * c
    for i in 1:2:length(x)
        @. result += x[i] * x[i+1]
    end
    return result
end

const BARS = collect("▁▂▃▄▅▆▇█")

function unicode_histogram(data, nbins = 12)
    # @show data
    f = fit(Histogram, data, nbins = nbins)  # nbins: more like a guideline than a rule, really
    # scale weights between 1 and 8 (length(BARS)) to fit the indices in BARS
    # eps is needed so indices are in the interval [0, 8) instead of [0, 8] which could
    # result in indices 0:8 which breaks things
    scaled = f.weights .* (length(BARS) / maximum(f.weights) - eps())
    indices = floor.(Int, scaled) .+ 1
    return join((BARS[i] for i in indices))
end

function precis(df::DataFrame; digits = 2, depth = Inf, alpha = 0.11)
    d = DataFrame()
    d.param = names(df)
    d.mean = mean.(eachcol(df))
    d.std = std.(eachcol(df))
    d[:, "5.5%"] = quantile.(eachcol(df), alpha/2)
    d[:, "50%"] = quantile.(eachcol(df), 0.5)
    d[:, "94.5%"] = quantile.(eachcol(df), 1 - alpha/2)
    u = unicode_histogram.(eachcol(df))
    d.histogram = unicode_histogram.(eachcol(df))

    for col in ["mean", "std", "5.5%", "50%", "94.5%"]
        d[:, col] .= round.(d[:, col], digits = digits)
    end

    pretty_table(d, nosubheader = true, vlines = [0, 1, 7])
end


text:
"""
A simple mer count struct.

MerCount is a simple struct that binds a mer value to a count of the number
of times it has been observed.

This type, (sorted) vectors of them, and some additional utility methods, can
form the basic building blocks of the higher-level mer counting functionality.
This struct can also be an eltype of kmer hashes or more involved specialized
types that store counts of Kmers.

!!! note
    The count is stored as an UInt8 because often once the count is more than
    255 we hardly care anymore.
"""
struct MerCount{M<:AbstractMer}
    mer::M
    count::UInt8
    function MerCount{M}(mer::M, count::Integer) where {M<:AbstractMer}
        return new(mer, convert(UInt8, min(typemax(UInt8), count)))
    end
end

"Shorthand for `MerCount{DNAMer{K}}`"
const DNAMerCount{K} = MerCount{DNAMer{K}}

"Shorthand for `MerCount{RNAMer{K}}`"
const RNAMerCount{K} = MerCount{RNAMer{K}}

"Get the mer from a `MerCount`."
@inline mer(x::MerCount{<:AbstractMer}) = x.mer

"Get the count from a `MerCount`."
@inline freq(x::MerCount{<:AbstractMer}) = x.count

"Get the count from a `MerCount`, and convert it to type R."
@inline freq(::Type{R}, x::MerCount{<:AbstractMer}) where {R<:Real} = convert(R, freq(x))

"""
    Base.isless(x::MerCount{M}, y::MerCount{M}) where {M<:AbstractMer}

Check if a MerCount x is less than MerCount y.

!!! note
    The comparison only considers the mer, not the count.
"""
Base.isless(x::MerCount{M}, y::MerCount{M}) where {M<:AbstractMer} = mer(x) < mer(y)

"""
    Base.:(>)(x::MerCount{M}, y::MerCount{M}) where {M<:AbstractMer}

Check if a MerCount x is greater than MerCount y.

!!! note
    The comparison only considers the mer, not the count.
"""
Base.:(>)(x::MerCount{M}, y::MerCount{M}) where {M<:AbstractMer} = mer(x) > mer(y)

function merge(x::MerCount{M}, y::MerCount{M}) where {M<:AbstractMer}
    return MerCount{M}(x.mer, UInt16(x.count) + UInt16(y.count))
end

function Base.show(io::IO, mfreq::MerCount{<:AbstractMer})
    print(io, mer(mfreq), " occurs ", freq(mfreq), " times")
end

function collect_mers(count_mode::CountMode, input::BioSequences.EveryMerIterator)
    return collect(count_mode(x) for x in input)
end

function collect_mers(::Type{M}, count_mode::CountMode, input::BioSequence) where {M<:AbstractMer}
    return collect_mers(count_mode, each(M, input))
end

collect_mers(v::Vector{M}) where {M<:AbstractMer} = v

"""
    unsafe_collapse_into_counts!(result::Vector{MerCount{M}}, mers::Vector{M}) where {M<:AbstractMer}

!!! warning
    This method is internal and marked as unsafe because it assumes that the
    `mers` input vector is already sorted.
"""
function unsafe_collapse_into_counts!(result::Vector{MerCount{M}}, mers::Vector{M}) where {M<:AbstractMer}
    wi = 1
    ri = 1
    stop = lastindex(mers) + 1
    empty!(result)
    @inbounds while ri < stop
        ci = one(UInt16)
        while (ri += 1) < stop && mers[ri] == mers[wi]
            ci = ci + one(UInt16)
        end
        push!(result, MerCount{M}(mers[wi], ci)) # TODO: See about removing this push!.
        wi = ri
    end
    return result
end

"""
    collapse_into_counts!(result::Vector{MerCount{M}}, mers::Vector{M}) where {M<:AbstractMer}

Build a vector of sorted `MerCount`s from a Vector of a mer type.

This is a basic kernel function used for any higher level and more complex
kmer counting procedures.

This is like `collapse_into_counts`, except it's first argument is a `result`
vector that is cleared and filled with the result.

!!! note
    The input vector `mers` will be sorted by this method.
"""
function collapse_into_counts!(result::Vector{MerCount{M}}, mers::Vector{M}) where {M<:AbstractMer}
    sort!(mers)
    return unsafe_collapse_into_counts!(result, mers)
end

"""
    collapse_into_counts(mers::Vector{M}) where {M<:AbstractMer}

Build a vector of sorted `MerCount`s from a Vector of a mer type.

This is a basic kernel function used for any higher level and more complex
kmer counting procedures.

!!! note
    The input vector `mers` will be sorted by this method.
"""
function collapse_into_counts!(mers::Vector{M}) where {M<:AbstractMer}
    return collapse_into_counts!(Vector{MerCount{M}}(), mers)
end

function Vector{MerCount{M}}(mode::CountMode, input::BioSequences.EveryMerIterator{M}) where {M<:AbstractMer}
    mers = collect_mers(mode, input)
    return collapse_into_counts!(mers)
end

function Vector{MerCount{M}}(mode::CountMode, input::Any) where {M<:AbstractMer}
    mers = collect_mers(M, mode, input)
    return collapse_into_counts!(mers)
end

"""
    unsafe_merge_into!(a::Vector{MerCount{M}}, b::Vector{MerCount{M}}) where {M<:AbstractMer}

Merge the `MerCount`s from vector `b` into the vector `a`.

!!! warning
    This method is marked as unsafe as it assumes both of the input vectors `a`
    and `b` are already sorted.
"""
function unsafe_merge_into!(a::Vector{MerCount{M}}, b::Vector{MerCount{M}}) where {M<:AbstractMer}
    a_i = firstindex(a)
    a_end = lastindex(a) + 1
    b_i = b_i2 = firstindex(b)
    b_end = lastindex(b) + 1

    # Merge, accumulating counts on `a`.
    @inbounds while b_i < b_end
        # Move cursor a_i forward so long as mers in a are less than the mer
        # in b[b_i]
        while a_i < a_end && a[a_i] < b[b_i]
            a_i = a_i + 1
        end
        # Check if the mers at the a_i and b_i cursors are the same. If so,
        # the count in b needs to be merged into the count in a.
        if a_i < a_end && mer(a[a_i]) == mer(b[b_i])
            # Combine entries
            a[a_i] = merge(a[a_i], b[b_i])
            a_i = a_i + 1
            b_i = b_i + 1
        end
        while b_i < b_end && (a_i == a_end || b[b_i] < a[a_i])
            b[b_i2] = b[b_i]
            b_i2 = b_i2 + 1
            b_i = b_i + 1
        end
    end
    # Shrink `b` to the size of the remaining contents.
    resize!(b, b_i2 - 1)

    # Expand `a` to allow the insertion of unique values in `b`.
    oldsize = length(a)
    resize!(a, oldsize + length(b))
    r_a = oldsize

    # Merge-sort from the bottom into `a`.
    wr_a = lastindex(a)
    rend_a = firstindex(a)
    r_b = lastindex(b)
    r_end_b = firstindex(b)
    @inbounds while wr_a >= rend_a
        if r_b >= r_end_b && (r_a < rend_a || b[r_b] > a[r_a])
            a[wr_a] = b[r_b]
            r_b = r_b - 1
        else
            a[wr_a] = a[r_a]
            r_a = r_a - 1
        end
        wr_a = wr_a - 1
    end
    empty!(b)
    return a
end

"""
    merge_into!(a::Vector{MerCount{M}}, b::Vector{MerCount{M}}) where {M<:AbstractMer}

Merge the `MerCount`s from vector `b` into the vector `a`.

!!! note
    This will sort the input vectors `a` and `b`.
"""
function merge_into!(a::Vector{MerCount{M}}, b::Vector{MerCount{M}}) where {M<:AbstractMer}
    sort!(a)
    sort!(b)
    return unsafe_merge_into!(a, b)
end

function unsafe_collapse!(freqs::Vector{MerCount{M}}) where {M<:AbstractMer}
    wi = 1
    ri = 1
    pi = 1
    stop = lastindex(freqs) + 1
    @inbounds while ri < stop
        ci = one(UInt16)
        while (ri += 1) < stop && mer(freqs[ri]) == mer(freqs[pi])
            ci = ci + one(UInt16)
        end
        freqs[wi] = MerCount{M}(mer(freqs[pi]), ci)
        pi = ri
        wi = wi + 1
    end
    resize!(freqs, wi - 1)
    return freqs
end

collapse!(freqs::Vector{MerCount{M}}) where {M<:AbstractMer} = collapse_sorted!(sort!(freqs))


text:
if VERSION >= v"0.7.0-beta.0"
    using Test
    using Dates
    using Printf
else
    using Base.Test
end
using NCDatasets

sz = (4,5)
filename = tempname()
#filename = "/tmp/test-6.nc"

# The mode "c" stands for creating a new file (clobber)
NCDatasets.Dataset(filename,"c") do ds

    # define the dimension "lon" and "lat"
    ds.dim["lon"] = sz[1]
    ds.dim["lat"] = sz[2]

    v = NCDatasets.defVar(ds,"small",Float64,("lon","lat"))
    @test_throws NCDatasets.NetCDFError v[:] = zeros(sz[1]+1,sz[2])
    #@test_throws NCDatasets.NetCDFError v[1:sz[1],1:sz[2]] = zeros(sz[1]+1,sz[2])
    @test_throws NCDatasets.NetCDFError v[sz[1]+1,1] = 1
    @test_throws NCDatasets.NetCDFError v[-1,1] = 1

    # variables
    for T in [UInt8,Int8,UInt16,Int16,UInt32,Int32,UInt64,Int64,Float32,Float64,
              Char,String]
    #for T in []
        local data
        data, scalar_data =
            if T == String
                [Char(i+60) * Char(j+60) for i = 1:sz[1], j = 1:sz[2]], "abcde"
            else
                [T(i+2*j) for i = 1:sz[1], j = 1:sz[2]], T(100)
            end

        v = NCDatasets.defVar(ds,"var-$T",T,("lon","lat"))
        v[:,:] = data
        @test v[:,:] == data[:,:]

        # issue #33
        @test Array(v) == data

        @test v[2,:] == data[2,:]

        @test v[:,3] == data[:,3]

        @test v[2,3] == data[2,3]

        # ignore extra index
        @test v[2,3,1,1] == data[2,3,1,1]

        # ignore extra index
        @test v[2:3,3,1,1] == data[2:3,3,1,1]

        @test v[[1,2,3],:] == data[[1,2,3],:]

        # write scalar,
        v.var[:,:] = scalar_data
        @test all(v.var[:,:][:] .== scalar_data)

    end
end

# quick interface
Dataset(filename,"c") do ds
    data = Int32[i+3*j for i = 1:sz[1], j = 1:sz[2]]
    defVar(ds,"temp",data,("lon","lat"), attrib = [
        "units" => "degree_Celsius",
        "long_name" => "Temperature"
    ])
    @test ds["temp"][:] == data
    @test eltype(ds["temp"].var) == Int32
    @test ds.dim["lon"] == sz[1]
    @test ds.dim["lat"] == sz[2]

    # test Union{Missing,T}
    defVar(ds,"foo",[missing,1.,2.],("dim",), fillvalue = -9999.)
    @test isequal(ds["foo"][:], [missing,1.,2.])

    defVar(ds,"scalar",123.)
    @test ds["scalar"][:] == 123.

    # load in-place
    data2 = similar(data)
    NCDatasets.load!(ds["temp"].var,data2,:,:)
    @test data2 == data

    data2 = zeros(eltype(data),sz[1],2)
    NCDatasets.load!(ds["temp"].var,data2,:,1:2)
    @test data2 == data[:,1:2]
end

# issue 23
# return type using CartesianIndex

if VERSION >= v"0.7.0-beta.0"
    filename = tempname()
    ds = Dataset(filename, "c");
    ds.dim["lon"] = 5;
    ds.dim["lat"] = 10;
    ds.dim["time"] = Inf;

    ncvar = defVar(ds, "var", Int64, ("lon", "lat", "time"));

    nt = 25;
    data = reshape(1:5*10*nt, 5, 10, nt);
    ncvar[:,:,1:nt] = data;
    close(ds);

    ds = Dataset(filename);
    start = 1;
    all(data[CartesianIndex(1, 1), start:end] .== ds["var"][CartesianIndex(1, 1), start:end])
    data11 = ds["var"][CartesianIndex(1, 1), start:end]
    close(ds)

    @test typeof(data11[1]) == Int64
    rm(filename)
end


text:
module PoolUtils

using ..CUDA

using TimerOutputs

using Printf

export MEMDEBUG

const MEMDEBUG = ccall(:jl_is_memdebug, Bool, ())

## locking

export NonReentrantLock, @spinlock, @lock

const var"@lock" = Base.var"@lock"

# a simple non-reentrant lock that errors when trying to reenter on the same task
struct NonReentrantLock <: Threads.AbstractLock
  rl::ReentrantLock
  NonReentrantLock() = new(ReentrantLock())
end

function Base.lock(nrl::NonReentrantLock)
  @assert !islocked(nrl.rl) || nrl.rl.locked_by !== current_task()
  lock(nrl.rl)
end

function Base.trylock(nrl::NonReentrantLock)
  @assert !islocked(nrl.rl) || nrl.rl.locked_by !== current_task()
  trylock(nrl.rl)
end

Base.unlock(nrl::NonReentrantLock) = unlock(nrl.rl)

# a safe way to acquire locks from finalizers, where we can't wait (which switches tasks)
macro spinlock(l, ex)
  quote
    temp = $(esc(l))
    while !trylock(temp)
      ccall(:jl_cpu_pause, Cvoid, ())
      # Temporary solution before we have gc transition support in codegen.
      ccall(:jl_gc_safepoint, Cvoid, ())
      # we can't yield here
    end
    try
      $(esc(ex))
    finally
      unlock(temp)
    end
  end
end

## block of memory

export Block, INVALID, AVAILABLE, ALLOCATED, FREED, iswhole

@enum BlockState begin
    INVALID
    AVAILABLE
    ALLOCATED
    FREED
end

mutable struct Block
    buf::Mem.DeviceBuffer # base allocation
    sz::Int               # size into it
    off::Int              # offset into it

    state::BlockState
    prev::Union{Nothing,Block}
    next::Union{Nothing,Block}

    Block(buf, sz; off=0, state=INVALID, prev=nothing, next=nothing) =
        new(buf, sz, off, state, prev, next)
end

Base.sizeof(block::Block) = block.sz
Base.pointer(block::Block) = pointer(block.buf) + block.off

iswhole(block::Block) = block.prev === nothing && block.next === nothing

function Base.show(io::IO, block::Block)
    fields = [@sprintf("%p", Int(pointer(block)))]
    push!(fields, Base.format_bytes(sizeof(block)))
    push!(fields, "$(block.state)")
    block.off != 0 && push!(fields, "offset=$(block.off)")
    block.prev !== nothing && push!(fields, "prev=Block(offset=$(block.prev.off))")
    block.next !== nothing && push!(fields, "next=Block(offset=$(block.next.off))")

    print(io, "Block(", join(fields, ", "), ")")
end

## timing

export @pool_timeit

const to = TimerOutput()

macro pool_timeit(args...)
    TimerOutputs.timer_expr(CUDA, true, :($PoolUtils.to), args...)
end

end


text:
int main() {
    if(true);
    else return 0;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
# Autogenerated wrapper script for Blosc_jll for x86_64-linux-gnu-cxx03
export libblosc

using Zlib_jll
using Zstd_jll
using Lz4_jll
JLLWrappers.@generate_wrapper_header("Blosc")
JLLWrappers.@declare_library_product(libblosc, "libblosc.so.1")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, Zstd_jll, Lz4_jll)
    JLLWrappers.@init_library_product(
        libblosc,
        "lib/libblosc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using StatsPlots, Revise
includet("multiple_good_market.jl")

τ_truth = 0.172 #0global_estimate()
simulation = mc_simulation(10000, 1000)
println("Mean: ", mean(simulation; dims=1))
println("SD: ", std(simulation; dims = 1))
density(simulation[:, 1:2], label=["Estimated Marginal Policy Effect" "Estimated Direct Effect"], linestyle = [:solid :dash])
vline!([τ_truth], label= "True Marginal Policy Effect", linestyle = :dot, linecolor=:black, linewidth=2.0)
savefig("simulation.pdf")

max_stars_repo_name: evanmunro/market-equilibrium

max_issues_repo_name: evanmunro/market-equilibrium

max_forks_repo_name: evanmunro/market-equilibrium


text:
using CairoMakie, Random
Random.seed!(2)
t = 0.3:0.3:3π
my_markers = [:circle, :rect, :utriangle, :dtriangle, :diamond, :pentagon,
    :cross, :xcross]
with_theme(theme_dark()) do
    fig, ax, = stem(t, 1.5exp.(-t/5).*cos.(t); color = 1:length(t),
        colormap = :linear_wyor_100_45_c55_n256, stemcolor = 1:length(t),
        stemcolormap = :linear_wcmr_100_45_c42_n256,
        figure = (; resolution = (600,400)))
    stem!(t .+ 0.15, -cos.(t) ./ t .+ 0.25; color = :transparent,
        stemwidth = 0.5, marker = :rect, markersize = 10, strokewidth = 1,
        strokecolor = :white)
    stem!(1:8, 1.5randn(8); marker = my_markers,
        color = tuple.(to_colormap(:mk_12, 8), 0.65), stemlinestyle = :dash,
        stemcolor = to_colormap(:mk_12, 8), markersize = 15*rand(8) .+ 10,
        strokewidth = 1.5, strokecolor = to_colormap(:mk_12, 8))
    hidedecorations!(ax; grid = false)
    save("stem.png", fig, px_per_unit = 2)
    display(fig)
end


text:
"Supertype of lazy grids, defined in terms of other grids."
abstract type LazyGrid{T,N} <: AbstractGrid{T,N}
end

"Grid defined in terms of a single `supergrid`."
abstract type SimpleLazyGrid{T,N} <: AbstractGrid{T,N}
end

supergrid(g::SimpleLazyGrid) = g.grid
supergrid(g::SimpleLazyGrid, I...) = supergrid(g)[I...]

"Grid defined in terms of multiple grids."
abstract type CompositeGrid{T,N} <: LazyGrid{T,N}
end

components(grid::CompositeGrid) = grid.grids


text:
using DayCounts, Dates
using Test

# based on tests from OpenGamma Strata
# https://github.com/OpenGamma/Strata/blob/efaa0be0d08dc61c23692e7b86df101ea0fb7223/modules/basics/src/test/java/com/opengamma/strata/basics/date/DayCountTest.java
@testset "Actual365Fixed" begin
    dc = DayCounts.Actual365Fixed()
    @test yearfrac(Date(2011,12,28), Date(2012, 2,28), dc) == (4 + 58) / 365
    @test yearfrac(Date(2011,12,28), Date(2012, 2,29), dc) == (4 + 59) / 365
    @test yearfrac(Date(2011,12,28), Date(2012, 3, 1), dc) == (4 + 60) / 365
    @test yearfrac(Date(2011,12,28), Date(2016, 2,28), dc) == (4 + 366 + 365 + 365 + 365 + 58) / 365
    @test yearfrac(Date(2011,12,28), Date(2016, 2,29), dc) == (4 + 366 + 365 + 365 + 365 + 59) / 365
    @test yearfrac(Date(2011,12,28), Date(2016, 3, 1), dc) == (4 + 366 + 365 + 365 + 365 + 60) / 365
    @test yearfrac(Date(2012, 2,28), Date(2012, 3,28), dc) == 29 / 365
    @test yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc) == 28 / 365
    @test yearfrac(Date(2012, 3, 1), Date(2012, 3,28), dc) == 27 / 365

    # zeros
    @test yearfrac(Date(2011,12,28), Date(2011,12,28), dc) == 0
    @test yearfrac(Date(2011,12,31), Date(2011,12,31), dc) == 0
    @test yearfrac(Date(2012, 2,28), Date(2012, 2,28), dc) == 0
    @test yearfrac(Date(2012, 2,29), Date(2012, 2,29), dc) == 0

    # reflection
    @test yearfrac(Date(2012, 2,28), Date(2011,12,28), dc) == -yearfrac(Date(2011,12,28), Date(2012, 2,28), dc)
    @test yearfrac(Date(2012, 3,28), Date(2012, 2,29), dc) == -yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc)
end

@testset "Actual360" begin
    dc = DayCounts.Actual360()
    @test yearfrac(Date(2011,12,28), Date(2012, 2,28), dc) == (4 + 58) / 360
    @test yearfrac(Date(2011,12,28), Date(2012, 2,29), dc) == (4 + 59) / 360
    @test yearfrac(Date(2011,12,28), Date(2012, 3, 1), dc) == (4 + 60) / 360
    @test yearfrac(Date(2011,12,28), Date(2016, 2,28), dc) == (4 + 366 + 365 + 365 + 365 + 58) / 360
    @test yearfrac(Date(2011,12,28), Date(2016, 2,29), dc) == (4 + 366 + 365 + 365 + 365 + 59) / 360
    @test yearfrac(Date(2011,12,28), Date(2016, 3, 1), dc) == (4 + 366 + 365 + 365 + 365 + 60) / 360
    @test yearfrac(Date(2012, 2,28), Date(2012, 3,28), dc) == 29 / 360
    @test yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc) == 28 / 360
    @test yearfrac(Date(2012, 3, 1), Date(2012, 3,28), dc) == 27 / 360

    # zeros
    @test yearfrac(Date(2011,12,28), Date(2011,12,28), dc) == 0
    @test yearfrac(Date(2011,12,31), Date(2011,12,31), dc) == 0
    @test yearfrac(Date(2012, 2,28), Date(2012, 2,28), dc) == 0
    @test yearfrac(Date(2012, 2,29), Date(2012, 2,29), dc) == 0

    # reflection
    @test yearfrac(Date(2012, 2,28), Date(2011,12,28), dc) == -yearfrac(Date(2011,12,28), Date(2012, 2,28), dc)
    @test yearfrac(Date(2012, 3,28), Date(2012, 2,29), dc) == -yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc)
end

@testset "ActualActualISDA" begin
    dc = DayCounts.ActualActualISDA()
    @test yearfrac(Date(2011,12,28), Date(2012, 2,28), dc) == 4 / 365 + 58 / 366
    @test yearfrac(Date(2011,12,28), Date(2012, 2,29), dc) == 4 / 365 + 59 / 366
    @test yearfrac(Date(2011,12,28), Date(2012, 3, 1), dc) == 4 / 365 + 60 / 366
    @test yearfrac(Date(2011,12,28), Date(2016, 2,28), dc) == 4 / 365 + 58 / 366 + 4
    @test yearfrac(Date(2011,12,28), Date(2016, 2,29), dc) == 4 / 365 + 59 / 366 + 4
    @test yearfrac(Date(2011,12,28), Date(2016, 3, 1), dc) == 4 / 365 + 60 / 366 + 4
    @test yearfrac(Date(2012, 2,28), Date(2012, 3,28), dc) == 29 / 366
    @test yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc) == 28 / 366
    @test yearfrac(Date(2012, 3, 1), Date(2012, 3,28), dc) == 27 / 366

    # zeros
    @test yearfrac(Date(2011,12,28), Date(2011,12,28), dc) == 0
    @test yearfrac(Date(2011,12,31), Date(2011,12,31), dc) == 0
    @test yearfrac(Date(2012, 2,28), Date(2012, 2,28), dc) == 0
    @test yearfrac(Date(2012, 2,29), Date(2012, 2,29), dc) == 0

    # reflection
    @test yearfrac(Date(2012, 2,28), Date(2011,12,28), dc) == -yearfrac(Date(2011,12,28), Date(2012, 2,28), dc)
    @test yearfrac(Date(2012, 3,28), Date(2012, 2,29), dc) == -yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc)
end

@testset "Thirty360" begin
    dc = DayCounts.Thirty360()
    # usual rule: ((d2-d1) + (m2-m1)*30 + (y2-y1)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2012, 2,28), dc) == ((28-28) + (2-12)*30 + (2012-2011)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2012, 2,29), dc) == ((29-28) + (2-12)*30 + (2012-2011)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2012, 3, 1), dc) == ((1-28) + (3-12)*30 + (2012-2011)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2016, 2,28), dc) == ((28-28) + (2-12)*30 + (2016-2011)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2016, 2,29), dc) == ((29-28) + (2-12)*30 + (2016-2011)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2016, 3, 1), dc) == ((1-28) + (3-12)*30 + (2016-2011)*360)/360
    @test yearfrac(Date(2012, 2,28), Date(2012, 3,28), dc) == ((28-28) + (3-2)*30 + (2012-2012)*360)/360
    @test yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc) == ((28-29) + (3-2)*30 + (2012-2012)*360)/360
    @test yearfrac(Date(2012, 3, 1), Date(2012, 3,28), dc) == ((28-1) + (3-3)*30 + (2012-2012)*360)/360

    @test yearfrac(Date(2012, 5,29), Date(2013, 8,29), dc) == ((29-29) + (8-5)*30 + (2013-2012)*360)/360
    @test yearfrac(Date(2012, 5,29), Date(2013, 8,30), dc) == ((30-29) + (8-5)*30 + (2013-2012)*360)/360
    @test yearfrac(Date(2012, 5,29), Date(2013, 8,31), dc) == ((31-29) + (8-5)*30 + (2013-2012)*360)/360
    @test yearfrac(Date(2012, 5,30), Date(2013, 8,29), dc) == ((29-30) + (8-5)*30 + (2013-2012)*360)/360
    @test yearfrac(Date(2012, 5,30), Date(2013, 8,30), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360
    @test yearfrac(Date(2012, 5,30), Date(2013, 8,31), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360 # exception
    @test yearfrac(Date(2012, 5,31), Date(2013, 8,29), dc) == ((29-30) + (8-5)*30 + (2013-2012)*360)/360 # exception
    @test yearfrac(Date(2012, 5,31), Date(2013, 8,30), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360 # exception
    @test yearfrac(Date(2012, 5,31), Date(2013, 8,31), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360 # exception

    # zeros
    @test yearfrac(Date(2011,12,28), Date(2011,12,28), dc) == 0
    @test yearfrac(Date(2011,12,31), Date(2011,12,31), dc) == 0
    @test yearfrac(Date(2012, 2,28), Date(2012, 2,28), dc) == 0
    @test yearfrac(Date(2012, 2,29), Date(2012, 2,29), dc) == 0

    # reflection
    @test yearfrac(Date(2012, 2,28), Date(2011,12,28), dc) == -yearfrac(Date(2011,12,28), Date(2012, 2,28), dc)
    @test yearfrac(Date(2012, 3,28), Date(2012, 2,29), dc) == -yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc)
end

@testset "ThirtyE360" begin
    dc = DayCounts.ThirtyE360()
    # usual rule: ((d2-d1) + (m2-m1)*30 + (y2-y1)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2012, 2,28), dc) == ((28-28) + (2-12)*30 + (2012-2011)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2012, 2,29), dc) == ((29-28) + (2-12)*30 + (2012-2011)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2012, 3, 1), dc) == ((1-28) + (3-12)*30 + (2012-2011)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2016, 2,28), dc) == ((28-28) + (2-12)*30 + (2016-2011)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2016, 2,29), dc) == ((29-28) + (2-12)*30 + (2016-2011)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2016, 3, 1), dc) == ((1-28) + (3-12)*30 + (2016-2011)*360)/360
    @test yearfrac(Date(2012, 2,28), Date(2012, 3,28), dc) == ((28-28) + (3-2)*30 + (2012-2012)*360)/360
    @test yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc) == ((28-29) + (3-2)*30 + (2012-2012)*360)/360
    @test yearfrac(Date(2012, 3, 1), Date(2012, 3,28), dc) == ((28-1) + (3-3)*30 + (2012-2012)*360)/360

    @test yearfrac(Date(2012, 5,29), Date(2013, 8,29), dc) == ((29-29) + (8-5)*30 + (2013-2012)*360)/360
    @test yearfrac(Date(2012, 5,29), Date(2013, 8,30), dc) == ((30-29) + (8-5)*30 + (2013-2012)*360)/360
    @test yearfrac(Date(2012, 5,29), Date(2013, 8,31), dc) == ((30-29) + (8-5)*30 + (2013-2012)*360)/360 # exception
    @test yearfrac(Date(2012, 5,30), Date(2013, 8,29), dc) == ((29-30) + (8-5)*30 + (2013-2012)*360)/360
    @test yearfrac(Date(2012, 5,30), Date(2013, 8,30), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360
    @test yearfrac(Date(2012, 5,30), Date(2013, 8,31), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360 # exception
    @test yearfrac(Date(2012, 5,31), Date(2013, 8,29), dc) == ((29-30) + (8-5)*30 + (2013-2012)*360)/360 # exception
    @test yearfrac(Date(2012, 5,31), Date(2013, 8,30), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360 # exception
    @test yearfrac(Date(2012, 5,31), Date(2013, 8,31), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360 # exception

    # zeros
    @test yearfrac(Date(2011,12,28), Date(2011,12,28), dc) == 0
    @test yearfrac(Date(2011,12,31), Date(2011,12,31), dc) == 0
    @test yearfrac(Date(2012, 2,28), Date(2012, 2,28), dc) == 0
    @test yearfrac(Date(2012, 2,29), Date(2012, 2,29), dc) == 0

    # reflection
    @test yearfrac(Date(2012, 2,28), Date(2011,12,28), dc) == -yearfrac(Date(2011,12,28), Date(2012, 2,28), dc)
    @test yearfrac(Date(2012, 3,28), Date(2012, 2,29), dc) == -yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc)
end

@testset "ThirtyE360ISDA" begin
    @testset "maturity not end of Feb" begin
        dc = DayCounts.ThirtyE360ISDA(Date(2012, 8, 31))
        @test yearfrac(Date(2011,12,28), Date(2012, 2,28), dc) == ((28-28) + (2-12)*30 + (2012-2011)*360)/360
        @test yearfrac(Date(2011,12,28), Date(2012, 2,29), dc) == ((30-28) + (2-12)*30 + (2012-2011)*360)/360 # exception
        @test yearfrac(Date(2011,12,28), Date(2012, 3, 1), dc) == ((1-28) + (3-12)*30 + (2012-2011)*360)/360
        @test yearfrac(Date(2011,12,28), Date(2016, 2,28), dc) == ((28-28) + (2-12)*30 + (2016-2011)*360)/360
        @test yearfrac(Date(2011,12,28), Date(2016, 2,29), dc) == ((30-28) + (2-12)*30 + (2016-2011)*360)/360 # exception
        @test yearfrac(Date(2011,12,28), Date(2016, 3, 1), dc) == ((1-28) + (3-12)*30 + (2016-2011)*360)/360
        @test yearfrac(Date(2012, 2,28), Date(2012, 3,28), dc) == ((28-28) + (3-2)*30 + (2012-2012)*360)/360
        @test yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc) == ((28-30) + (3-2)*30 + (2012-2012)*360)/360 # exception
        @test yearfrac(Date(2012, 3, 1), Date(2012, 3,28), dc) == ((28-1) + (3-3)*30 + (2012-2012)*360)/360

        @test yearfrac(Date(2012, 5,29), Date(2013, 8,29), dc) == ((29-29) + (8-5)*30 + (2013-2012)*360)/360
        @test yearfrac(Date(2012, 5,29), Date(2013, 8,30), dc) == ((30-29) + (8-5)*30 + (2013-2012)*360)/360
        @test yearfrac(Date(2012, 5,29), Date(2013, 8,31), dc) == ((30-29) + (8-5)*30 + (2013-2012)*360)/360 # exception
        @test yearfrac(Date(2012, 5,30), Date(2013, 8,29), dc) == ((29-30) + (8-5)*30 + (2013-2012)*360)/360
        @test yearfrac(Date(2012, 5,30), Date(2013, 8,30), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360
        @test yearfrac(Date(2012, 5,30), Date(2013, 8,31), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360 # exception
        @test yearfrac(Date(2012, 5,31), Date(2013, 8,29), dc) == ((29-30) + (8-5)*30 + (2013-2012)*360)/360 # exception
        @test yearfrac(Date(2012, 5,31), Date(2013, 8,30), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360 # exception
        @test yearfrac(Date(2012, 5,31), Date(2013, 8,31), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360 # exception

        # zeros
        @test yearfrac(Date(2011,12,28), Date(2011,12,28), dc) == 0
        @test yearfrac(Date(2011,12,31), Date(2011,12,31), dc) == 0
        @test yearfrac(Date(2012, 2,28), Date(2012, 2,28), dc) == 0
        @test yearfrac(Date(2012, 2,29), Date(2012, 2,29), dc) == 0

        # reflection
        @test yearfrac(Date(2012, 2,28), Date(2011,12,28), dc) == -yearfrac(Date(2011,12,28), Date(2012, 2,28), dc)
        @test yearfrac(Date(2012, 3,28), Date(2012, 2,29), dc) == -yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc)
    end

    @testset "maturity end of Feb" begin
        dc = DayCounts.ThirtyE360ISDA(Date(2012, 2,29)) # leap year, end of month
        @test yearfrac(Date(2011,12,28), Date(2012, 2,28), dc) == ((28-28) + (2-12)*30 + (2012-2011)*360)/360
        @test yearfrac(Date(2011,12,28), Date(2012, 2,29), dc) == ((29-28) + (2-12)*30 + (2012-2011)*360)/360
        @test yearfrac(Date(2011,12,28), Date(2012, 3, 1), dc) == ((1-28) + (3-12)*30 + (2012-2011)*360)/360
        @test yearfrac(Date(2011,12,28), Date(2016, 2,28), dc) == ((28-28) + (2-12)*30 + (2016-2011)*360)/360
        @test yearfrac(Date(2011,12,28), Date(2016, 2,29), dc) == ((30-28) + (2-12)*30 + (2016-2011)*360)/360 # exception
        @test yearfrac(Date(2011,12,28), Date(2016, 3, 1), dc) == ((1-28) + (3-12)*30 + (2016-2011)*360)/360
        @test yearfrac(Date(2012, 2,28), Date(2012, 3,28), dc) == ((28-28) + (3-2)*30 + (2012-2012)*360)/360
        @test yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc) == ((28-30) + (3-2)*30 + (2012-2012)*360)/360 # exception
        @test yearfrac(Date(2012, 3, 1), Date(2012, 3,28), dc) == ((28-1) + (3-3)*30 + (2012-2012)*360)/360

        @test yearfrac(Date(2012, 5,29), Date(2013, 8,29), dc) == ((29-29) + (8-5)*30 + (2013-2012)*360)/360
        @test yearfrac(Date(2012, 5,29), Date(2013, 8,30), dc) == ((30-29) + (8-5)*30 + (2013-2012)*360)/360
        @test yearfrac(Date(2012, 5,29), Date(2013, 8,31), dc) == ((30-29) + (8-5)*30 + (2013-2012)*360)/360 # exception
        @test yearfrac(Date(2012, 5,30), Date(2013, 8,29), dc) == ((29-30) + (8-5)*30 + (2013-2012)*360)/360
        @test yearfrac(Date(2012, 5,30), Date(2013, 8,30), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360
        @test yearfrac(Date(2012, 5,30), Date(2013, 8,31), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360 # exception
        @test yearfrac(Date(2012, 5,31), Date(2013, 8,29), dc) == ((29-30) + (8-5)*30 + (2013-2012)*360)/360 # exception
        @test yearfrac(Date(2012, 5,31), Date(2013, 8,30), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360 # exception
        @test yearfrac(Date(2012, 5,31), Date(2013, 8,31), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360 # exception

        # zeros
        @test yearfrac(Date(2011,12,28), Date(2011,12,28), dc) == 0
        @test yearfrac(Date(2011,12,31), Date(2011,12,31), dc) == 0
        @test yearfrac(Date(2012, 2,28), Date(2012, 2,28), dc) == 0
        @test yearfrac(Date(2012, 2,29), Date(2012, 2,29), dc) == 0

        # reflection
        @test yearfrac(Date(2012, 2,28), Date(2011,12,28), dc) == -yearfrac(Date(2011,12,28), Date(2012, 2,28), dc)
        @test yearfrac(Date(2012, 3,28), Date(2012, 2,29), dc) == -yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc)
    end
end

@testset "ActualActualICMA" begin
    dc = DayCounts.ActualActualICMA(Date(2011,7,1):Month(6):Date(2012,7,1))

    # zeros
    @test yearfrac(Date(2011,12,28), Date(2011,12,28), dc) == 0
    @test yearfrac(Date(2011,12,31), Date(2011,12,31), dc) == 0
    @test yearfrac(Date(2012, 2,28), Date(2012, 2,28), dc) == 0
    @test yearfrac(Date(2012, 2,29), Date(2012, 2,29), dc) == 0

    # reflection
    @test yearfrac(Date(2012, 2,28), Date(2011,12,28), dc) == -yearfrac(Date(2011,12,28), Date(2012, 2,28), dc)
    @test yearfrac(Date(2012, 3,28), Date(2012, 2,29), dc) == -yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc)
end

@testset "ActuaActual ISDA Memo" begin
    # these test cases are all derived from
    # The Actual/Actual Day Count Fraction, ISDA
    # Paper for use with the ISDA Market Conventions Survey - 3rd June 1999
    # https://www.isda.org/a/pIJEE/The-Actual-Actual-Day-Count-Fraction-1999.pdf

    @testset "Semi-annual payment" begin
        dc_isda = DayCounts.ActualActualISDA()
        dc_icma = DayCounts.ActualActualICMA(Date(2003,11,1):Month(6):Date(2004,5,1))

        @test yearfrac(Date(2003,11,1), Date(2004,5,1), dc_isda) == (61/365 + 121/366)
        @test yearfrac(Date(2003,11,1), Date(2004,5,1), dc_icma) == 182/(2*182)
    end

    @testset "short first" begin
        dc_isda = DayCounts.ActualActualISDA()
        dc_icma = DayCounts.ActualActualICMA(Date(1998,7,1):Month(12):Date(2000,7,1))

        @test yearfrac(Date(1999,2,1), Date(1999,7,1), dc_isda) == 150/365
        @test yearfrac(Date(1999,2,1), Date(1999,7,1), dc_icma) ≈ 150/365 atol=eps()

        @test yearfrac(Date(1999,7,1), Date(2000,7,1), dc_isda) == 184/365 + 182/366
        @test yearfrac(Date(1999,7,1), Date(2000,7,1), dc_icma) == 366/366

        @test yearfrac(Date(1999,2,1), Date(2000,7,1), dc_isda) == (150 + 184)/365 + 182/366
        @test yearfrac(Date(1999,2,1), Date(2000,7,1), dc_icma) ≈ 150/365 + 366/366 atol=eps()
    end

    @testset "long first" begin
        dc_isda = DayCounts.ActualActualISDA()
        dc_icma = DayCounts.ActualActualICMA(Date(2002,7,15):Month(6):Date(2004,1,15))

        @test yearfrac(Date(2002,8,15), Date(2003,7,15), dc_isda) == 334/365
        @test yearfrac(Date(2002,8,15), Date(2003,7,15), dc_icma) == 181/(181*2) + 153/(184*2)

        @test yearfrac(Date(2003,7,15),Date(2004,1,15), dc_isda) == 170/365 + 14/366
        @test yearfrac(Date(2003,7,15), Date(2004,1,15), dc_icma) == 184/(184*2)

        @test yearfrac(Date(2002,8,15), Date(2004,1,15), dc_isda) ≈ (334+170)/365 + 14/366 rtol=eps()
        @test yearfrac(Date(2002,8,15), Date(2004,1,15), dc_icma) == 181/(181*2) + (153+184)/(184*2)
    end

    @testset "short final" begin
        dc_isda = DayCounts.ActualActualISDA()
        dc_icma = DayCounts.ActualActualICMA(Date(1999,7,30):Month(6):Date(2000,7,30))

        @test yearfrac(Date(1999,7,30), Date(2000,1,30), dc_isda) == 155/365 + 29/366
        @test yearfrac(Date(1999,7,30), Date(2000,1,30), dc_icma) == 184/(184*2)

        @test yearfrac(Date(2000,1,30), Date(2000,6,30), dc_isda) == 152/366
        @test yearfrac(Date(2000,1,30), Date(2000,6,30), dc_icma) == 152/(182*2)

        @test yearfrac(Date(1999,7,30), Date(2000,6,30), dc_isda) == 155/365 + (29+152)/366
        @test yearfrac(Date(1999,7,30), Date(2000,6,30), dc_icma) == 184/(184*2) + 152/(182*2)
    end

    @testset "long final" begin
        dc_isda = DayCounts.ActualActualISDA()
        # TODO: need a way to express "last day of month" schedules.
        dc_icma = DayCounts.ActualActualICMA(Date(1999,5,31):Month(3):Date(2000,5,31))

        @test yearfrac(Date(1999,11,30), Date(2000,4,30), dc_isda) == 32/365 + 120/366
        @test yearfrac(Date(1999,11,30), Date(2000,4,30), dc_icma) == 91/(91*4) + 61/(92*4)
    end
end

@testset "Thirty360Excel" begin
    dc = DayCounts.Thirty360Excel()
    # usual rule: ((d2-d1) + (m2-m1)*30 + (y2-y1)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2012, 2,28), dc) == ((28-28) + (2-12)*30 + (2012-2011)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2012, 2,29), dc) == ((29-28) + (2-12)*30 + (2012-2011)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2012, 3, 1), dc) == ((1-28) + (3-12)*30 + (2012-2011)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2016, 2,28), dc) == ((28-28) + (2-12)*30 + (2016-2011)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2016, 2,29), dc) == ((29-28) + (2-12)*30 + (2016-2011)*360)/360
    @test yearfrac(Date(2011,12,28), Date(2016, 3, 1), dc) == ((1-28) + (3-12)*30 + (2016-2011)*360)/360
    @test yearfrac(Date(2012, 2,28), Date(2012, 3,28), dc) == ((28-28) + (3-2)*30 + (2012-2012)*360)/360
    @test yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc) == ((28-30) + (3-2)*30 + (2012-2012)*360)/360 # exception
    @test yearfrac(Date(2012, 3, 1), Date(2012, 3,28), dc) == ((28-1) + (3-3)*30 + (2012-2012)*360)/360

    @test yearfrac(Date(2012, 5,29), Date(2013, 8,29), dc) == ((29-29) + (8-5)*30 + (2013-2012)*360)/360
    @test yearfrac(Date(2012, 5,29), Date(2013, 8,30), dc) == ((30-29) + (8-5)*30 + (2013-2012)*360)/360
    @test yearfrac(Date(2012, 5,29), Date(2013, 8,31), dc) == ((31-29) + (8-5)*30 + (2013-2012)*360)/360
    @test yearfrac(Date(2012, 5,30), Date(2013, 8,29), dc) == ((29-30) + (8-5)*30 + (2013-2012)*360)/360
    @test yearfrac(Date(2012, 5,30), Date(2013, 8,30), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360
    @test yearfrac(Date(2012, 5,30), Date(2013, 8,31), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360 # exception
    @test yearfrac(Date(2012, 5,31), Date(2013, 8,29), dc) == ((29-30) + (8-5)*30 + (2013-2012)*360)/360 # exception
    @test yearfrac(Date(2012, 5,31), Date(2013, 8,30), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360 # exception
    @test yearfrac(Date(2012, 5,31), Date(2013, 8,31), dc) == ((30-30) + (8-5)*30 + (2013-2012)*360)/360 # exception

    # zeros
    @test yearfrac(Date(2011,12,28), Date(2011,12,28), dc) == 0
    @test yearfrac(Date(2011,12,31), Date(2011,12,31), dc) == 0
    @test yearfrac(Date(2012, 2,28), Date(2012, 2,28), dc) == 0
    @test yearfrac(Date(2012, 2,29), Date(2012, 2,29), dc) == 0

    # reflection
    @test yearfrac(Date(2012, 2,28), Date(2011,12,28), dc) == -yearfrac(Date(2011,12,28), Date(2012, 2,28), dc)
    @test yearfrac(Date(2012, 3,28), Date(2012, 2,29), dc) == -yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc)
end

@testset "ActualActualExcel" begin
    dc = DayCounts.ActualActualExcel()
    # calculated from Excel sheet
    @test yearfrac(Date(2011,12,28), Date(2012, 2,28), dc) == 62/365
    @test yearfrac(Date(2011,12,28), Date(2012, 2,29), dc) == 63/366
    @test yearfrac(Date(2011,12,28), Date(2012, 3, 1), dc) == 64/366
    @test yearfrac(Date(2011,12,28), Date(2016, 2,28), dc) == 1523/365.3333333333333
    @test yearfrac(Date(2011,12,28), Date(2016, 2,29), dc) == 1524/365.3333333333333
    @test yearfrac(Date(2011,12,28), Date(2016, 3, 1), dc) == 1525/365.3333333333333
    @test yearfrac(Date(2012, 2,28), Date(2012, 3,28), dc) == 29/366
    @test yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc) == 28/366
    @test yearfrac(Date(2012, 3, 1), Date(2012, 3,28), dc) == 27/366

    @test yearfrac(Date(2012, 2,28), Date(2013, 2,28), dc) == 366/366
    @test yearfrac(Date(2012, 2,29), Date(2013, 2,28), dc) == 365/366
    @test yearfrac(Date(2012, 3, 1), Date(2013, 2,28), dc) == 364/365
    @test yearfrac(Date(2012, 2,28), Date(2013, 3, 1), dc) == 367/365.5
    @test yearfrac(Date(2012, 2,29), Date(2013, 3, 1), dc) == 366/365.5
    @test yearfrac(Date(2012, 3, 1), Date(2013, 3, 1), dc) == 365/365

    @test yearfrac(Date(2011, 2,28), Date(2012, 2,28), dc) == 365/365
    @test yearfrac(Date(2011, 3, 1), Date(2012, 2,28), dc) == 364/365
    @test yearfrac(Date(2011, 2,28), Date(2012, 2,29), dc) == 366/365.5
    @test yearfrac(Date(2011, 3, 1), Date(2012, 2,29), dc) == 365/366
    @test yearfrac(Date(2011, 2,28), Date(2012, 3, 1), dc) == 367/365.5
    @test yearfrac(Date(2011, 3, 1), Date(2012, 3, 1), dc) == 366/366

    # zeros
    @test yearfrac(Date(2011,12,28), Date(2011,12,28), dc) == 0
    @test yearfrac(Date(2011,12,31), Date(2011,12,31), dc) == 0
    @test yearfrac(Date(2012, 2,28), Date(2012, 2,28), dc) == 0
    @test yearfrac(Date(2012, 2,29), Date(2012, 2,29), dc) == 0

    # reflection
    @test yearfrac(Date(2012, 2,28), Date(2011,12,28), dc) == -yearfrac(Date(2011,12,28), Date(2012, 2,28), dc)
    @test yearfrac(Date(2012, 3,28), Date(2012, 2,29), dc) == -yearfrac(Date(2012, 2,29), Date(2012, 3,28), dc)
end

include("excel.jl")


text:
using ParameterizedFunctions, OrdinaryDiffEq, Test

f = @ode_def begin
    dGut    = -Ka*Gut
    dCent   =  Ka*Gut - CL*(Cent/Vc) -Q*(Cent/Vc) + Q*(Periph/Vp)
    dPeriph =  Q*(Cent/Vc)  - Q*(Periph/Vp)
    dResp   =  Kin*(1-(IMAX*(Cent/Vc)^γ/(IC50^γ+(Cent/Vc)^γ)))  - Kout*Resp
end Ka CL Vc Q Vp Kin Kout IC50 IMAX γ

p = (Ka = 2.769674683505292, CL = 0.8315548845614537, Vc = 45.87339297906569,
      Q = 1.4544520084387496, Vp = 9.460961128900388, Kin = 6.360418843010778,
      Kout = 3.524234094210465, IC50 = 0.4483840673123028, IMAX = 1.0, γ = 2.0,)

u0 = [0.35793622925128044,
      0.7123008259214687,
      0.10685666303615937,
      0.00]
prob = ODEProblem(f, u0, (0,26.), p)
sol = solve(prob, AutoTsit5(Rodas5(), dtfac=2))

@test sol([21.2])[1][4] ≈ 1.803 atol = 1e-2

#=
using Plots
plot(sol, vars=4, lab="dtfac=1")
sol = solve(prob, AutoTsit5(Rodas5(), dtfac=2))
plt1 = plot!(sol, vars=4, lab="dtfac=2")
plt2 = plot(sol.t, sol.alg_choice, lab="alg_choice")
plot(plt1, plt2, layout=(2,1))
plot!(;title="AutoTsit5(Rodas5())",dpi=200)
=#


text:
# Julia code
# Copyright (c) 2022 the SLICOTMath.jl developers
# Portions extracted from SLICOT-Reference distribution:
# Copyright (c) 2002-2020 NICONET e.V.

function run_mb04yd(datfile, io=stdout)
    ZERO = 0.0e0
    NIN = 5
    NOUT = 6
    MMAX = 20
    NMAX = 20
    MNMIN = min( MMAX, NMAX )
    LDU = MMAX
    LDV = NMAX
    LDWORK = 6*MNMIN - 5
    Q = Array{Float64,1}(undef, MNMIN)
    E = Array{Float64,1}(undef, MNMIN-1)
    U = Array{Float64,2}(undef, LDU,MNMIN)
    V = Array{Float64,2}(undef, LDV,MNMIN)
    INUL = Array{BlasBool,1}(undef, MNMIN)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    M = parse(BlasInt, vs[1])
    N = parse(BlasInt, vs[2])
    THETA = parse(Float64, replace(vs[3],'D'=>'E'))
    RANK = parse(BlasInt, vs[4])
    TOL = parse(Float64, replace(vs[5],'D'=>'E'))
    RELTOL = parse(Float64, replace(vs[6],'D'=>'E'))
    JOBU =  vs[7][1]
    JOBV =  vs[8][1]
    MINMN = min( M, N )
if ( M<0 || M>MMAX )
elseif ( N<0 || N>NMAX )
elseif ( RANK>MINMN )
elseif ( RANK<0 && THETA<ZERO )
else

    vs = String[]
    _isz = MINMN
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    Q[1:_isz] .= parsex.(Float64, vs)

    vs = String[]
    _isz = MINMN-1
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    E[1:_isz] .= parsex.(Float64, vs)

    RANK1 = RANK
    LJOBUU = LSAME( JOBU, 'U' )
    LJOBVU = LSAME( JOBV, 'U' )
    # CHECKME: partly translated condition
    if  LJOBUU

    vs = String[]
    _isz,_jsz = (M,MINMN)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       U[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    end
    # CHECKME: partly translated condition
    if  LJOBVU

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

    INUL .= false
    # CHECKME: partly translated condition
    if  LJOBUU||LJOBVU

    vs = String[]
    _isz = MINMN
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    INUL[1:_isz] .= parsex.(Bool, vs)

    end
# interp call 1

    RANK, THETA, INFO, IWARN = SLICOT.mb04yd!(JOBU, JOBV, M, N, RANK, THETA, Q, E, U, V, INUL, TOL, RELTOL, LDWORK)
    println(io, "RANK = $RANK")
    println(io, "THETA = $THETA")
    @test INFO == 0
    INFO == 0 || return
    println(io, "IWARN = $IWARN")

if ( INFO!=0 )
else
if ( IWARN!=0 )
end # if
if ( !LSAME( JOBV, 'N' ) )
# interp output 1
    println(io,"V:")
    _nc = MINMN
    _nr = N
    show(io,"text/plain",V[1:_nr,1:_nc])
    println(io,)

end # if
if ( !LSAME( JOBU, 'N' ) )
# interp output 2
    println(io,"U:")
    _nc = MINMN
    _nr = M
    show(io,"text/plain",U[1:_nr,1:_nc])
    println(io,)

end # if
end # if
end # if
    close(f)
end # run_X()


text:
using ArgParse
using MagneticResonanceSignals

function rda_diff(in1file, in2file, outfile; anonymize=true)
    in1header,in1data = load_rda(in1file)
    in2header,in2data = load_rda(in2file)
    outheader = copy(in1header)
    if anonymize
        outheader["PatientName"] = "xxxxxxxxx"
        outheader["PatientID"] = "xxx"
        outheader["PatientBirthDate"] = "20000101"
    end
    save_rda(outfile, outheader, in1data .- in2data)
end

function megapress_argparse(arg_strings; extra_settings...)
    settings = ArgParseSettings(;extra_settings...)
    @add_arg_table settings begin
        "--output", "-o"
            help = "Output file name.  Default is autogenerated based on the input file names."
        "input_files"
            required=true
            nargs=2
            help = "Two input .rda files"
    end
    parse_args(arg_strings, settings)
end

# To test
# megapress_argparse(["a", "b", "c"], exc_handler=ArgParse.debug_handler)

parsed_args = megapress_argparse(ARGS)

rda_diff(parsed_args["input_files"][1], parsed_args["input_files"][2], parsed_args["output"])

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
using Test, MLToolkit
using Distributions: Beta
using Statistics: mean
using StatsFuns: logit, logistic
using Flux: gpu

@testset "Bernoulli" begin
    n_randtests = 3
    n_samples = 5_000
    d, n = 2, 10
    atol = 0.015

    @testset "Bernoulli" begin
        for _ = 1:n_randtests
            p = rand(d, n) |> gpu
            b = Bernoulli(p)
            @test mean(b) == p
            xs = [rand(b) for _ = 1:n_samples]
            @test mean(xs) ≈ mean(b) atol=atol * d * n

            @test logpdf.(Bernoulli.(p), xs[1]) ≈ logpdf(b, xs[1])

            p2 = rand(d, n) |> gpu
            b2 = Bernoulli(p2)
            kl_mc = mean(logpdf.(b, xs) - logpdf.(b2, xs))
            @test kldiv(b, b2) ≈ kl_mc atol=atol * d * n
        end

        @warn "`rand(b::Bernoulli)` is not tested."
    end

    # @testset "BatchBernoulliLogit" begin
    #     for _ = 1:n_randtests
    #         p = Matrix{FT}(rand(Beta(1.0, 1.0), 1, 1))
    #
    #         b = Bernoulli(p[1,1])
    #         bbl = BatchBernoulliLogit{AT}(logit.(p))
    #
    #         x = rand(b, n)
    #         @test cpucopy(vec(logpdf(bbl, AT{FT,2}(reshape(x, 1, n))))) ≈ logpdf.(b, x) atol=5ATOL
    #     end
    # end
end


text:
# validation of properties and code

function check_method(job::CompilerJob)
    isa(job.f, Core.Builtin) && throw(KernelError(job, "function is not a generic function"))

    # get the method
    ms = Base.methods(job.f, job.tt)
    isempty(ms)   && throw(KernelError(job, "no method found"))
    length(ms)!=1 && throw(KernelError(job, "no unique matching method"))
    m = first(ms)

    # kernels can't return values
    if job.kernel
        rt = Base.return_types(job.f, job.tt)[1]
        if rt != Nothing
            throw(KernelError(job, "kernel returns a value of type `$rt`",
                """Make sure your kernel function ends in `return`, `return nothing` or `nothing`.
                   If the returned value is of type `Union{}`, your Julia code probably throws an exception.
                   Inspect the code with `@device_code_warntype` for more details."""))
        end
    end

    return
end

function check_invocation(job::CompilerJob, entry::LLVM.Function)
    # make sure any non-isbits arguments are unused
    real_arg_i = 0
    sig = Base.signature_type(job.f, job.tt)::Type
    for (arg_i,dt) in enumerate(sig.parameters)
        isghosttype(dt) && continue
        real_arg_i += 1

        if !isbitstype(dt)
            param = parameters(entry)[real_arg_i]
            if !isempty(uses(param))
                throw(KernelError(job, "passing and using non-bitstype argument",
                    """Argument $arg_i to your kernel function is of type $dt.
                       That type is not isbits, and such arguments are only allowed when they are unused by the kernel."""))
            end
        end
    end

    return
end

## IR validation

const IRError = Tuple{String, StackTraces.StackTrace, Any} # kind, bt, meta

struct InvalidIRError <: Exception
    job::CompilerJob
    errors::Vector{IRError}
end

const RUNTIME_FUNCTION = "call to the Julia runtime"
const UNKNOWN_FUNCTION = "call to an unknown function"
const POINTER_FUNCTION = "call through a literal pointer"

function Base.showerror(io::IO, err::InvalidIRError)
    print(io, "InvalidIRError: compiling $(signature(err.job)) resulted in invalid LLVM IR")
    for (kind, bt, meta) in err.errors
        print(io, "\nReason: unsupported $kind")
        if meta != nothing
            if kind == RUNTIME_FUNCTION || kind == UNKNOWN_FUNCTION || kind == POINTER_FUNCTION
                print(io, " (call to ", meta, ")")
            end
        end
        Base.show_backtrace(io, bt)
    end
    return
end

function check_ir(job, args...)
    errors = check_ir!(job, IRError[], args...)
    unique!(errors)
    if !isempty(errors)
        throw(InvalidIRError(job, errors))
    end

    return
end

function check_ir!(job, errors::Vector{IRError}, mod::LLVM.Module)
    for f in functions(mod)
        check_ir!(job, errors, f)
    end

    return errors
end

function check_ir!(job, errors::Vector{IRError}, f::LLVM.Function)
    for bb in blocks(f), inst in instructions(bb)
        if isa(inst, LLVM.CallInst)
            check_ir!(job, errors, inst)
        end
    end

    return errors
end

const special_fns = (
    # PTX intrinsics
    "vprintf", "__assertfail", "malloc", "free",
    # libdevice
    "__nvvm_reflect",
    # libcudevrt
    "cudaDeviceSynchronize", "cudaGetParameterBufferV2", "cudaLaunchDeviceV2",
    "cudaCGGetIntrinsicHandle", "cudaCGSynchronize"
)

const libjulia = Ref{Ptr{Cvoid}}(C_NULL)

function check_ir!(job, errors::Vector{IRError}, inst::LLVM.CallInst)
    dest = called_value(inst)
    if isa(dest, LLVM.Function)
        fn = LLVM.name(dest)

        # detect calls to undefined functions
        if isdeclaration(dest) && intrinsic_id(dest) == 0 && !(fn in special_fns)
            # figure out if the function lives in the Julia runtime library
            if libjulia[] == C_NULL
                paths = filter(Libdl.dllist()) do path
                    name = splitdir(path)[2]
                    startswith(name, "libjulia")
                end
                libjulia[] = Libdl.dlopen(first(paths))
            end

            bt = backtrace(inst)
            if Libdl.dlsym_e(libjulia[], fn) != C_NULL
                push!(errors, (RUNTIME_FUNCTION, bt, LLVM.name(dest)))
            else
                push!(errors, (UNKNOWN_FUNCTION, bt, LLVM.name(dest)))
            end
        end
    elseif isa(dest, InlineAsm)
        # let's assume it's valid ASM
    elseif isa(dest, ConstantExpr)
        # detect calls to literal pointers
        # FIXME: can we detect these properly?
        # FIXME: jl_apply_generic and jl_invoke also have such arguments
        if occursin("inttoptr", string(dest))
            # extract the literal pointer
            ptr_arg = first(operands(dest))
            @compiler_assert isa(ptr_arg, ConstantInt) job
            ptr_val = convert(Int, ptr_arg)
            ptr = Ptr{Cvoid}(ptr_val)

            # look it up in the Julia JIT cache
            bt = backtrace(inst)
            frames = ccall(:jl_lookup_code_address, Any, (Ptr{Cvoid}, Cint,), ptr, 0)
            if length(frames) >= 1
                @compiler_assert length(frames) == 1 job frames=frames
                fn, file, line, linfo, fromC, inlined, ip = last(frames)
                push!(errors, (POINTER_FUNCTION, bt, fn))
            else
                fn, file, line, linfo, fromC, inlined, ip = last(frames)
                push!(errors, (POINTER_FUNCTION, bt, nothing))
            end
        end
    end

    return errors
end


text:
# Test time for push! and unshift!
# They are supposed to be amortized constant time

function art(m)
    a = Array(Any,m)
    for i in 1:m
        a[i] = i
    end
    a
end

function partloop(n,m)
    s = 0
    t = @elapsed  for i in 1:n
     a = art(m)
        x = pop!(a)
        push!(a,3)
        s += 1
    end
    t
end

function sartloop(n,m)
    s = 0
    t = @elapsed  for i in 1:n
     a = art(m)
        x = shift!(a)
        unshift!(a,3)
        s += 1
    end
    t
end


text:
#=
							Even Fibonacci numbers

	Each  new term in the Fibonacci sequence is generated by adding the previous
	two terms. By starting with 1 and 2, the first 10 terms will be:

					1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...

	By  considering  the  terms  in  the  Fibonacci sequence whose values do not
	exceed four million, find the sum of the even-valued terms.

								Answer: 4613732
=#

function fibonacci(limit::Int)
    array = [1, 2]
    flag = true
    i = 2

    while flag
        new_value = array[i]+array[i-1]

        if new_value < limit
            push!(array, new_value)
            i += 1
        else
            flag = false
        end
    end

    return array
end

println(sum(filter(iseven, fibonacci(4_000_000))))

max_stars_repo_name: Fazendaaa/project-euler

max_issues_repo_name: Fazendaaa/project-euler

max_forks_repo_name: Fazendaaa/project-euler


text:
module customZipMover

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/CustomZipMover" CustomZipMover(x::Integer, y::Integer, percentage::Integer=100, lineColor::String="663931", lineLightColor::String="9b6157", directory::String="objects/zipmover", isCore::Bool=false, coldLineColor::String="006bb3", coldLineLightColor::String="0099ff", tint::String="ffffff", showLine::Bool=true, fillMiddle::Bool=true, bloomAlpha::Number=1.0, bloomRadius::Number=6.0, rainbow::Bool=false)

const placements = Ahorn.PlacementDict(
	"Custom Zip Mover (Slow, Frost Helper)" => Ahorn.EntityPlacement(
        CustomZipMover,
        "rectangle",
        Dict{String, Any}(),
        function(entity)
            entity.data["nodes"] = [(Int(entity.data["x"]) + Int(entity.data["width"]) + 8, Int(entity.data["y"]))]
			entity.data["percentage"] = 50
			entity.data["lineColor"] = "006bb3"
			entity.data["lineLightColor"] = "0099ff"
			entity.data["directory"] = "objects/FrostHelper/customZipMover/redcog/cold"
        end
    ),
	"Custom Zip Mover (Fast, Frost Helper)" => Ahorn.EntityPlacement(
        CustomZipMover,
        "rectangle",
        Dict{String, Any}(),
        function(entity)
            entity.data["nodes"] = [(Int(entity.data["x"]) + Int(entity.data["width"]) + 8, Int(entity.data["y"]))]
			entity.data["percentage"] = 200
			entity.data["lineColor"] = "e62e00"
			entity.data["lineLightColor"] = "ff5c33"
			entity.data["directory"] = "objects/FrostHelper/customZipMover/redcog"
        end
    ),
	"Custom Zip Mover (Core, Frost Helper)" => Ahorn.EntityPlacement(
        CustomZipMover,
        "rectangle",
        Dict{String, Any}(),
        function(entity)
            entity.data["nodes"] = [(Int(entity.data["x"]) + Int(entity.data["width"]) + 8, Int(entity.data["y"]))]
			entity.data["percentage"] = 200
			entity.data["lineColor"] = "e62e00"
			entity.data["lineLightColor"] = "ff5c33"
			entity.data["directory"] = "objects/FrostHelper/customZipMover/redcog"
			entity.data["isCore"] = true
			entity.data["coldLineColor"] = "006bb3"
			entity.data["coldLineLightColor"] = "0099ff"
        end
    ),
	"Custom Zip Mover (Custom, Frost Helper)" => Ahorn.EntityPlacement(
        CustomZipMover,
        "rectangle",
        Dict{String, Any}(),
        function(entity)
            entity.data["nodes"] = [(Int(entity.data["x"]) + Int(entity.data["width"]) + 8, Int(entity.data["y"]))]
			entity.data["percentage"] = 100
        end
    )
)

Ahorn.nodeLimits(entity::CustomZipMover) = 1, 1

Ahorn.minimumSize(entity::CustomZipMover) = 16, 16
Ahorn.resizable(entity::CustomZipMover) = true, true

function Ahorn.selection(entity::CustomZipMover)
    x, y = Ahorn.position(entity)
    nx, ny = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return [Ahorn.Rectangle(x, y, width, height), Ahorn.Rectangle(nx + floor(Int, width / 2) - 5, ny + floor(Int, height / 2) - 5, 10, 10)]
end

ropeColor = (102, 57, 49) ./ 255
frame = "objects/zipmover/block"
light = "objects/zipmover/light01"
#customCog = String(entity.data["directory"]) + "/cog"

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::CustomZipMover, room::Maple.Room)
    x, y = Ahorn.position(entity)
    nx, ny = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

	#lightSprite = Ahorn.getSprite(get(entity.data, "directory", "objects/zipmover")+"/light01", "Gameplay")
	lightSprite = Ahorn.getSprite(light, "Gameplay")
    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    cx, cy = x + width / 2, y + height / 2
    cnx, cny = nx + width / 2, ny + height / 2

    length = sqrt((x - nx)^2 + (y - ny)^2)
    theta = atan(cny - cy, cnx - cx)

    Ahorn.Cairo.save(ctx)

    Ahorn.translate(ctx, cx, cy)
    Ahorn.rotate(ctx, theta)
	ropeColor = Ahorn.argb32ToRGBATuple(parse(Int, get(entity.data, "lineColor", "663931"), base=16))[1:3] ./ 255
	realRopeColor = (ropeColor..., 1.0)
    Ahorn.setSourceColor(ctx, realRopeColor)
    Ahorn.set_antialias(ctx, 1)
    Ahorn.set_line_width(ctx, 1);

    # Offset for rounding errors
    Ahorn.move_to(ctx, 0, 4 + (theta <= 0))
    Ahorn.line_to(ctx, length, 4 + (theta <= 0))

    Ahorn.move_to(ctx, 0, -4 - (theta > 0))
    Ahorn.line_to(ctx, length, -4 - (theta > 0))

    Ahorn.stroke(ctx)

    Ahorn.Cairo.restore(ctx)

    # Ahorn.drawRectangle(ctx, x + 2, y + 2, width - 4, height - 4, (0.0, 0.0, 0.0, 1.0))
	rawTint = Ahorn.argb32ToRGBATuple(parse(Int, get(entity.data, "tint", "ffffff"), base=16))[1:3] ./ 255
	realTint = (rawTint..., 1.0)
	color = get(entity.data, "color", "Custom")
    if color == "Red"
		Ahorn.drawRectangle(ctx, x + 2, y + 2, width - 4, height - 4, (230.0, 46.0, 0.0, 1.0))
        Ahorn.drawSprite(ctx, "objects/FrostHelper/customZipMover/redcog/cog.png", cnx, cny)
	elseif color == "Core"
		Ahorn.drawRectangle(ctx, x + 2, y + 2, width - 4, height - 4, (230.0, 46.0, 0.0, 1.0))
        Ahorn.drawSprite(ctx, "objects/FrostHelper/customZipMover/redcog/cog.png", cnx, cny)
    elseif color == "Blue"
		Ahorn.drawRectangle(ctx, x + 2, y + 2, width - 4, height - 4, (0.0, 107.0, 179.0, 1.0))
        Ahorn.drawSprite(ctx, "objects/FrostHelper/customZipMover/redcog/cold/cog.png", cnx, cny)
	elseif color == "Black"
		Ahorn.drawRectangle(ctx, x + 2, y + 2, width - 4, height - 4, (0.0, 0.0, 0.0, 0.3))
		Ahorn.drawSprite(ctx, "objects/FrostHelper/customZipMover/blackcog/cog.png", cnx, cny)
	elseif color == "Custom"
		texture = get(entity.data, "directory", "objects/zipmover")
		realTexture = "$texture/cog.png"
		Ahorn.drawRectangle(ctx, x + 2, y + 2, width - 4, height - 4, (0.0, 0.0, 0.0, 1.0))
		Ahorn.drawSprite(ctx, realTexture, cnx, cny, tint = realTint)
    else
		Ahorn.drawRectangle(ctx, x + 2, y + 2, width - 4, height - 4, (0.0, 0.0, 0.0, 1.0))
        Ahorn.drawSprite(ctx, "objects/zipmover/cog.png", cnx, cny)
    end

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y, 8, 0, 8, 8, tint = realTint)
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y + height - 8, 8, 16, 8, 8, tint = realTint)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, x, y + (i - 1) * 8, 0, 8, 8, 8, tint = realTint)
        Ahorn.drawImage(ctx, frame, x + width - 8, y + (i - 1) * 8, 16, 8, 8, 8, tint = realTint)
    end

    Ahorn.drawImage(ctx, frame, x, y, 0, 0, 8, 8, tint = realTint)
    Ahorn.drawImage(ctx, frame, x + width - 8, y, 16, 0, 8, 8, tint = realTint)
    Ahorn.drawImage(ctx, frame, x, y + height - 8, 0, 16, 8, 8, tint = realTint)
    Ahorn.drawImage(ctx, frame, x + width - 8, y + height - 8, 16, 16, 8, 8, tint = realTint)

    Ahorn.drawImage(ctx, lightSprite, x + floor(Int, (width - lightSprite.width) / 2), y, tint = realTint)
end

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper