text:
#===============================================================================

    optical_calibration.jl

    uses prep_normalize() and get_k() with default parameters
    to check validity of calibration data

===============================================================================#

import DataStructures.OrderedDict
import Memento.debug

#===============================================================================
    function definition >>
===============================================================================#

## called by dispatch()
function act(
    ::Type{Val{optical_cal}},
    req                 ::Associative;
    out_format          ::OutputFormat = pre_json_output,

    ## remove MySql dependency
    #
    # db_conn::MySQL.MySQLHandle,
    # exp_id::Integer, ## not used for computation
    # calib_info::Union{Integer,OrderedDict}; ## really used
)
    debug(logger, "at act(::Type{Val{optical_cal}})")

    ## remove MySql dependency
    # calib_info_ori = calib_info
    # calib_info_dict = ensure_ci(db_conn, calib_info_ori)
    # print_v(
    #     println, verbose,
    #     "original: ", calib_info_ori,
    #     "dict: ", calib_info_dict
    # )

    ## get calibration data and use default analysis parameters
    @get_calibration_data_from_req(optical_cal)
    calibration_args = CalibrationParameters()
    #
    ## check validity of data for normalization
    wells = try
        prep_normalize(calibration_data)[2]
    catch err
        return fail(logger, err; bt = true) |> out(out_format)
    end
    #
    ## if there are 2 or more channels then
    ## check validity of data for deconvolution
    if isa(calibration_data, CalibrationData{DualChannel,<: Real})
        result_k = try
            get_k(
                calibration_data,
                calibration_args,
				SVector{length(wells)}(eachindex(wells)),
                wells)
        catch err
            return fail(logger, err; bt = true) |> out(out_format)
        end
        (length(result_k.inv_note) > 0) &&
            return fail(logger, result_k.inv_note) |> out(out_format)
    end ## if
    #
    return OrderedDict(:valid => true) |> out(out_format)
end ## act(::Type{Val{optical_cal}})

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
@enum(
    FNode,

    # leaf nodes
    NEWLINE,
    SEMICOLON,
    WHITESPACE,
    PLACEHOLDER,
    NOTCODE,
    INLINECOMMENT,
    TRAILINGCOMMA,
    TRAILINGSEMICOLON,
    INVERSETRAILINGSEMICOLON,
    KEYWORD,
    LITERAL,
    OPERATOR,
    PUNCTUATION,
    IDENTIFIER,
    # non-leaf nodes
    MacroBlock,
    MacroCall,
    MacroStr,
    Macroname,
    GlobalRefDoc,
    TupleN,
    RefN,
    ModuleN,
    Unary,
    Binary,
    Chain,
    Comparison,
    Conditional,
    Where,
    Vect,
    Braces,
    Brackets,
    Curly,
    Call,
    Parameters,
    Kw,
    Vcat,
    Hcat,
    TypedVcat,
    TypedHcat,
    Row,
    BracesCat,
    TypedComprehension,
    Comprehension,
    Generator,
    Filter,
    Flatten,
    For,
    While,
    If,
    Begin,
    Try,
    Quote,
    Do,
    Let,
    Block,
    BareModule,
    TopLevel,
    StringN,
    Macro,
    FunctionN,
    Struct,
    Mutable,
    Primitive,
    Abstract,
    Return,
    Local,
    Outer,
    Global,
    Const,
    Import,
    Export,
    Using,
    File,
    Quotenode,
    Unknown,
    As,
    NonStdIdentifier, # i.e. var"##iv369"
)

@enum(NestBehavior, AllowNest, AlwaysNest, NeverNest, NeverNestNode)

struct Metadata
    op_kind::Tokens.Kind
    op_dotted::Bool
end

"""
Formatted Syntax Tree
"""
mutable struct FST
    typ::FNode

    # Start and end lines of the node
    # in the original source file.
    startline::Int
    endline::Int

    indent::Int
    len::Int
    val::Union{Nothing,AbstractString}
    nodes::Union{Nothing,Vector{FST}}
    ref::Union{Nothing,Ref{CSTParser.EXPR}}
    nest_behavior::NestBehavior

    # Extra margin caused by parent nodes.
    # i.e. `(f(arg))`
    #
    # `f(arg)` would have `extra_margin` = 1
    # due to `)` after `f(arg)`.
    extra_margin::Int
    line_offset::Int

    metadata::Union{Nothing,Metadata}
end

FST(typ::FNode, cst::CSTParser.EXPR, indent::Int) =
    FST(typ, -1, -1, indent, 0, nothing, FST[], Ref(cst), AllowNest, 0, -1, nothing)

FST(typ::FNode, indent::Int) =
    FST(typ, -1, -1, indent, 0, nothing, FST[], nothing, AllowNest, 0, -1, nothing)

function FST(
    typ::FNode,
    cst::CSTParser.EXPR,
    line_offset::Int,
    startline::Int,
    endline::Int,
    val::AbstractString,
)
    FST(
        typ,
        startline,
        endline,
        0,
        length(val),
        val,
        nothing,
        Ref(cst),
        AllowNest,
        0,
        line_offset,
        nothing,
    )
end

function FST(
    typ::FNode,
    line_offset::Int,
    startline::Int,
    endline::Int,
    val::AbstractString,
)
    FST(
        typ,
        startline,
        endline,
        0,
        length(val),
        val,
        nothing,
        nothing,
        AllowNest,
        0,
        line_offset,
        nothing,
    )
end

@inline function Base.setindex!(fst::FST, node::FST, ind::Int)
    fst.len -= fst.nodes[ind].len
    fst.nodes[ind] = node
    fst.len += node.len
end
@inline Base.getindex(fst::FST, inds...) = fst.nodes[inds...]
@inline Base.lastindex(fst::FST) = length(fst.nodes)
@inline Base.firstindex(fst::FST) = 1
@inline Base.length(fst::FST) = fst.len
@inline function Base.iterate(fst::FST, state = 1)
    if state > length(fst.nodes)
        return nothing
    end
    return fst.nodes[state], state + 1
end

@inline function Base.insert!(fst::FST, ind::Int, node::FST)
    insert!(fst.nodes, ind, node)
    fst.len += node.len
    return
end

@inline Newline(; length = 0, nest_behavior = AllowNest) =
    FST(NEWLINE, -1, -1, 0, length, "\n", nothing, nothing, nest_behavior, 0, -1, nothing)
@inline Semicolon() =
    FST(SEMICOLON, -1, -1, 0, 1, ";", nothing, nothing, AllowNest, 0, -1, nothing)
@inline TrailingComma() =
    FST(TRAILINGCOMMA, -1, -1, 0, 0, "", nothing, nothing, AllowNest, 0, -1, nothing)
@inline TrailingSemicolon() =
    FST(TRAILINGSEMICOLON, -1, -1, 0, 0, "", nothing, nothing, AllowNest, 0, -1, nothing)
@inline InverseTrailingSemicolon() = FST(
    INVERSETRAILINGSEMICOLON,
    -1,
    -1,
    0,
    1,
    ";",
    nothing,
    nothing,
    AllowNest,
    0,
    -1,
    nothing,
)
@inline Whitespace(n) =
    FST(WHITESPACE, -1, -1, 0, n, " "^n, nothing, nothing, AllowNest, 0, -1, nothing)
@inline Placeholder(n) =
    FST(PLACEHOLDER, -1, -1, 0, n, " "^n, nothing, nothing, AllowNest, 0, -1, nothing)
@inline Notcode(startline, endline) =
    FST(NOTCODE, startline, endline, 0, 0, "", nothing, nothing, AllowNest, 0, -1, nothing)
@inline InlineComment(line) =
    FST(INLINECOMMENT, line, line, 0, 0, "", nothing, nothing, AllowNest, 0, -1, nothing)

@inline must_nest(fst::FST) = fst.nest_behavior === AlwaysNest
@inline cant_nest(fst::FST) = fst.nest_behavior === NeverNest
@inline can_nest(fst::FST) = fst.nest_behavior === AllowNest

@inline is_leaf(cst::CSTParser.EXPR) = cst.args === nothing
@inline is_leaf(fst::FST) = fst.nodes === nothing

@inline is_punc(cst::CSTParser.EXPR) = CSTParser.ispunctuation(cst)
@inline is_punc(fst::FST) = fst.typ === PUNCTUATION

@inline is_end(x::CSTParser.EXPR) = x.head === :END && x.val == "end"
@inline is_end(x::FST) = x.typ === KEYWORD && x.val == "end"

@inline is_colon(x::FST) = x.typ === OPERATOR && x.val == ":"

@inline is_comma(fst::FST) = fst.typ === TRAILINGCOMMA || (is_punc(fst) && fst.val == ",")
@inline is_comment(fst::FST) = fst.typ === INLINECOMMENT || fst.typ === NOTCODE

@inline is_circumflex_accent(x::CSTParser.EXPR) =
    CSTParser.isoperator(x) && endswith(CSTParser.valof(x), "^")
@inline is_fwdfwd_slash(x::CSTParser.EXPR) =
    CSTParser.isoperator(x) && endswith(CSTParser.valof(x), "//")

# TODO: fix this
function is_multiline(fst::FST)
    if fst.typ === StringN && fst.endline > fst.startline
        return true
    elseif fst.typ === Vcat && fst.endline > fst.startline
        return true
    elseif fst.typ === TypedVcat && fst.endline > fst.startline
        return true
    elseif fst.typ === MacroStr && fst.endline > fst.startline
        return true
    end
    false
end

@inline is_macrocall(fst::FST) = fst.typ === MacroCall || fst.typ === MacroBlock

function is_macrodoc(cst::CSTParser.EXPR)
    cst.head === :macrocall &&
        length(cst) == 4 &&
        CSTParser.isidentifier(cst[1]) &&
        cst[1].val == "@doc" &&
        (is_str_or_cmd(cst[3]) || is_macrostr(cst[3]))
end

function is_macrodoc(fst::FST)
    fst.typ === GlobalRefDoc && return true
    fst.typ === MacroBlock &&
        fst[1].typ === Macroname &&
        fst[1][1].val == "@doc" &&
        return true
    return false
end

function is_macrostr(cst::CSTParser.EXPR)
    cst.head === :macrocall || return false
    length(cst) > 2 || return false

    local n::Union{Nothing,CSTParser.EXPR}
    n = if CSTParser.isidentifier(cst[1])
        cst[1]
    elseif CSTParser.is_getfield_w_quotenode(cst[1])
        # last quotenode
        qn = cst[1][end]
        # identifier
        qn[end]
    else
        nothing
    end
    n !== nothing || return false

    is_macrostr_identifier(n) || return false

    return is_str_or_cmd(cst[3])
end

function is_macrostr_identifier(cst::CSTParser.EXPR)
    CSTParser.isidentifier(cst) || return false
    # handles odd cases like """M.var"@f";"""
    cst.val !== nothing || return false

    m = match(r"@(.*)_(str|cmd)", cst.val)
    m !== nothing || return false
    val = m.captures[1]

    # r"hello" is parsed as @r_str"hello"
    # if it was originally r"hello" then
    # the number of bytes in the matched string (between @ and _)
    # will be == the span of the node.
    #
    # In this example the span == 1.
    ncodeunits(val) == cst.span
end

function is_call(cst::CSTParser.EXPR)
    t = CSTParser.is_func_call(cst)
    t !== nothing && t && return is_opener(cst[2])
    return false
end

function is_if(cst::CSTParser.EXPR)
    cst.head === :if && cst[1].head === :IF && return true
    cst.head === :elseif && cst[1].head === :ELSEIF && return true
    return false
end

function is_conditional(cst::CSTParser.EXPR)
    CSTParser.isconditional(x)
end

function is_colon_call(cst::CSTParser.EXPR)
    CSTParser.isoperator(cst) && cst.val == ":" && return true
    cst.head == :call && is_colon_call(cst.args[1]) && return true
    return false
end

function is_custom_leaf(fst::FST)
    fst.typ === NEWLINE ||
        fst.typ === SEMICOLON ||
        fst.typ === WHITESPACE ||
        fst.typ === PLACEHOLDER ||
        fst.typ === NOTCODE ||
        fst.typ === INLINECOMMENT ||
        fst.typ === TRAILINGCOMMA ||
        fst.typ === TRAILINGSEMICOLON ||
        fst.typ === INVERSETRAILINGSEMICOLON
end

function is_nothing(cst::CSTParser.EXPR)
    CSTParser.is_nothing(cst) && return true
    cst.val === nothing && cst.args === nothing && return true
    return false
end

"""
Returns whether the first unignored parent of `cst` matches the criteria determined by
`valid`, which is a function that returns a boolean. `ignore` can be used to ignore certain
parent nodes if desired, also a function which returns a boolean.
"""
function parent_is(cst::CSTParser.EXPR, valid; ignore = _ -> false)
    p = cst.parent
    p === nothing && return false
    while p !== nothing && ignore(p)
        p = p.parent
    end
    valid(p)
end

contains_comment(nodes::Vector{FST}) = findfirst(is_comment, nodes) !== nothing
function contains_comment(fst::FST)
    is_leaf(fst) && return false
    contains_comment(fst.nodes)
end

# TODO: Remove once this is fixed in CSTParser.
# https://github.com/julia-vscode/CSTParser.jl/issues/108
function get_args(cst::CSTParser.EXPR)
    if cst.head === :macrocall ||
       cst.head === :typed_vcat ||
       cst.head === :ref ||
       cst.head === :curly ||
       is_call(cst) ||
       cst.head === :typed_comprehension
        return get_args(collect(cst)[2:end])
    elseif cst.head === :where
        # get the arguments in B of `A where B`
        return get_args(collect(cst)[3:end])
    elseif cst.head === :braces ||
           cst.head === :vcat ||
           cst.head === :bracescat ||
           cst.head === :tuple ||
           cst.head === :vect ||
           cst.head === :brackets ||
           cst.head === :parameters ||
           cst.head === :comprehension
        return get_args(cst.args)
    end
    return get_args(cst.args)
end

function get_args(args::Union{Vector{Any},Vector{CSTParser.EXPR}})
    args0 = CSTParser.EXPR[]
    for a in args
        if CSTParser.ispunctuation(a) || CSTParser.is_nothing(a)
            continue
        end
        if a.head === :parameters
            push!(args0, get_args(a)...)
        else
            push!(args0, a)
        end
    end
    args0
end

function get_args(args::Nothing)
    return CSTParser.EXPR[]
end

@inline n_args(x) = length(get_args(x))

"""
    add_node!(
        t::FST,
        n::FST,
        s::State;
        join_lines::Bool = false,
        max_padding::Int = -1,
        override_join_lines_based_on_source::Bool = false,
    )

Appends `n` to `t`.

- `join_lines` if `false` a NEWLINE node will be added and `n` will appear
on the next line, otherwise it will appear on the same line as the previous
node (when printing).
- `max_padding` >= 0 indicates margin of `t` should be based on whether the margin
of `n` + `max_padding` is greater than the current margin of `t`. Otherwise the margin
`n` will be added to `t`.
- `override_join_lines_based_on_source` is only used when `join_lines_based_on_source` option is `true`. In which
`n` is added to `t` as if `join_lines_based_on_source` was false.
"""
function add_node!(
    t::FST,
    n::FST,
    s::State;
    join_lines::Bool = false,
    max_padding::Int = -1,
    override_join_lines_based_on_source::Bool = false,
)
    if n.typ === SEMICOLON
        join_lines = true
        loc =
            s.offset > length(s.doc.text) && t.typ === TopLevel ?
            cursor_loc(s, s.offset - 1) : cursor_loc(s)
        for l = t.endline:loc[1]
            if has_semicolon(s.doc, l)
                n.startline = l
                n.endline = l
                break
            end
        end

        # If there's no semicolon, treat it
        # as formatter node
        if n.startline == -1
            t.len += length(n)
            n.startline = t.endline
            n.endline = t.endline
            push!(t.nodes, n)
            return
        end
    elseif n.typ === TRAILINGCOMMA
        en = t.nodes[end]
        if en.typ === Generator ||
           en.typ === Filter ||
           en.typ === Flatten ||
           en.typ === MacroCall ||
           en.typ === MacroBlock ||
           en.typ === SEMICOLON
            # don't add trailing comma in these cases
        elseif is_comma(en) && t.typ === TupleN && n_args(t.ref[]) == 1
            # preserve comma
        elseif s.opts.trailing_comma === nothing
        elseif !s.opts.trailing_comma
            if is_comma(en)
                t[end] = Whitespace(0)
            end
        elseif is_comma(en)
            t[end] = n
        else
            t.len += length(n)
            n.startline = t.endline
            n.endline = t.endline
            push!(t.nodes, n)
        end
        return
    elseif n.typ === NOTCODE
        n.indent = s.indent
        push!(t.nodes, n)
        return
    elseif n.typ === INLINECOMMENT
        push!(t.nodes, n)
        return
    elseif is_custom_leaf(n)
        t.len += length(n)
        n.startline = t.endline
        n.endline = t.endline
        push!(t.nodes, n)
        return
    end

    if n.typ === Block && length(n) == 0
        push!(t.nodes, n)
        return
    elseif n.typ === Parameters
        # unpack Parameters arguments into the parent node
        if t.ref !== nothing && n_args(t.ref[]) == n_args(n.ref[])
            # There are no arguments prior to params
            # so we can remove the initial placeholder.
            idx = findfirst(n -> n.typ === PLACEHOLDER, t.nodes)
            idx !== nothing && (t[idx] = Whitespace(0))
        end
        add_node!(t, Semicolon(), s)

        if length(n.nodes) > 0
            nws = 1
            if (t.typ === Curly || t.typ === Where || t.typ === BracesCat) &&
               !s.opts.whitespace_typedefs
                nws = 0
            end
            multi_arg = t.ref !== nothing && n_args(t.ref[]) > 0
            multi_arg ? add_node!(t, Placeholder(nws), s) : add_node!(t, Whitespace(nws), s)
        end
        for nn in n.nodes
            add_node!(t, nn, s, join_lines = true)
        end
        return
    elseif s.opts.import_to_using && n.typ === Import && t.typ !== MacroBlock
        usings = import_to_usings(n, s)
        if length(usings) > 0
            for (i, nn) in enumerate(usings)
                add_node!(t, nn, s; join_lines = false, max_padding = 0)
            end
            return
        end
    elseif n.typ === Binary && n[1].typ === Binary && n[1][end].typ === Where
        # normalize FST representation for Where
        binaryop_to_whereop!(n, s)
    end

    if length(t.nodes) == 0
        t.startline = n.startline
        t.endline = n.endline
        t.len += length(n)
        t.line_offset = n.line_offset
        push!(t.nodes, n)
        return
    end

    # if `max_padding` >= 0 `n` should appear on the next line
    # even if it's contrary on the original source.
    if s.opts.join_lines_based_on_source &&
       !override_join_lines_based_on_source &&
       max_padding == -1 &&
       !(
           is_comma(n) ||
           is_block(t) ||
           t.typ === FunctionN ||
           t.typ === Macro ||
           is_typedef(t) ||
           t.typ === ModuleN ||
           t.typ === BareModule ||
           is_end(n)
       )
        # join based on position in original file
        join_lines = t.endline == n.startline
    end

    if !is_prev_newline(t.nodes[end])
        current_line = t.endline
        notcode_startline = current_line + 1
        notcode_endline = n.startline - 1
        nt = t.nodes[end].typ

        if notcode_startline <= notcode_endline
            # If there are comments in between node elements
            # nesting is forced in an effort to preserve them.

            rm_block_nl =
                s.opts.remove_extra_newlines &&
                t.typ !== ModuleN &&
                (n.typ === Block || is_end(n))

            nest = true
            if remove_empty_notcode(t) || rm_block_nl
                nest = false
                for l = notcode_startline:notcode_endline
                    if hascomment(s.doc, l)
                        nest = true
                        break
                    end
                end
            end

            t.nest_behavior = AlwaysNest

            # If the previous node type is WHITESPACE - reset it.
            # This fixes cases similar to the one shown in issue #51.
            nt === WHITESPACE && (t.nodes[end] = Whitespace(0))

            if hascomment(s.doc, current_line)
                add_node!(t, InlineComment(current_line), s)
            end

            if nt !== PLACEHOLDER
                add_node!(t, Newline(nest_behavior = AlwaysNest), s)
            elseif hascomment(s.doc, current_line) && nt === PLACEHOLDER
                # swap PLACEHOLDER (will be NEWLINE) with INLINECOMMENT node
                idx = length(t.nodes)
                t.nodes[idx-1], t.nodes[idx] = t.nodes[idx], t.nodes[idx-1]
            end

            if nest
                add_node!(t, Notcode(notcode_startline, notcode_endline), s)
                add_node!(t, Newline(nest_behavior = AlwaysNest), s)
            end
        elseif !join_lines
            if hascomment(s.doc, current_line) && current_line != n.startline
                add_node!(t, InlineComment(current_line), s)
            end
            add_node!(t, Newline(nest_behavior = AlwaysNest), s)
        elseif nt === PLACEHOLDER &&
               current_line != n.startline &&
               hascomment(s.doc, current_line)
            t.nest_behavior = AlwaysNest
            add_node!(t, InlineComment(current_line), s)
            # swap PLACEHOLDER (will be NEWLINE) with INLINECOMMENT node
            idx = length(t.nodes)
            t.nodes[idx-1], t.nodes[idx] = t.nodes[idx], t.nodes[idx-1]
        elseif nt === WHITESPACE &&
               hascomment(s.doc, current_line) &&
               current_line != n.startline
            # rely on the whitespace tracked for the inline comment
            t.nodes[end] = Whitespace(0)
            add_node!(t, InlineComment(current_line), s)
            add_node!(t, Newline(nest_behavior = AlwaysNest), s)
        end
    end

    if n.startline < t.startline || t.startline == -1
        t.startline = n.startline
    end
    if n.endline > t.endline || t.endline == -1
        t.endline = n.endline
    end

    if !join_lines && is_end(n)
        # end keyword isn't useful w.r.t margin lengths
    elseif t.typ === StringN
        # The length of this node is the length of
        # the longest string. The length of the string is
        # only considered "in the positive" when it's past
        # the hits the initial """ offset, i.e. `t.indent`.
        t.len = max(t.len, n.indent + length(n) - t.indent)
    elseif is_multiline(n)
        is_iterable(t) && n_args(t.ref[]) > 1 && (t.nest_behavior = AlwaysNest)
        t.len += length(n)
    elseif max_padding >= 0
        t.len = max(t.len, length(n) + max_padding)
    else
        t.len += length(n)
    end
    push!(t.nodes, n)
    nothing
end

@inline function is_prev_newline(fst::FST)
    if fst.typ === NEWLINE
        return true
    elseif is_leaf(fst) || length(fst.nodes) == 0
        return false
    end
    is_prev_newline(fst[end])
end

"""
    `length_to(x::FST, ntyps; start::Int = 1)`

Returns the length to any node type in `ntyps` based off the `start` index.
"""
@inline function length_to(fst::FST, ntyps; start::Int = 1)
    fst.typ in ntyps && return 0, true
    is_leaf(fst) && return length(fst), false
    len = 0
    for i = start:length(fst.nodes)
        l, found = length_to(fst.nodes[i], ntyps)
        len += l
        found && return len, found
    end
    return len, false
end

@inline is_closer(fst::FST) =
    fst.typ === PUNCTUATION && (fst.val == "}" || fst.val == ")" || fst.val == "]")
@inline is_closer(cst::CSTParser.EXPR) =
    cst.head === :RBRACE || cst.head === :RPAREN || cst.head === :RSQUARE

@inline is_opener(fst::FST) =
    fst.typ === PUNCTUATION && (fst.val == "{" || fst.val == "(" || fst.val == "[")
@inline is_opener(cst::CSTParser.EXPR) =
    cst.head === :LBRACE || cst.head === :LPAREN || cst.head === :LSQUARE

function is_iterable(x::CSTParser.EXPR)
    x.head === :tuple && return true
    x.head === :vect && return true
    x.head === :vcat && return true
    x.head === :braces && return true
    is_call(x) && return true
    x.head === :curly && return true
    x.head === :comprehension && return true
    x.head === :typed_comprehension && return true
    x.head === :macrocall && return true
    x.head === :brackets && return true
    x.head === :ref && return true
    x.head === :typed_vcat && return true
    x.head === :import && return true
    x.head === :using && return true
    x.head === :export && return true
    return false
end

function is_iterable(x::FST)
    is_named_iterable(x) && return true
    is_unnamed_iterable(x) && return true
    is_import_expr(x) && return true
    return false
end

function is_unnamed_iterable(x::FST)
    x.typ === TupleN && return true
    x.typ === Vect && return true
    x.typ === Vcat && return true
    x.typ === Braces && return true
    x.typ === Comprehension && return true
    x.typ === Brackets && return true
    return false
end

function is_named_iterable(x::FST)
    x.typ === Call && return true
    x.typ === Curly && return true
    x.typ === TypedComprehension && return true
    x.typ === MacroCall && return true
    x.typ === RefN && return true
    x.typ === TypedVcat && return true
    return false
end

function is_import_expr(x::FST)
    x.typ === Import && return true
    x.typ === Using && return true
    x.typ === Export && return true
    return false
end

"""
Returns whether `fst` can be an iterable argument. For example in
the case of a function call, which is of type `Call`:

```julia
(a, b, c; k1=v1)
```

This would return `true` for `a`, `b`, `c` and `k1=v1` and `false` for all other nodes.
"""
function is_iterable_arg(fst::FST)
    fst.typ === PUNCTUATION && return false
    fst.typ === KEYWORD && return false
    fst.typ === OPERATOR && return false
    is_custom_leaf(fst) && return false
    return true
end

function is_comprehension(x::CSTParser.EXPR)
    x.head === :comprehension && return true
    x.head === :typed_comprehension && return true
    return false
end

function is_comprehension(x::FST)
    x.typ === Comprehension && return true
    x.typ === TypedComprehension && return true
    return false
end

function is_block(x::CSTParser.EXPR)
    is_if(x) && return true
    x.head === :do && return true
    x.head === :try && return true
    (x.head === :block && length(x) > 1 && x[1].head == :BEGIN) && return true
    x.head === :for && return true
    x.head === :while && return true
    x.head === :let && return true
    (x.head === :quote && x[1].head == :QUOTE) && return true
    return false
end

function is_block(x::FST)
    x.typ === If && return true
    x.typ === Do && return true
    x.typ === Try && return true
    x.typ === Begin && return true
    x.typ === For && return true
    x.typ === While && return true
    x.typ === Let && return true
    x.typ === Quote && x[1].val == "quote" && return true
    return false
end

function is_typedef(fst::FST)
    fst.typ === Struct && return true
    fst.typ === Mutable && return true
    fst.typ === Primitive && return true
    fst.typ === Abstract && return true
    return false
end

function is_opcall(x::CSTParser.EXPR)
    is_binary(x) && return true
    x.head === :comparison && return true
    is_chain(x) && return true
    # is_unary(x) && return true
    # Brackets are often mixed with operators
    # so kwargs are propagated through its related
    # functions
    x.head === :brackets && return true
    return false
end

function is_opcall(x::FST)
    x.typ === Binary && return true
    x.typ === Comparison && return true
    x.typ === Chain && return true
    # x.typ === Unary && return true
    # Brackets are often mixed with operators
    # so kwargs are propagated through its related
    # functions
    x.typ === Brackets && return true
    return false
end

# Generator typ
# (x for x in 1:10)
# (x for x in 1:10 if x % 2 == 0)
function is_gen(x::CSTParser.EXPR)
    x.head === :generator && return true
    x.head === :filter && return true
    x.head === :flatten && return true
    # x.head === :brackets && return true
    return false
end

function is_gen(x::FST)
    x.typ === Generator && return true
    x.typ === Filter && return true
    x.typ === Flatten && return true
    # x.typ === Brackets && return true
    return false
end

function is_binary(x::CSTParser.EXPR)
    (CSTParser.isbinarycall(x) || CSTParser.isbinarysyntax(x)) &&
        length(x) == 3 &&
        CSTParser.isoperator(x[2])
end

function is_unary(x::CSTParser.EXPR)
    CSTParser.isunarycall(x) || CSTParser.isunarysyntax(x)
end

function is_chain(x::CSTParser.EXPR)
    CSTParser.ischainedcall(x)
end

function is_assignment(x::FST)
    if x.typ === Binary && is_assignment(op_kind(x))
        return true
    end

    if (
        x.typ === Const ||
        x.typ === Local ||
        x.typ === Global ||
        x.typ === Outer ||
        x.typ === MacroBlock
    ) && is_assignment(x[end])
        return true
    end

    return false
end

is_assignment(kind::Tokens.Kind) = CSTParser.precedence(kind) == CSTParser.AssignmentOp
function is_assignment(cst::CSTParser.EXPR)
    op = get_binary_op(cst)
    op === nothing && return false
    precedence(op) == CSTParser.AssignmentOp
end
is_assignment(::Nothing) = false

function is_pairarrow(cst::CSTParser.EXPR)
    op = get_binary_op(cst)
    op === nothing && return false
    CSTParser.is_pairarrow(op)
end

function precedence(cst::CSTParser.EXPR)
    CSTParser.isoperator(cst) || return 0
    val = CSTParser.isdotted(cst) ? cst.val[2:end] : cst.val
    CSTParser.get_prec(val)
end

function is_function_or_macro_def(cst::CSTParser.EXPR)
    CSTParser.defines_function(cst) && return true
    cst.head === :macro && return true
    cst.head === :where && return true
    return false
end

function nest_block(cst::CSTParser.EXPR)
    is_if(cst) && return true
    cst.head === :do && return true
    cst.head === :try && return true
    cst.head === :for && return true
    cst.head === :while && return true
    cst.head === :let && return true
    return false
end

function remove_empty_notcode(fst::FST)
    is_iterable(fst) && return true
    fst.typ === Binary && return true
    fst.typ === Conditional && return true
    fst.typ === Comparison && return true
    fst.typ === Chain && return true
    return false
end

function get_binary_op(cst::CSTParser.EXPR)
    if cst.head == :call
        return cst.args[1]
    elseif length(cst) == 3 && CSTParser.isoperator(cst[2])
        return cst[2]
    elseif CSTParser.isoperator(cst.head)
        return cst.head
    else
        return nothing
    end
end

"""
`cst` is assumed to be a single child node. Returns true if the node is of the syntactic form `{...}, [...], or (...)`.
"""
function unnestable_node(cst::CSTParser.EXPR)
    cst.head === :tuple && return true
    cst.head === :vect && return true
    cst.head === :braces && return true
    cst.head === :bracescat && return true
    cst.head === :comprehension && return true
    cst.head === :brackets && return true
    return false
end

function is_binaryop_nestable(::S, cst::CSTParser.EXPR) where {S<:AbstractStyle}
    CSTParser.defines_function(cst) && is_unary(cst[1]) && return true
    if is_assignment(cst) || is_pairarrow(cst)
        return !is_str_or_cmd(cst[3])
    end
    true
end

function nest_rhs(cst::CSTParser.EXPR)::Bool
    if CSTParser.defines_function(cst)
        rhs = cst[3]
        rhs.head === :block && (rhs = rhs[1])
        return nest_block(rhs)
    end
    false
end

function op_kind(cst::CSTParser.EXPR)
    if is_binary(cst) || cst.head === :comparison || is_chain(cst)
        return tokenize(cst[2].val)
    elseif is_unary(cst)
        op = CSTParser.isoperator(cst[1]) ? cst[1] : cst[2]
        return tokenize(op.val)
    end
    return nothing
end

function op_kind(fst::FST)
    if fst.typ === OPERATOR
        return fst.metadata.op_kind
    elseif is_opcall(fst)
        idx = findfirst(n -> n.typ === OPERATOR, fst.nodes)
        idx === nothing && return nothing
        return fst[idx].metadata.op_kind
    end
    return nothing
end

function _valid_parent_node_for_standalone_circuit(n::CSTParser.EXPR)
    n.head === :brackets && return false
    n.head === :macrocall && return false
    n.head === :return && return false
    is_if(n) && return false
    n.head === :block && is_assignment(n.parent) && return false
    is_binary(n) && is_assignment(n) && return false
    return true
end
_valid_parent_node_for_standalone_circuit(n::Nothing) = true

function _ignore_node_for_standalone_circuit(n::CSTParser.EXPR)
    n.head === :brackets && return false
    is_if(n) && return false
    n.head === :block && return false
    n.head === :macrocall && return false
    n.head === :return && return false
    is_binary(n) && is_assignment(n) && return false
    return true
end

"""
    is_standalone_shortcircuit(cst::CSTParser.EXPR)

Returns `true` if the `cst` is a short-circuit expression (uses `&&`, `||`)
and is *standalone*, meaning it's not directly associated with another statement or
expression.

### Examples

```julia
# this IS a standalone short-circuit
a && b

# this IS NOT a standalone short-circuit
if a && b
end

# this IS NOT a standalone short-circuit
var = a && b

# this IS NOT a standalone short-circuit
@macro a && b

# operation inside parenthesis IS NOT a standalone short-circuit
# operation outside parenthesis IS a standalone short-circuit
(a && b) && c
```
"""
function is_standalone_shortcircuit(cst::CSTParser.EXPR)
    val = cst[2].val
    (val == "&&" || val == "||") || return false

    return parent_is(
        cst,
        _valid_parent_node_for_standalone_circuit,
        ignore = _ignore_node_for_standalone_circuit,
    )
end

"""
    eq_to_in_normalization!(fst::FST, always_for_in::Bool)
    eq_to_in_normalization!(fst::FST, always_for_in::Nothing)

Transforms

```julia
for i = iter body end

=>

for i in iter body end
```

AND

```julia
for i in 1:10 body end

=>

for i = 1:10 body end
```

`always_for_in=nothing` disables this normalization behavior.

- https://github.com/domluna/JuliaFormatter.jl/issues/34
"""
function eq_to_in_normalization!(fst::FST, always_for_in::Bool)
    if fst.typ === Binary
        idx = findfirst(n -> n.typ === OPERATOR, fst.nodes)
        idx === nothing && return
        op = fst[idx]

        if always_for_in
            op.val = "in"
            op.len = length(op.val)
            return
        end

        if op.val == "=" && op_kind(fst[end]) !== Tokens.COLON
            op.val = "in"
            op.len = length(op.val)
        elseif op.val == "in" && op_kind(fst[end]) === Tokens.COLON
            op.val = "="
            op.len = length(op.val)
        end
    elseif fst.typ === Block || fst.typ === Brackets
        for n in fst.nodes
            eq_to_in_normalization!(n, always_for_in)
        end
    end
end

eq_to_in_normalization!(::FST, ::Nothing) = nothing

ext: jl


text:
using Base.Test
using PValueAdjust

newPvalues = padjust([0.05, 0.03, 0.01, 0.5], BenjaminiHochberg)
@test_approx_eq newPvalues [0.06666666667, 0.06, 0.04, 0.5]

@test_approx_eq padjust([1, 0.5, 0.1, 0.001, 0.0001, 0.002], BenjaminiHochberg) [1.0, 0.6, 0.15, 0.003, 0.0006, 0.004]
@test_approx_eq padjust([0.05, 0.05, 0.05], BenjaminiHochberg) [0.05, 0.05, 0.05]

@test all(padjust([1,1,1], BenjaminiHochberg) .<= 1)

@test_throws DomainError padjust([-1], BenjaminiHochberg)
@test_throws DomainError padjust([2], BenjaminiHochberg)

@test padjust([0.5], BenjaminiHochberg) == [0.5]


text:
@testset "CST element" begin
    # test stiffness matrix
    dim = 2
    material = LinearElasticity{2}(K=1.25, G=0.4)

    ip = Lagrange{dim, RefTetrahedron, 1}()
    qr = QuadratureRule{dim, RefTetrahedron}(1)
    cv = CellVectorValues(qr, ip)

    xe = [Vec(6., -4.), Vec(5., 3.), Vec(-4., -1.)]

    thickness = 2.0

    ke = Matrix{Float64}(undef, 6, 6)
    fe = Vector{Float64}(undef, 6)

    qr_face = QuadratureRule{dim-1, RefTetrahedron}(1)
    fv = FaceVectorValues(qr_face, ip)

    nodes = Node.(xe)
    cells = [Triangle((1,2,3))]
    grid = Grid(cells, nodes)
    addfaceset!(grid, "Γ", Set((FaceIndex(1,1),)))

    tₚ = Vec(1.0, 1.0)

    element_routine!(Primal(), ke, fe, cv, xe, material, thickness, fv, grid, 1, tₚ, "Γ")

    ke_reference = [
        0.9095 -0.7433 -0.2179 0.4259 -0.6915 0.3174
        -0.7433 2.2515 -0.1575 -2.3239 0.9007 0.0724
        -0.2179 -0.1575 0.8366 0.6194 -0.6187 -0.4619
        0.4259 -2.3239 0.6194 2.7154 -1.0453 -0.3915
        -0.6915 0.9007 -0.6187 -1.0453 1.3102 0.1445
        0.3174 0.0724 -0.4619 -0.3915 0.1445 0.3192
        ]
    @test round.(ke, digits=4) ≈ ke_reference

    fe_reference = tₚ*thickness*norm(xe[2]-xe[1]) / 2
    @test fe[1:4] ≈ vcat(fe_reference, fe_reference)
end


text:
#=
reverser.jl
based on
https://stackoverflow.com/questions/27411401/julia-reverse-n-dimensional-arrays
=#

export reverser

"""
    y = reverser(x, dims)
reverse array along specified dimensions (or all if unspecified)
"""
function reverser(x::AbstractArray, dims::AbstractVector{<:Int})
    y = copy(x)
    for d in dims
        y = reverse(y, dims=d)
    end
    return y
end

reverser(x::AbstractArray) = reverser(x, 1:ndims(x)) # all dimensions
reverser(x::AbstractArray, d::Int) = reverser(x, [d])


text:
using PyPlot
using JLD2
using DelimitedFiles

n = 30
@load "data/dat0-1.jld2" rcvv U sigma tau rcvi rcvj φ ψ
rcvv0 = rcvv
sigma0 = sigma
tau0 = tau
U0 = U
Phi0 = φ
Psi0 = ψ
rcvi0 = rcvi
rcvj0 = rcvj

Uk = zeros(size(U0,1), n, n)
Phik = zeros(size(U0,1), n, n)
Psik = zeros(size(U0,1), n, n)

@load "data/dat0-4.jld2" rcvv U sigma tau rcvi rcvj φ ψ
rcvv1 = rcvv
sigma1 = sigma[1:end-1, 1:end-1]
tau1 = tau[1:end-1, 1:end-1]
φ1 = φ
ψ1 = ψ

U1 = U
Uk[:,1:n÷2,1:n÷2] = U1
Phik[:,1:n÷2,1:n÷2] = φ1
Psik[:,1:n÷2,1:n÷2] = ψ1

rcvi1 = rcvi
rcvj1 = rcvj

@load "data/dat1-4.jld2" rcvv U sigma tau rcvi rcvj φ ψ
rcvv2 = rcvv
sigma2 = sigma[1:end-1, 2:end]
tau2 = tau[1:end-1, 2:end]
U2 = U
Uk[:,1:n÷2,n÷2+1:end] = U2
rcvi2 = rcvi
rcvj2 = rcvj
φ2 = φ
ψ2 = ψ

Phik[:,1:n÷2,n÷2+1:end] = φ2
Psik[:,1:n÷2,n÷2+1:end] = ψ2

@load "data/dat2-4.jld2" rcvv U sigma tau rcvi rcvj φ ψ
rcvv3 = rcvv
sigma3 = sigma[2:end, 1:end-1]
tau3 = tau[2:end, 1:end-1]
U3 = U
Uk[:,n÷2+1:end, 1:n÷2] = U3
rcvi3 = rcvi
rcvj3 = rcvj
φ3 = φ
ψ3 = ψ

Phik[:,n÷2+1:end, 1:n÷2] = φ3
Psik[:,n÷2+1:end, 1:n÷2] = ψ3

@load "data/dat3-4.jld2" rcvv U sigma tau rcvi rcvj φ ψ
rcvv4 = rcvv
sigma4 = sigma[2:end, 2:end]
tau4 = tau[2:end, 2:end]
U4 = U
Uk[:,n÷2+1:end, n÷2+1:end] = U4
rcvi4 = rcvi
rcvj4 = rcvj
φ4 = φ
ψ4 = ψ

Phik[:,n÷2+1:end, n÷2+1:end] = φ4
Psik[:,n÷2+1:end, n÷2+1:end] = ψ4

@show maximum(abs.(U0 - Uk))

S = [sigma1 sigma2; sigma3 sigma4]
T = [tau1 tau2; tau3 tau4]

@show maximum(abs.(sigma0-S))
@show maximum(abs.(tau0-T))

rcvv_mpi = [rcvv1 rcvv2]

idx0 = @. (rcvj0 - 1) * n + rcvi0
rcvv0_ref = reshape(U0, size(U0,1), :)[:, idx0]

idx1 = @. (rcvj1 - 1) * (n÷2) + rcvi1
rcvv1_ref = reshape(U1, size(U1,1), :)[:, idx1]

@show maximum(abs.(rcvv_mpi - rcvv0))

# writedlm(stdout, Uk[8,1:15,16:end])

# writedlm(stdout, U0[8,1:15,16:end])
# close("all")
# pcolormesh(rcvv_mpi, cmap="gray")
# colorbar()
# savefig("receiver_mpi.png")


text:
# # Turn to Gold

# ### Turn blocks to gold by hitting them
# Demonstration of pollBlockHits()
# the function pollBlockHits() returns an Array of hit Events
# A hit event is a tuple, say `((5, 4, 5), 2, 17)`
# `(5, 4, 5)` are the coordinates to the block which was hit
# `2` is the face of the block which was hit
# `17` is the entityId of the Entity which was responsible for the hit

using Dates
using PiCraft

clearEvents() # Clear all previous events
t = now() # Mark the current time

post("60 seconds to turn blocks to gold.")

while (now() - t).value/1000 < 60.0 # Continue the loop if 60 seconds has not elapsed
    for hit in pollBlockHits()
        setBlock(hit[1], Block(41)) # hit[1] are the coordinates to the Block which was hit
    end
end

# Post a message to screen when done
post("Time Over!")


text:
import DataFrames
import Random
import StatsBase

"""
"""
function split_data(
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        split::Real,
        )
    result = split_data(
        Random.GLOBAL_RNG,
        features_df,
        labels_df,
        split,
        )
    return result
end

"""
"""
function split_data(
        rng::Random.AbstractRNG,
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        split::Real,
        )
    #
    if !(0 < split < 1)
        error("split must be >0 and <1")
    end
    if size(features_df, 1) != size(labels_df, 1)
        error("features_df and labels_df do not have the same number of rows")
    end
    #
    num_rows = size(features_df, 1)
    num_partition_1 = round(Int, split * num_rows)
    num_partition_2 = num_rows - num_partition_1
    #
    allrows = convert(Array, 1:num_rows)
    partition_1_rows = StatsBase.sample(
        rng,
        allrows,
        num_partition_1;
        replace = false,
        )
    partition_2_rows = setdiff(allrows, partition_1_rows)
    #
    partition_1_features_df = features_df[partition_1_rows, :]
    partition_2_features_df = features_df[partition_2_rows, :]
    #
    partition_1_labels_df = labels_df[partition_1_rows, :]
    partition_2_labels_df = labels_df[partition_2_rows, :]
    #
    return partition_1_features_df,
        partition_1_labels_df,
        partition_2_features_df,
        partition_2_labels_df
end


text:
include("../ranger.jl")

using Base.Math: deg2rad

using .Ranger

const RGeo = Ranger.Geometry
const RMath = Ranger.Math
const RAnim = Ranger.Animation

println("Basic Ranger test")

p = RGeo.Point{Float64}()
println("p: $p")

aft = RMath.AffineTransform{Float64}()
println("aft: $aft")

RMath.set!(aft, 1.1, 2.2, 3.3, 4.4, 5.5, 6.6)
println("aft: $aft")

RMath.to_identity!(aft)
println("aft: $aft")

RMath.translate!(aft, 5.0, 10.0)
println("aft: $aft")

v = RGeo.Point{Float64}(1.0, 2.0)
println("v: $v")
RMath.transform!(aft, v)

println("v: $v")

println("transform by components")
cv = RMath.transform!(aft, 5.0, 5.0)

println("cv: $cv")

println("--------------------------")
println("rotate +x vector")
v = RGeo.Point{Float64}(1.0, 0.0)
println("new v: $v")
rott = RMath.AffineTransform{Float64}()
radians = deg2rad(45.0);

RMath.make_rotate!(rott, radians)
RMath.transform!(rott, v)
println("v rotated: $v")

rv = RMath.transform!(rott, 1.0, 0.0)
println("rv rotated: $rv")

println("--------------------------")
p0 = RGeo.Point{Float64}(5.0, 5.0)
p1 = RGeo.Point{Float64}(2.5, 10.0)
p2 = RGeo.Point{Float64}(7.5, 10.0)

aabb = RGeo.AABB{Float64}()
println("aabb: $aabb")
RGeo.set!(aabb, p0, p1, p2)
println("aabb: $aabb")

println("--------------------------")
verts = [p0, p1, p2]
# println(verts)

aabb = RGeo.AABB{Float64}()
println("A Verts aabb: $aabb")
RGeo.set!(aabb, verts)
println("B Verts aabb: $aabb")

println("intersection rectangle --------------------------")
intersect = RGeo.Rectangle{Float64}()
rectA = RGeo.Rectangle{Float64}(0.0, 0.0, 10.0, 10.0)
rectB = RGeo.Rectangle{Float64}(5.0, 5.0, 20.0, 20.0)

RGeo.intersect!(intersect, rectA, rectB)

println(intersect)

println("intersects True rectangle --------------------------")
intersects = RGeo.intersects(rectA, rectB)

println(intersects)

println("intersects False rectangle --------------------------")
rectC = RGeo.Rectangle{Float64}(0.0, 0.0, 10.0, 10.0)
rectD = RGeo.Rectangle{Float64}(15.0, 15.0, 25.0, 25.0)
intersects = RGeo.intersects(rectC, rectD)

println(intersects)

println("bounds rectangle --------------------------")
bounds = RGeo.Rectangle{Float64}()
RGeo.bounds!(bounds, rectA, rectB)

println(bounds)

println("contains True point --------------------------")
contains = RGeo.contains_point(rectA, p0)

println(contains)

println("contains False point --------------------------")
p4 = RGeo.Point{Float64}(25.0, 25.0)
contains = RGeo.contains_point(rectA, p4)

println(contains)

println("motion ----------------------------------------")
motion = RAnim.AngularMotion{Float64}()
RAnim.set!(motion, 1.0, 2.0, 0.0)

value = RAnim.interpolate!(motion, 0.0)
println(value)
value = RAnim.interpolate!(motion, 0.1)
println(value)
value = RAnim.interpolate!(motion, 0.5)
println(value)
value = RAnim.interpolate!(motion, 1.0)
println(value)

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
# Visually inspect ContinuousSpace `nearby_ids`
using Agents
using AgentsPlots
AgentsPlots.Plots.pyplot()

mutable struct Ag <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
end

ε = 0.2
s = ContinuousSpace((10, 10), ε)
model = ABM(Ag, s)
for i in 1:800
    add_agent!(model)
end

r = 2
function ac(a)
    nids = nearby_ids(model[agentid], model, r, exact=false)
    if a.id == agentid
        return :red
 elseif a.id in nids
        return :blue
    else
        return :yellow
    end
end

function circleShape(h, k, r)
    θ = LinRange(0, 2π, 500)
    h .+ r*sin.(θ), k .+ r .* cos.(θ)
end

# %% --- INEXACT ---
agentid = rand(1:800)
a = model[agentid]
δ = Agents.distance_from_cell_center(a.pos, model)
grid_r = (r+δ)/model.space.spacing
focal_cell = Agents.pos2cell(a.pos, model)
allcells = Agents.grid_space_neighborhood(CartesianIndex(focal_cell), model, grid_r)
search_region = [(a .* s.spacing) .- (s.spacing/2) for a in allcells]

p = plotabm(model, as=4, ac=ac, grid = (:both, :dot, 1, 0.9), xticks=(0:s.spacing:s.extent[1]), yticks=(0:s.spacing:s.extent[2]), size=(1000, 1000))
scatter!(search_region; markershape=:square, markersize=8, markerstrokewidth = 0, markeralpha = 0.2, markercolor=:grey)

plot!(p, circleShape(a.pos[1], a.pos[2], r), seriestype = :shape, lw = 0.5,
      c = :blue, linecolor = :black, legend = false, fillalpha = 0.2, aspect_ratio=1)
# Add all possible mirrors
plot!(p, circleShape(a.pos[1]-s.extent[1], a.pos[2], r), seriestype = :shape, lw = 0.5,
      c = :blue, linecolor = :black, legend = false, fillalpha = 0.2, aspect_ratio=1)
plot!(p, circleShape(a.pos[1]-s.extent[1], a.pos[2]-s.extent[2], r), seriestype = :shape, lw = 0.5,
      c = :blue, linecolor = :black, legend = false, fillalpha = 0.2, aspect_ratio=1)
plot!(p, circleShape(a.pos[1], a.pos[2]-s.extent[2], r), seriestype = :shape, lw = 0.5,
      c = :blue, linecolor = :black, legend = false, fillalpha = 0.2, aspect_ratio=1)
plot!(p, circleShape(a.pos[1]+s.extent[1], a.pos[2], r), seriestype = :shape, lw = 0.5,
      c = :blue, linecolor = :black, legend = false, fillalpha = 0.2, aspect_ratio=1)
plot!(p, circleShape(a.pos[1]+s.extent[1], a.pos[2]+s.extent[2], r), seriestype = :shape, lw = 0.5,
      c = :blue, linecolor = :black, legend = false, fillalpha = 0.2, aspect_ratio=1)
plot!(p, circleShape(a.pos[1], a.pos[2]+s.extent[2], r), seriestype = :shape, lw = 0.5,
      c = :blue, linecolor = :black, legend = false, fillalpha = 0.2, aspect_ratio=1)
# Truncate view
xlims!(0, s.extent[1])
ylims!(0, s.extent[2])
title!(p, "inexact")

# %% --- EXACT ---
function act(a)
    if a.id == agentid
        return :red
    elseif a.id in nearby_ids(model[agentid], model, r, exact=true)
        return :blue
    else
        return :yellow
    end
end
p = plotabm(model, as=4, ac=act, grid = (:both, :dot, 1, 0.9), xticks=(0:s.spacing:s.extent[1]), yticks=(0:s.spacing:s.extent[2]), size=(1000, 1000))

focal_cell = Agents.pos2cell(a.pos, model)
grid_r_max = r < model.space.spacing ? 1.0 : r / model.space.spacing + 1
allcells = Agents.grid_space_neighborhood(CartesianIndex(focal_cell), model, grid_r_max)
search_region = [(a .* s.spacing) .- (s.spacing/2) for a in allcells]
scatter!(p, search_region; markershape=:square, markersize=8, markerstrokewidth = 0, markeralpha = 0.2, markercolor=:grey)

grid_r_certain = grid_r_max - 1.2*sqrt(2)
certain_cells = Agents.grid_space_neighborhood(CartesianIndex(focal_cell), model, grid_r_certain)
search_region = [(a .* s.spacing) .- (s.spacing/2) for a in certain_cells]
scatter!(p, search_region; markershape=:diamond, markersize=8, markerstrokewidth = 0, markeralpha = 0.2, markercolor=:green)

plot!(p, circleShape(a.pos[1], a.pos[2], r), seriestype = :shape, lw = 0.5,
      c = :blue, linecolor = :black, legend = false, fillalpha = 0.2, aspect_ratio=1)
# Add all possible mirrors
plot!(p, circleShape(a.pos[1]-s.extent[1], a.pos[2], r), seriestype = :shape, lw = 0.5,
      c = :blue, linecolor = :black, legend = false, fillalpha = 0.2, aspect_ratio=1)
plot!(p, circleShape(a.pos[1]-s.extent[1], a.pos[2]-s.extent[2], r), seriestype = :shape, lw = 0.5,
      c = :blue, linecolor = :black, legend = false, fillalpha = 0.2, aspect_ratio=1)
plot!(p, circleShape(a.pos[1], a.pos[2]-s.extent[2], r), seriestype = :shape, lw = 0.5,
      c = :blue, linecolor = :black, legend = false, fillalpha = 0.2, aspect_ratio=1)
plot!(p, circleShape(a.pos[1]+s.extent[1], a.pos[2], r), seriestype = :shape, lw = 0.5,
      c = :blue, linecolor = :black, legend = false, fillalpha = 0.2, aspect_ratio=1)
plot!(p, circleShape(a.pos[1]+s.extent[1], a.pos[2]+s.extent[2], r), seriestype = :shape, lw = 0.5,
      c = :blue, linecolor = :black, legend = false, fillalpha = 0.2, aspect_ratio=1)
plot!(p, circleShape(a.pos[1], a.pos[2]+s.extent[2], r), seriestype = :shape, lw = 0.5,
      c = :blue, linecolor = :black, legend = false, fillalpha = 0.2, aspect_ratio=1)
# Truncate view
xlims!(0, s.extent[1])
ylims!(0, s.extent[2])
title!(p, "exact")

# 3D

mutable struct Ag3<: AbstractAgent
    id::Int
    pos::NTuple{3,Float64}
end
ε = 1.0
s = ContinuousSpace((10, 10, 10), ε)
model = ABM(Ag3, s)
for i in 1:800
    add_agent!(model)
end

function Sphere(O, r)
N = 32
u = LinRange(0, 2π, N)
v = LinRange(0, π, N)
x = cos.(u) * sin.(v)'
y = sin.(u) * sin.(v)'
z = repeat(cos.(v)',outer=[N, 1])
(x.*r.+O[1],y.*r.+O[2],z.*r.+O[3])
end

r = 2
agentid = rand(1:800)
a = model[agentid]
δ = Agents.distance_from_cell_center(a.pos, model)
grid_r = (r+δ)/model.space.spacing
focal_cell = Agents.pos2cell(a.pos, model)
allcells = Agents.grid_space_neighborhood(CartesianIndex(focal_cell), model, grid_r)
search_region = [(a .* s.spacing) .- (s.spacing/2) for a in allcells]

p = plotabm(model, as=4, ac=act,
            grid = (:both, :dot, 1, 0.9), xticks=(0:s.spacing:s.extent[1]), yticks=(0:s.spacing:s.extent[2]), zticks=(0:s.spacing:s.extent[3]),
            size=(1000, 1000, 1000))
#scatter!(search_region; markershape=:square, markersize=8, markerstrokewidth = 0, markeralpha = 0.2, markercolor=:grey)

plot!(p, Sphere(a.pos, r), seriestype = :shape, lw = 0.5,
      c = :blue, linecolor = :black, legend = false, fillalpha = 0.2, aspect_ratio=1)
# Truncate view
xlims!(0, s.extent[1])
ylims!(0, s.extent[2])
zlims!(0, s.extent[3])

plot!(p, camera=(20,50))

plot!(p, camera=(0,90))

plot!(p, camera=(90,0))

plot!(p, camera=(50,50))


text:
module TestProgressMacro

using ProgressLogging
using Test

let i = 0, x
    x = @progress for _ = 1:100
        i += 1
    end
    @test i == 100
    @test x == nothing
end

let i = 0, r = -50:10:50, x
    x = @progress for _ in r
        i += 1
    end
    @test i == 11
    @test x == nothing
end

let i = 0, x
    x = @progress "named" for _ = 1:100
        i += 1
    end
    @test i == 100
    @test x == nothing
end

let i = 0, j = 0, x
    x = @progress for _ = 1:10, __ = 1:20
        i += 1
    end
    @test i == 200
    @test x == nothing
end

let i = 0, j = 0, x
    bar = "bar"
    x = @progress "foo $bar" for _ = 1:10
        i += 1
    end
    @test i == 10
    @test x == nothing
end

let x, y
    x = @progress y = [i + 3j for i = 1:3, j = 1:4]
    @test y == reshape(4:15, 3, 4)
    @test x == y
end

let a = [], x
    x = @progress for i = 1:3, j in [-5, -2, -1, 8]
        j > 0 && continue
        push!(a, (i, j))
        i > 1 && break
    end
    @test a == [(1, -5), (1, -2), (1, -1), (2, -5)]
    @test x == nothing
end

end  # module


text:
module TopMetaTables
using  ...Dispatch.Structure
using  ...Dispatch.TableManipulation
export TopMetaTable, MetaTableNotFoundError, init_module_table!, init_metatable!, get_metatable,
       import_metatable!, export_metatable!, metatable_exports

struct TopMetaTable
  values::ObjectIdDict
  TopMetaTable() = new(ObjectIdDict())
end

struct MetaModuleTable
  values  :: Dict{Symbol, MetaMethodTable}
  exports :: Set{Symbol}
  MetaModuleTable() = new(Dict{Symbol, MetaMethodTable}(), Set{Symbol}())
end

Base.getindex(A::TopMetaTable,    mod::Module) = A.values[mod]
Base.getindex(A::MetaModuleTable, sym::Symbol) = A.values[sym]

Base.setindex!(A::TopMetaTable,    val, mod::Module) = A.values[mod] = val
Base.setindex!(A::MetaModuleTable, val, sym::Symbol) = A.values[sym] = val

Base.haskey(A::TopMetaTable,    mod::Module) = haskey(A.values, mod)
Base.haskey(A::MetaModuleTable, sym::Symbol) = haskey(A.values, sym)

struct MetaTableNotFoundError <: Exception
  key::Symbol
  mod::Module
end

Base.showerror(io::IO, err::MetaTableNotFoundError) =
  println(STDERR, ("Unable to find metatable $(err.key) in $(err.mod)."))

function init_module_table!(toptable, mod)
  if !haskey(toptable, mod)
    toptable[mod] = MetaModuleTable()
  end
end

function init_metatable!(toptable, mod, key, name)
  init_module_table!(toptable, mod)
  if !haskey(toptable.values[mod], key)
    toptable[mod][key] = MetaMethodTable(name)
  end
end

function import_metatable!(toptable, tablekey, from, to=current_module())
  init_module_table!(toptable, from)
  init_module_table!(toptable, to)
  !haskey(toptable[from], tablekey) && return warn("Unable to find metatable $(tablekey) in $(from).")
   haskey(toptable[to],   tablekey) && return warn("Ignoring conflicting import of $(tablekey) in $(to).")

  toptable[to][tablekey] = toptable[from][tablekey]
end

function export_metatable!(toptable, tablekey, from=current_module())
  init_module_table!(toptable, from)
  push!(toptable[from].exports, tablekey)
end

function metatable_exports(toptable, mod::Module)
  init_module_table!(toptable, mod)
  toptable[mod].exports
end

function get_metatable(toptable, mod, key)
  init_module_table!(toptable, mod)
  !haskey(toptable[mod], key) && throw(MetaTableNotFoundError(key,mod))

  toptable.values[mod][key]
end

end


text:
const tradestats = [foo, bar]
const rets = [.001, .020, -.231, .051]


text:
using DataFrames
using LinearAlgebra
using Statistics: std, var, mean, cov
using HighFrequencyCovariance
using Random
using Test

brownian_corr_matrix = Hermitian([1.0 0.75 0.5 0.0;
                                  0.0 1.0 0.5 0.25;
                                  0.0 0.0 1.0 0.25;
                                  0.0 0.0 0.0 1.0])
assets = [:BARC, :HSBC, :VODL, :RYAL]
twister = MersenneTwister(1)
ts, true_covar, micro_noise, update_rates = generate_random_path(4, 200000;
          brownian_corr_matrix = brownian_corr_matrix, assets = assets, vols = [0.02,0.03,0.04,0.05],
          twister = deepcopy(twister), syncronous = true)

preav_estimate1 = preaveraged_covariance(ts, assets)
covar = covariance(preav_estimate1,60)

spacing = 5
time_grid = collect(minimum(ts.df[:,ts.time]):spacing:maximum(ts.df[:,ts.time]))

eigenvalues, eigenvectors = eigen(covar)

dd = latest_value(ts, time_grid; assets = assets)

function make_factor(dd, cols, weights)
    return Matrix(dd[:,cols]) * weights
end
dd[!,:fac1] = make_factor(dd, assets, eigenvectors[:,1])
dd[!,:fac2] = make_factor(dd, assets, eigenvectors[:,2])
dd[!,:fac3] = make_factor(dd, assets, eigenvectors[:,3])
dd[!,:fac4] = make_factor(dd, assets, eigenvectors[:,4])

#using DataFrames
#using CSV
#using Plots
#using Lathe
using GLM
#using Statistics
#using StatsPlots
#using MLBase

fm = @formula(BARC  ghvy8y7cu  ~ Adult_Mortality)
linearRegressor = lm(fm, dd)

eigenvalues, eigenvectors = eigen(covar)
reconstituted = construct_matrix_from_eigen(eigenvalues, eigenvectors)

# factors


text:
using HypothesisTests, Test
using HypothesisTests: default_tail

@testset "Mann-Whitney" begin
@testset "Basic exact test" begin
	@test abs(pvalue(ExactMannWhitneyUTest([1:10;], [2.1:2:21;])) - 0.0232) <= 1e-4
	@test abs(pvalue(ExactMannWhitneyUTest([2.1:2:21;], [1:10;])) - 0.0232) <= 1e-4
	@test abs(pvalue(ExactMannWhitneyUTest([1.5:10:100;], [2.1:2:21;])) - 0.0068) <= 1e-4
	@test abs(pvalue(ExactMannWhitneyUTest([2.1:2:21;], [1.5:10:100;])) - 0.0068) <= 1e-4
	@test default_tail(ExactMannWhitneyUTest([1:10;], [2.1:2:21;])) == :both
	show(IOBuffer(), ExactMannWhitneyUTest([1:10;], [2.1:2:21;]))
end

@testset "Exact with ties" begin
	@test abs(pvalue(ExactMannWhitneyUTest([1:10;], [1:10;])) - 1) <= 1e-4
	@test abs(pvalue(ExactMannWhitneyUTest([1:10;], [2:11;])) - 0.5096) <= 1e-4
	@test abs(pvalue(ExactMannWhitneyUTest([2:11;], [1:10;])) - 0.5096) <= 1e-4
	@test abs(pvalue(ExactMannWhitneyUTest([1:10;], [1:5; ones(5)])) - 0.0057) <= 1e-4
	@test abs(pvalue(ExactMannWhitneyUTest([1:5; ones(5)], [1:10;])) - 0.0057) <= 1e-4
	show(IOBuffer(), ExactMannWhitneyUTest([1:10;], [1:10;]))
end

@testset "Exact with ties and unequal lengths" begin
	@test abs(pvalue(ExactMannWhitneyUTest([1:10;], [2:2:24;])) - 0.0118) <= 1e-4
	@test abs(pvalue(ExactMannWhitneyUTest([2:2:24;], [1:10;])) - 0.0118) <= 1e-4
	show(IOBuffer(), ExactMannWhitneyUTest([1:10;], [2:2:24;]))
end

@testset "Approximate test" begin
	@test abs(pvalue(ApproximateMannWhitneyUTest([1:10;], [2.1:2:21;])) - 0.0257) <= 1e-4
	@test abs(pvalue(ApproximateMannWhitneyUTest([2.1:2:21;], [1:10;])) - 0.0257) <= 1e-4
	@test abs(pvalue(ApproximateMannWhitneyUTest([1.5:10:100;], [2.1:2:21;])) - 0.0091) <= 1e-4
	@test abs(pvalue(ApproximateMannWhitneyUTest([2.1:2:21;], [1.5:10:100;])) - 0.0091) <= 1e-4
	@test default_tail(ApproximateMannWhitneyUTest([1:10;], [2.1:2:21;])) == :both
	show(IOBuffer(), ApproximateMannWhitneyUTest([1:10;], [2.1:2:21;]))
end

@testset "Approximate with ties" begin
	@test abs(pvalue(ApproximateMannWhitneyUTest([1:10;], [1:10;])) - 1) <= 1e-4
	@test abs(pvalue(ApproximateMannWhitneyUTest([1:10;], [2:11;])) - 0.4948) <= 1e-4
	@test abs(pvalue(ApproximateMannWhitneyUTest([2:11;], [1:10;])) - 0.4948) <= 1e-4
	@test abs(pvalue(ApproximateMannWhitneyUTest([1:10;], [1:5; ones(5)])) - 0.0076) <= 1e-4
	@test abs(pvalue(ApproximateMannWhitneyUTest([1:5; ones(5)], [1:10;])) - 0.0076) <= 1e-4
	show(IOBuffer(), ApproximateMannWhitneyUTest([1:10;], [1:10;]))
end

@testset "Tests for automatic selection" begin
	@test abs(pvalue(MannWhitneyUTest([1:10;], [2.1:2:21;])) - 0.0232) <= 1e-4
	@test abs(pvalue(MannWhitneyUTest([1:10;], [2:11;])) - 0.4948) <= 1e-4
	show(IOBuffer(), MannWhitneyUTest([1:10;], [2.1:2:21;]))
end

@testset "Issue #39" begin
	@test abs(pvalue(ExactMannWhitneyUTest(Float32[1:10;], Float32[2:11;])) - 0.5096) <= 1e-4
end
end


text:
using IfsFractal

@everywhere function draworbits(rng::MersenneTwister,pchunk_size::Int)
	canvas=IfsFractal.FractalCanvas(512,512,Complex(-0.5,0.0),0.25)
	for j=1:pchunk_size
		c=IfsFractal.randcomplex(rng)
		orbit=IfsFractal.mandelbrot(c)
		if orbit.escaped && length(orbit.zs) > 1
			IfsFractal.draw!(canvas,orbit)
		end
	end
	return canvas
end

function main2()
	rng=MersenneTwister(1)
	const num_orbits=100000
	const chunk_size=100000
	np=nprocs()
	canvases=Array{FractalCanvas,1}()
	# for p=1:np
	# 	push!(canvases,FractalCanvas(512,512,zero(Complex),0.25))
	# end
	rngs=randjump(rng,np)
	final_canvas=FractalCanvas(512,512,Complex(-0.5,0.0),0.25)

	for i=1:num_orbits/chunk_size
		# @time cs=[randcomplex(rng,2.0) for j=1:chunk_size]
		# start_cs_idx=1
		futures=Array{RemoteRef}(np)
		@time begin
		for p=1:np
			proc_chunk_size=trunc(Int,chunk_size/np)
			if p==np
				proc_chunk_size=chunk_size-((np-1)*proc_chunk_size)
			end
			futures[p]=@spawn draworbits(rngs[p],proc_chunk_size)
			# start_cs_idx+=proc_chunk_size
		end
		println("Spawned.")
		end
		@time begin
		for p=1:np
			push!(canvases,fetch(futures[p]))
			@show minimum(canvases[p].data),maximum(canvases[p].data),mean(canvases[p].data)
		end
		end
		# proc_chunk_size=trunc(Int,chunk_size/nprocs())
		# @time cs=[randcomplex(rng,2.0) for j=1:chunk_size]
		# future
		s=Array{RemoteRef}(nprocs())
		# @time begin
		# @everywhere function genorbits(pcs::Array{Complex{Float64}})
		# 	orbits=Array{IfsFractal.Orbit}(length(pcs))
		# 	for j=1:length(pcs)
		# 		orbits[j]=IfsFractal.mandelbrot(pcs[j])
		# 	end
		# 	orbits
		# end
		# start_cs_idx=1
		# for p=1:nprocs()
		# 	if p==nprocs()
		# 		proc_chunk_size=chunk_size-((nprocs()-1)*proc_chunk_size)
		# 	end
		# 	# @show typeof(cs[start_cs_idx:start_cs_idx+proc_chunk_size-1])
		# 	futures[p]=@spawnat p genorbits(cs[start_cs_idx:start_cs_idx+proc_chunk_size-1])
		# 	start_cs_idx+=proc_chunk_size
		# end
		# orbits=Array{Orbit,1}()
		# for p=1:nprocs()
		# 	porbits=fetch(futures[p])
		# 	# @show porbits
		# 	append!(orbits,porbits)
		# end
		# end
		# @show size(orbits)

		# orbits=Array{Orbit}(chunk_size)
		# orbit_ids=SharedArray(Int,chunk_size)
		# for j=1:chunk_size
		# 	orbit_ids[j]=j
		# end
		# @time cs=[randcomplex(rng,2.0) for j=1:chunk_size]
		# @time orbits=pmap(mandelbrot,cs)
		# @time begin
		# orbits=@sync @parallel (vcat) for j=1:chunk_size
		# 	mandelbrot(cs[j])
		# end
		# end
		# @time begin
		# for j=1:chunk_size
		# 	c=randcomplex(rng,2.0)
		# 	orbits[j]=mandelbrot(c,50000)
		# end
		# end
		# # @show orbits
		# @time begin
		# for j=1:chunk_size
		# 	if orbits[j].escaped
		# 		draw!(canvas,orbits[j])
		# 	end
		# end
		# end
		println("$(trunc(Int,i*chunk_size))/$num_orbits")
		for p=1:np
			# @show maximum(canvases[p].data),mean(canvases[p].data)
			final_canvas.data[:,:]+=canvases[p].data[:,:]
		end
	end
	# reduce
	f=open("canvas.dat","w")
	serialize(f,final_canvas)
	close(f)
	render(final_canvas,"out.4.png.bias0.05.png",0.05)
	render(final_canvas,"out.4.png.bias0.1.png",0.1)
	render(final_canvas,"out.4.png.bias0.2.png",0.2)
	render(final_canvas,"out.4.png.bias0.3.png",0.3)
end
main2()

using Colors
function main3()
	f=open("canvas.dat","r")
	canvas=deserialize(f)
	close(f)
	black=RGB{Float64}(0.0,0.0,0.0)
	white=RGB{Float64}(1.0,1.0,1.0)
	render(canvas,"out.4.png.bias0.05.png",0.05,black,white)
	render(canvas,"out.4.png.bias0.1.png",0.1,black,white)
	render(canvas,"out.4.png.bias0.2.png",0.2,black,white)
	render(canvas,"out.4.png.bias0.3.png",0.3,black,white)
end
# main3()

max_stars_repo_name: jblondin/Fractal

max_issues_repo_name: jblondin/Fractal

max_forks_repo_name: jblondin/Fractal


text:
"""
# module Phonon

# Examples

```jldoctest
julia>
```
"""
module Phonon

include("Matdyn.jl")

end


text:
function irisws(cha::String,
                 d0::String,
                 d1::String,
                fmt::String,
               opts::String,
                 to::Int64,
                  v::Integer,
                  w::Bool)

  # init
  parse_err = false
  parsable = false
  fname = ""

  # parse channel string cha
  c = (parse_chstr(cha, ',', false, false)[1,:])[1:min(end,4)]
  if isempty(c[3])
    c[3] = "--"
  end
  ID = join([c[1], c[2], strip(c[3],'-'), c[4]], '.')

  # Build query url
  url = "http://service.iris.edu/irisws/timeseries/1/query?" *
          build_stream_query(c, d0, d1) * "&scale=" * (fmt == "miniseed" ? "AUTO" : "1.0") * "&format=" * fmt
  v > 0 && println(url)
  req_info_str = datareq_summ("IRISWS data", ID, d0, d1)
  # see CHANGELOG, 2020-05-28

  # Do request
  (R, parsable) = get_http_req(url, req_info_str, to)
  if parsable
    if w
      fname = savereq(R, fmt, ID, d0)
    end

    if fmt == "sacbl"
      Ch = read_sac_stream(IOBuffer(R), false, false)
    elseif fmt == "miniseed"
      S = SeisData()
      parsemseed!(S, IOBuffer(R), KW.nx_add, KW.nx_add, true, v)
      Ch = S[1]
    elseif fmt == "geocsv"
      S = SeisData()
      read_geocsv_tspair!(S, IOBuffer(R))
      Ch = S[1]
    else
      # other parsers not supported
      parse_err = true
      Ch = SeisChannel(id = string("XX.FMT..001"),
                       misc = Dict{String,Any}(
                         "url" => url,
                         "raw" => read(IOBuffer(R))
                       )
                     )
      note!(Ch, "unparseable format; raw bytes in :misc[\"raw\"]")
    end

  else
    parse_err = true
    Ch = SeisChannel(id = string("XX.FAIL..001"),
                     misc = Dict{String,Any}(
                       "url" => url,
                       "msg" => String(read(IOBuffer(R)))
                       )
                     )
    note!(Ch, "request failed; response in :misc[\"msg\"]")
  end
  setfield!(Ch, :src, url)
  note!(Ch, "+source ¦ " * url)

  # fill :id and empty fields if no parse_err
  if parse_err == false
    setfield!(Ch, :id, ID)
    if isempty(Ch.name)
      Ch.name = deepcopy(ID)
    end
    unscale!(Ch)        # see CHANGELOG, 2020-05-28
    Ch.loc = GeoLoc()   # see CHANGELOG, 2020-05-28
  end
  if parsable && w
    push!(Ch.notes, string(timestamp(), " ¦ write ¦ get_data(\"IRIS\", ... w=true) ¦ wrote raw download to file ", fname))
  end
  return parse_err, Ch
end

# Programming note: if this method is the default, and S is only modified
# within the for loop, then S is copied to a local scope and the newly-added
# data are deleted upon return
function IRISget!(S::GphysData,
                  C::Array{String, 1},
                 d0::String,
                 d1::String,
                fmt::String,
               opts::String,
                 to::Int64,
                  v::Integer,
                  w::Bool)

  parse_err = false
  if fmt == "mseed"
    fmt = "miniseed"
  elseif fmt == "sac"
    fmt = "sacbl"
  end
  U = SeisData()
  K = size(C, 1)
  v > 0 && println("IRISWS data request begins...")
  for k = 1:K
    (p, Ch) = irisws(C[k], d0, d1, fmt, opts, to, v, w)
    push!(U, Ch)
    parse_err = max(parse_err, p)
  end
  append!(S, U)
  return parse_err
end


text:
using EventTimeSeries
using Test

@testset "test neighbor iterator" begin
    @test EventTimeSeries.neighbors([1,2,3]) |> collect == [(1, 2), (2, 3)]
end

@testset "test splice of timeseries with unique and equal tags" begin
    t1 = [1,3]
    v1 = [:a,:c]
    t2 = [2,4]
    v2 = [:b,:d]
    tag = :foo

    ts1 = EventTS(timestamps=t1, values=v1, tag=tag)
    ts2 = EventTS(timestamps=t2, values=v2, tag=tag)

    for (a,b) in ((ts1,ts2), (ts2,ts1))
        ts = splice(a, b)
        @test timestamps(ts) |> collect == [1,2,3,4]
        @test values(ts) |> collect == [:a,:b,:c,:d]
        @test ts.tag == tag
    end
end

@testset "test splice of timeseries with unique but unequal tags" begin
    t1 = [1,3]
    v1 = [:a,:c]
    t2 = [2,4]
    v2 = [:b,:d]
    tag1 = :foo
    tag2 = :bar

    ts1 = EventTS(timestamps=t1, values=v1, tag=tag1)
    ts2 = EventTS(timestamps=t2, values=v2, tag=tag2)

    for (a,b) in ((ts1,ts2), (ts2,ts1))
        ts = splice(a, b)
        @test timestamps(ts) |> collect == [1, 2, 3, 4]
        @test values(ts) |> collect == [:a, :b, :c, :d]
        @test ts.tag == [:foo, :bar, :foo, :bar]
    end
end

@testset "test splice of timeseries with both varying tags" begin
    t1 = [1,3]
    v1 = [:a,:c]
    t2 = [2,4]
    v2 = [:b,:d]
    tag1 = [:foo, :bar]
    tag2 = [:foo, :bar]

    ts1 = EventTS(timestamps=t1, values=v1, tag=tag1)
    ts2 = EventTS(timestamps=t2, values=v2, tag=tag2)

    for (a,b) in ((ts1,ts2), (ts2,ts1))
        ts = splice(a, b)
        @test timestamps(ts) |> collect == [1,2,3,4]
        @test values(ts) |> collect == [:a,:b,:c,:d]
        @test ts.tag == [:foo,:foo,:bar,:bar]
    end
end

@testset "test splice of timeseries with unique and varying tag" begin
    t1 = [1,3]
    v1 = [:a,:c]
    t2 = [2,4]
    v2 = [:b,:d]
    tag1 = [:foo, :bar]
    tag2 = :foobar

    ts1 = EventTS(timestamps=t1, values=v1, tag=tag1)
    ts2 = EventTS(timestamps=t2, values=v2, tag=tag2)

    for (a,b) in ((ts1,ts2), (ts2,ts1))
        ts = splice(a, b)
        @test timestamps(ts) |> collect == [1,2,3,4]
        @test values(ts) |> collect == [:a,:b,:c,:d]
        @test ts.tag == [:foo,:foobar,:bar,:foobar]
    end
end

@testset "splice -> split gives back result" begin
    t1 = [1,3]
    v1 = [:a,:c]
    tag1 = :A
    ts1 = EventTS(timestamps=t1, values=v1, tag=tag1)

    t2 = [2,4]
    v2 = [:b,:d]
    tag2 = :B
    ts2 = EventTS(timestamps=t2, values=v2, tag=tag2)

    ts1_, ts2_ = split(splice(ts1, ts2))

    @test ts1.timestamps == ts1_.timestamps
    @test ts1.tag == ts1_.tag
    @test ts1.values == ts1_.values
end

@testset "split -> splice gives back result" begin
    t = [1,2,3,4]
    v = rand(4)
    tag = [:a,:b,:a,:c]
    ts = EventTS(timestamps=t, values=v, tag=tag)
    ts_ = splice(split(ts)...)

    @test ts_.timestamps == ts.timestamps
    @test ts_.tag == ts.tag
    @test ts_.values == ts.values
end

@testset "drop_repeated 1" begin
    t1 = 0:2:10 |> collect
    t2 = 1:2:9 |>  collect
    tag = :tagtag
    val1 = 1:length(t1)
    val2 = 1:length(t2)

    ts1 = EventTS(timestamps=t1, values=val1, tag=tag)
    ts2 = EventTS(timestamps=t2, values=val2, tag=tag)
    tsd = drop_repeated(splice(ts1,ts2))

    @test ts1.timestamps == tsd.timestamps
    @test ts1.tag == tsd.tag
    @test ts1.values == tsd.values
end

@testset "drop_repeated 2" begin
    time = 1:6 |> collect
    tag = [i==4 ? :foo : :bar for (i,t) in enumerate(time)]
    values = [:a, :b, :b, :b, :d, :d]
    ts = EventTS(timestamps=time, tag=tag, values=values)
    ts1 = drop_repeated(ts) # keep_end true by default
    ts2 = drop_repeated(ts, keep_end=false)

    @test ts1.timestamps == [1,2,4,5,6]
    @test ts1.values == [:a, :b, :b, :d, :d]

    @test ts2.timestamps == [1,2,4,5]
    @test ts2.values == [:a, :b, :b, :d]
end

@testset "merge_tags fill forward" begin
    time   = [0, 1, 2, 3, 4, 5]
    tag    = [:A, :A,:B,:A,:B, :B]
    values = time

    ts  = EventTS(timestamps=time, values=values, tag=tag)
    ts_ = merge_tags(ts)

    @test ts_.values == [(0,nothing),
                         (1,nothing),
                         (1,2),
                         (3,2),
                         (3,4),
                         (3,5)]

    @test ts_.tag == (:A, :B)
end


text:
"""
    make_hermitian(arr)

Takes an array `arr` and appends rows, cols, ... if necessary
so that `arr` is a hermitian array which preserves Parseval's theorem.

# Examples
```jldoctest
julia> FFTResampling.make_hermitian([1.0 2.0])
1×3 Array{Float64,2}:
 0.5  2.0  0.5

julia> FFTResampling.make_hermitian([1.0 2.0; 3.0 4.0])
3×3 Array{Float64,2}:
 0.5  1.0  0.0
 1.5  4.0  1.5
 0.0  1.0  0.5

julia> FFTResampling.make_hermitian([1im 2.0; 3.0im 4.0im; 5.0 6.0im])
3×3 Array{Complex{Float64},2}:
 0.0+0.5im  2.0+0.0im  2.5-0.0im
 0.0+1.5im  0.0+4.0im  0.0-1.5im
 2.5+0.0im  0.0+6.0im  0.0-0.5im
```
"""
function make_hermitian(arr::AbstractArray{T, N}) where {T, N}
    # if the size is odd, we need to add a slice to that dimension
    # at this slice we then fill new values, so that `arr` is hermitian
    mf(x) = x % 2 == 0 ? x +1 : x
    size_new = map(mf, size(arr))
    # new array with new size
    arr_new = zeros(eltype(arr), size_new)

    # fill the new array with the old one
    same_ind = map(x -> 1:x, size(arr))
    arr_new[same_ind...] = arr

    # copy again to have a out array (might be necessary because of reassigning)
    arr_out = copy(arr_new)

    # now modify values at the edges to get hermitian property
    # and Parsveval's theorem correct
    @inbounds for d = 1:N
        # special case when singleton dimension
        if size(arr_new)[d] == 1 || size(arr)[d] % 2 == 1
            continue
        end

        # extract indices represent the left slice we need to copy
        extract_indices = slice_indices(arr_new, d, 1)
        y = arr_new[extract_indices...] ./ 2
        y_rev = reverse_all(y)

        # assign indices represent the right sliuce we need to assign the
        # conjugated and halved one
        assign_indices = slice_indices(arr_new, d, size(arr_out)[d])

        # assign the halved
        arr_out[extract_indices...] = y
        # assign the conjugated, halved one
        arr_out[assign_indices...] = conj.(y_rev)
    end

    return arr_out
end

"""
    get_indices_around_center(i_in, i_out)
A function which provides two output indices `i1` and `i2`
where `i2 - i1 = i_out`
The indices are chosen in a way that the set `i1:i2`
cuts the interval `1:i_in` in a way that the center frequency
stays at the center position.
Works for both odd and even indices
"""
function get_indices_around_center(i_in, i_out)
    if (mod(i_in, 2) == 0 && mod(i_out, 2) == 0
     || mod(i_in, 2) == 1 && mod(i_out, 2) == 1)
        x = (i_in - i_out) ÷ 2
        return 1 + x, i_in - x
    elseif mod(i_in, 2) == 1 && mod(i_out, 2) == 0
        x = (i_in - 1 - i_out) ÷ 2
        return 1 + x, i_in - x - 1
    elseif mod(i_in, 2) == 0 && mod(i_out, 2) == 1
        x = (i_in - (i_out - 1)) ÷ 2
        return 1 + x, i_in - (x - 1)
    end
end

"""
    center_extract(arr, new_size_array)
Extracts a center of an array.
`new_size_array` must be list of sizes indicating the output
size of each dimension. Centered means that a center frequency
stays at the center position. Works for even and uneven.
If `length(new_size_array) < length(ndims(arr))` the remaining dimensions
are untouched and copied.
# Examples
```jldoctest
julia> FFTResampling.center_extract([1 2; 3 4], [1])
1×2 Array{Int64,2}:
 3  4

julia> FFTResampling.center_extract([1 2; 3 4], [1, 1])
1×1 Array{Int64,2}:
 4

julia> FFTResampling.center_extract([1 2 3; 3 4 5; 6 7 8], [2 2])
2×2 Array{Int64,2}:
 1  2
 3  4
```
"""
function center_extract(arr::AbstractArray, new_size_array)
    new_size_array = collect(new_size_array)

    # we construct two lists
    # the reason is, that we don't change higher dimensions which are not
    # specified in new_size_array
    out_indices1 = [get_indices_around_center(size(arr)[x], new_size_array[x])
                    for x = 1:length(new_size_array)]

    out_indices1 = [x[1]:x[2] for x = out_indices1]

    # out_indices2 contains just ranges covering the full size of each dimension
    out_indices2 = [1:size(arr)[i] for i = (1 + length(new_size_array)):ndims(arr)]
    return arr[out_indices1..., out_indices2...]
end

"""
    center_set!(arr_large, arr_small)
Puts the `arr_small` central into `arr_large`.
The convention, where the center is, is the same as the definition
as for FFT based centered.
Function works both for even and uneven arrays.
# Examples
```jldoctest
julia> FFTResampling.center_set!([1, 1, 1, 1, 1, 1], [5, 5, 5])
6-element Array{Int64,1}:
 1
 1
 5
 5
 5
 1
```
"""
function center_set!(arr_large, arr_small)
    out_is = []
    for i = 1:ndims(arr_large)
        a, b = get_indices_around_center(size(arr_large)[i], size(arr_small)[i])
        push!(out_is, a:b)
    end

    #rest = ones(Int, ndims(arr_large) - 3)
    arr_large[out_is...] = arr_small

    return arr_large
end

"""
    center_pos(x)
Calculate the position of the center frequency.
Size of the array is `x`
# Examples
```jldoctest
julia> FFTResampling.center_pos(3)
2
julia> FFTResampling.center_pos(4)
3
```
"""
function center_pos(x::Integer)
    # integer division
    return div(x, 2) + 1
end

"""
    slice(arr, dim, index)

Return a `N` dimensional slice (where one dimensions has size 1) of the N-dimensional `arr` at the index position
`index` in the `dim` dimension of the array.
It holds `size(out)[dim] == 1`.

# Examples
```jldoctest
julia> x = [1 2 3; 4 5 6; 7 8 9]
3×3 Array{Int64,2}:
 1  2  3
 4  5  6
 7  8  9

julia> FFTResampling.slice(x, 1, 1)
1×3 view(::Array{Int64,2}, 1:1, :) with eltype Int64:
 1  2  3

julia> FFTResampling.slice(x, 2, 3)
3×1 view(::Array{Int64,2}, :, 3:3) with eltype Int64:
 3
 6
 9

```

"""
function slice(arr::AbstractArray{T, N}, dim::Integer, index::Integer) where {T, N}
    inds = slice_indices(arr, dim, index)
    return @view arr[inds...]
end

function slice_indices(arr::AbstractArray{T, N}, dim::Integer, index::Integer) where {T, N}
    f(x) = x[1] == dim ? (index:index) : (1:x[2])
    inds = map(f, enumerate(size(arr)))
    return inds
end

"""
    reverse_all(arr)

Reverse an array `arr` over all dimensions.

# Examples
```jldoctest
julia> FFTResampling.reverse_all([1 2 3 4 5])
1×5 Array{Int64,2}:
 5  4  3  2  1

julia> FFTResampling.reverse_all([1 2 3 4 5; 7 8 9 10 11])
2×5 Array{Int64,2}:
 11  10  9  8  7
  5   4  3  2  1

julia> FFTResampling.reverse_all([1; 2; 3; 4; 5])
5-element Array{Int64,1}:
 5
 4
 3
 2
 1
```
"""
function reverse_all(arr::AbstractArray)
    return arr[reverse_all_indices(arr)...]
end

function reverse_all_indices(arr)
    out = []
    for i in size(arr)
        push!(out, i:-1:1)
    end
    return out
end

function dft_1D(arr)
    N = length(arr)
    X = zeros(ComplexF64, N)

    for k = 0:N-1
        s = zero(X[1])
        for n = 0:N-1
            s += arr[n + 1]  * exp(-1im*2*π*k *n/ N)
        end
        X[k+1] = s
    end
    return X
end

"""
This function ensures the constraints a spectrum must have so that after
iffting the result is purely real.

For example, inspect `ffshift(fft(randn((6,6))))` and you see that the highest negative
frequencies obey a certain symmetry. This can be generalized to N dimensions.
If we cut the array so that it'll be even in the that d dimension (new_size[d] % 2 ==0)
we need to manually restore the symmetry to guarantee a purely real result.
"""
function add_high_frequencies(size_arr, arr_out_f, new_size, N)
    fix_mul = false
    for d = 1:N
        if new_size[d] % 2 == 1 || new_size[d] == size_arr[d] ||
            (size_arr[d] % 2 == 0 && new_size[d] == (size_arr[d] + 1))
            continue
        end

        # construct the slice containing the highest positive frequency
        inds_extract = []
        inds_assign = []

        for i = 1:N
            a,b = get_indices_around_center(size_arr[i], new_size[i])
            if i == d
                # b+1 is the highest positive frequency which
                # will be cut off (under certain conditions)
                iend = min(b+1, size_arr[i])
                push!(inds_extract, iend:iend)
                push!(inds_assign, a:a)
            else
                push!(inds_extract, a:min(b+1, size_arr[i]))
                push!(inds_assign, a:min(b+1, size_arr[i]))
            end
        end

        # add the highest positive frequency slice to the highest negative
        # but use the same arr_out_f
        # in that way we create a matrix which has the demanded symmetry of FFT to produce a real result
        arr_out_f[inds_assign...] = (arr_out_f[inds_assign...] .+ arr_out_f[inds_extract...])
    end

    return arr_out_f
end


text:
## import packages
using MAT,Plots,Dates,TimerOutputs,WriteVTK,DataFrames,CSV,
Statistics,DelimitedFiles,ProgressMeter

const USE_GPU=false;  # Use GPU? If this is set false, then no GPU needs to be available
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
    @init_parallel_stencil(CUDA,Float64,3);
else
    @init_parallel_stencil(Threads,Float64,3);
end
## 8th-order derivatives
@parallel function Dx_inn(in,out)
    @inn(out)=@d_xi_8(in);
    return nothing
end

@parallel function Dy_inn(in,out)
    @inn(out)=@d_yi_8(in);
    return nothing
end

@parallel function Dz_inn(in,out)
    @inn(out)=@d_zi_8(in);
    return nothing
end

@parallel_indices (iy,iz) function u_1_plus(in,out)
out[4:end-2,iy,iz]=in[:,iy,iz];
return nothing
end

@parallel_indices (iy,iz) function u_1_minus(in,out)
out[3:end-3,iy,iz]=in[:,iy,iz];
return nothing
end

@parallel_indices (ix,iz) function u_2_plus(in,out)
out[ix,4:end-2,iz]=in[ix,:,iz];
return nothing
end

@parallel_indices (ix,iz) function u_2_minus(in,out)
out[ix,3:end-3,iz]=in[ix,:,iz];
return nothing
end

@parallel_indices (ix,iy) function u_3_plus(in,out)
out[ix,iy,4:end-2]=in[ix,iy,:];
return nothing
end

@parallel_indices (ix,iy) function u_3_minus(in,out)
out[ix,iy,3:end-3]=in[ix,iy,:];
return nothing
end
##
nx=12;
ny=12;
nz=12;
A=@zeros(nx,ny,nz);
dtt1=@zeros(nx-5,ny,nz);
dtt2=@zeros(nx,ny-5,nz);
dtt3=@zeros(nx,ny,nz-5);
C=@zeros(nx,ny,nz);
A[5,5,5]=1;
## Dx_plus: place derivative in the later point
@parallel Dx_inn(A,dtt1);
@parallel (2:ny-1,2:nz-1) u_1_plus(dtt1,C);
## Dx_minus: place derivative in the previous point
@parallel Dx_inn(A,dtt1);
@parallel (2:ny-1,2:nz-1) u_1_minus(dtt1,C);
## Dy_plus: place derivative in the later point
@parallel Dy_inn(A,dtt2);
@parallel (2:nx-1,2:nz-1) u_2_plus(dtt2,C);
## Dy_minus: place derivative in the previous point
@parallel Dy_inn(A,dtt2);
@parallel (2:nx-1,2:nz-1) u_2_minus(dtt2,C);
## Dz_plus: place derivative in the later point
@parallel Dz_inn(A,dtt3);
@parallel (2:nx-1,2:ny-1) u_3_plus(dtt3,C);
## Dz_minus: place derivative in the previous point
@parallel Dz_inn(A,dtt3);
@parallel (2:nx-1,2:nz-1) u_3_minus(dtt3,C);
## 2nd-order derivatives
@parallel function Dx_inn(in,out)
    @inn(out)=@d_xi_1(in);
    return nothing
end

@parallel_indices (iy,iz) function u_1_plus(in,out)
out[2:end-1,iy,iz]=in[:,iy,iz];
return nothing
end

@parallel_indices (iy,iz) function u_1_minus(in,out)
out[1:end-2,iy,iz]=in[:,iy,iz];
return nothing
end
##
nx=12;
ny=12;
nz=12;
IND_accuracy=2;
A=@zeros(nx,ny,nz);
dtt1=@zeros(nx-IND_accuracy,ny,nz);
dtt2=@zeros(nx,ny-IND_accuracy,nz);
dtt3=@zeros(nx,ny,nz-IND_accuracy);
C=@zeros(nx,ny,nz);
A[5,5,5]=1;
## Dx_plus: place derivative in the later point
@parallel Dx_inn(A,dtt1);
@parallel (2:ny-1,2:nz-1) u_1_plus(dtt1,C);
## Dx_minus: place derivative in the previous point
@parallel Dx_inn(A,dtt1);
@parallel (2:ny-1,2:nz-1) u_1_minus(dtt1,C);
## 12th-order derivatives
@parallel function Dx_inn(in,out)
    @inn(out)=@d_xi_12(in);
    return nothing
end

@parallel function Dy_inn(in,out)
    @inn(out)=@d_yi_12(in);
    return nothing
end

@parallel function Dz_inn(in,out)
    @inn(out)=@d_zi_12(in);
    return nothing
end

@parallel_indices (iy,iz) function u_1_plus(in,out)
out[6:end-4,iy,iz]=in[:,iy,iz];
return nothing
end

@parallel_indices (iy,iz) function u_1_minus(in,out)
out[5:end-5,iy,iz]=in[:,iy,iz];
return nothing
end

@parallel_indices (ix,iz) function u_2_plus(in,out)
out[ix,6:end-4,iz]=in[ix,:,iz];
return nothing
end

@parallel_indices (ix,iz) function u_2_minus(in,out)
out[ix,5:end-5,iz]=in[ix,:,iz];
return nothing
end

@parallel_indices (ix,iy) function u_3_plus(in,out)
out[ix,iy,6:end-4]=in[ix,iy,:];
return nothing
end

@parallel_indices (ix,iy) function u_3_minus(in,out)
out[ix,iy,5:end-5]=in[ix,iy,:];
return nothing
end
##
nx=21;
ny=21;
nz=21;
A=@zeros(nx,ny,nz);
dtt1=@zeros(nx-9,ny,nz);
dtt2=@zeros(nx,ny-9,nz);
dtt3=@zeros(nx,ny,nz-9);
C=@zeros(nx,ny,nz);
A[11,11,11]=1;
## Dx_plus: place derivative in the later point
@parallel Dx_inn(A,dtt1);
@parallel (2:ny-1,2:nz-1) u_1_plus(dtt1,C);
## Dx_minus: place derivative in the previous point
@parallel Dx_inn(A,dtt1);
@parallel (2:ny-1,2:nz-1) u_1_minus(dtt1,C);


text:
export renewal_counting
"""
```julia
renewal_counting(n::Int, t, pdf, first_pdf = pdf) -> totalt, Qn, Pn
```
Recursively calculates the p.d.f. of the `n`-th event, `Qn`, happening at time
`totalt`, as
well as the probability to have `n` events, `Pn`, occuring up to time `totalt`.
# Inputs:
1. `n::Int` : The number of events you care about. Notice that everything
  starting from `k=1` up to `k=n` will be calculated (and returned).
2. `t::Vector` : The time vector versus which the `pdf` is calculated. Has to have a
  `step()` method defined (i.e. be equally-spaced).
3. `pdf::Vector` : The vector containing the histogram (prob. density function)
  of the
  interarrival times of the events. It should be integral-normalized versus `t`,
  and if it is not, this is taken care internally. It is expected that the
  histogram/p.d.f. is calculated with **left-closed intervals**! It is also
  expected that the p.d.f. has finite support, i.e. it is 0 for any time ∉ `t`.
4. `first_pdf::Vector` : [optional] The p.d.f. of the interarrival times of the
  *first* event, in
  case it is different from the general p.d.f.
# Outputs:
1. `totalt::Vector` : Time vector of length: `L = l*(n+1) - n` where `l = length(t)`
2. `Qn::Matrix` : A `L×n` matrix. Each column `j` is the **p.d.f.** of the
  `j`-th event versus `totalt`,  for `j in 1:n`
3. `Pn::Matrix` : A `L×(n+1)` matrix. Each column `j` is the **probability** to
  have a total of `j` events up to time `totalt` for `j in 0:n`.

Notice that the output `Qn` represents probability density function, while the
output `Pn` represents commulative probability function (*direct* probability).
"""
function renewal_counting(n::Int, t::AbstractVector, pdf::AbstractVector,
  first_pdf = pdf)

  if length()
  l = length(t)
  tmax = maximum(t)
  integral = trapezint(t, pdf)

  if !isapprox(integral, 1)
    warn("Integral of p.d.f. versus t was not normalized,
using a re-normalized copy instead...")
    f = deepcopy(pdf)
    f ./= integral
  else
    f = pdf
  end
  try
    stepp = step(t)
  catch
    throw(ArgumentError("Input `t` must be equally spaced (have `step()`)"))
  end

  #this is the time vector my functions will be functions of
  totalt = 0:stepp:ceil((n+1)*tmax, -Int(mag(stepp)))
  L = length(totalt) #should be l*n+1 - n because you only count 0.0 once
  append!(f, zeros(L - l)) # extend f for totalt
  # Calculate average event time:
  κ = trapezint(stepp, totalt.*f)
  W = zeros(totalt)
  # Calculate general survival function W:
  for i in 1:l
    W[i] = 1 - trapezint(t[1:i], f[1:i])
  end

  # Define the matrix P: each column a P_n, from n=0
  # P_n[i] is probability to have n events up to time t[i]
  P = zeros(eltype(t), L, n+1)
  # P_0, probability to have 0 events up to time t (identical to
  # survival function of first event)
  W = zeros(totalt)
  # Calculate general survival function W:
  for i in 1:l
    P[i, 1] = 1 - trapezint(t[1:i], first_pdf[1:i])
  end
  # Define the matrix Q: each column a Q_n from n=1
  # Qn is the p.d.f. of the n-th event happening at time t
  Q = zeros(eltype(t), L, n)
  # Calculate Q1: probability that 1st event happens at time t
  # equivalent to p.d.f. of the first event or to general p.d.f.
  # if they are the same
  Q1[1:l] .= first_pdf
  # Calculate Qn, Pn
  for j in 1:n
    # Qn:
    if j != 1
      reverse_Qnm1 = reverse(Q[:,j-1])
      for t_idx in 2:L #no reason to start from 1, gives 0 anyway
        Q[t_idx, j] =
        trapezint(totalt[1:t_idx], f[1:t_idx].*reverse_Qnm1[(L-t_idx+1):end])
      end
      qint = trapezint(totalt, Q[:, j])
      if !isapprox(qint, 1; rtol = 1e-14)
        Q[:, j] ./= qint
      end
    end

    # Pn:
    reverse_Qn = reverse(Q[:,j])
    for t_idx in 2:L #t_idx is the timepoint t I want to find the Pn(t)
      P[t_idx, j+1] =
      trapezint(totalt[1:t_idx], W[1:t_idx].*reverse_Qn[(L-t_idx+1):end])
    end
  end
  return totalt, Q, P
end


text:
function parse_line(l)
    parse.(Int, collect(l))
end

const data = reduce(hcat, parse_line.(readlines()))

neighbours(pos) = [
    CartesianIndex(Tuple(pos) .- n)
    for n ∈ ((-1, 0), (1, 0), (0, -1), (0, 1))
]

# Part One
sum(CartesianIndices(data)) do i
    neib = [get(data, neib, Inf) for neib ∈ neighbours(i)]
    all(
        data[i] < n for n in neib
    ) ? data[i] + 1 : 0
end |> println

function dfs(pos, seen::Set)
    curr_n = 1
    push!(seen, pos)
    for neib in neighbours(pos)
        neib ∉ CartesianIndices(data) && continue
        if neib ∉ seen
            curr_n += dfs(neib, seen)
        end
    end
    curr_n
end

function solve_2()
    seen_pos = Set(findall(==(9), data))
    remain_pos = setdiff(CartesianIndices(data), seen_pos)
    basins = []
    while !isempty(remain_pos)
        next_pos = pop!(remain_pos)
        push!(basins, dfs(next_pos, seen_pos))
        setdiff!(remain_pos, seen_pos)
    end
    partialsort!(basins, 3; rev=true)
    prod(basins[1:3])
end

# Part Two
solve_2() |> println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
"""
**SeisRead**
*Read seismic data in seis format. The format is comprised of three elements:*
* a text file with information about data extent, data and header file names
* a binary file containing data
* a binary file containing headers
**
**IN**
* filename
* group="all" ("some" or "gather")
* key=["imx","imy"]
* itrace=1
* ntrace=10000
**OUT**
* d: data as 2d array
* h: headers as 1d array
* extent: extent of the data (try _fieldnames(Extent)_ to see the information this contains)
*Credits: AS, 2015*
"""
function SeisRead(filename;group="all",key=["imx","imy"],itrace=1,ntrace=10000)

    filename_data = ParseDataName(filename)
    filename_headers = ParseHeaderName(filename)
    extent = ReadTextHeader(filename)
    #println(filename_data)
    stream_d = open(filename_data)
    dtype = ParseDataFormat(filename)
    dtype = dtype == "native_float" ? Float32 : Complex{Float32}
    esize = ParseDataESize(filename)
    total = convert(Int,filesize(stream_d)/esize)
    close(stream_d)
    @compat nhead = length(fieldnames(Header))
    curr = zeros(length(key),1)
    prev = 1*curr
    nx = extent.n2*extent.n3*extent.n4*extent.n5
    if filename_headers != "NULL"
	stream_h = open(filename_headers)
	if (group == "all")
	    ntrace = nx
	elseif (group == "gather")
	    j = 1
	    for j=itrace:nx
		h1 = GrabHeader(stream_h,j)
		for ikey=1:length(key)
		    curr[ikey] = getfield(h1,Symbol(key[ikey]))
				end
		if curr != prev && j > itrace
		    ntrace = j - itrace
		    break
		end
		prev = 1*curr
	    end
	    ntrace = j < nx ? j - itrace : j - itrace + 1
	else
	    ntrace = nx - itrace + 1 > ntrace ? ntrace : nx - itrace + 1
	end
	position_h = 4*nhead*(itrace-1)
	seek(stream_h,position_h)
	h1 = read(stream_h,Header32Bits,nhead*ntrace)
	h1 = reshape(h1,nhead,convert(Int64,ntrace))
	h = Header[]
	for j = 1 : ntrace
	    h = push!(h,BitsToHeader(h1[:,j]))
	end
	close(stream_h)
    else
	ntrace = nx > ntrace ? ntrace : nx
    end
	stream_d = open(filename_data)
    position_d = 4*extent.n1*(itrace-1)
    seek(stream_d,position_d)
    d = read(stream_d,dtype,extent.n1*ntrace)
    if group=="all"
	if extent.n5 == 1 && extent.n4 == 1 && extent.n3 == 1 && extent.n2 == 1
	    d = reshape(d,Int(extent.n1))
	elseif extent.n5 == 1 && extent.n4 == 1 && extent.n3 == 1
	    d = reshape(d,convert(Int,extent.n1),convert(Int,extent.n2))
	elseif extent.n5 == 1 && extent.n4 == 1
	    d = reshape(d,convert(Int,extent.n1),convert(Int,extent.n2),convert(Int,extent.n3))
	elseif extent.n5 == 1
	    d = reshape(d,convert(Int,extent.n1),convert(Int,extent.n2),convert(Int,extent.n3),convert(Int,extent.n4))
	else
	    d = reshape(d,convert(Int,extent.n1),convert(Int,extent.n2),convert(Int,extent.n3),convert(Int,extent.n4),convert(Int,extent.n5))
	end
    else
	d = reshape(d,convert(Int64,extent.n1),convert(Int64,ntrace))
	end
    close(stream_d)
    if filename_headers != "NULL"
	return d,h,extent
    else
	return d,extent
    end
end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
using Logging

metafmt(level, _module, group, id, file, line) =
    Logging.default_metafmt(level, nothing, group, id, nothing, nothing)

submodules = [
    "players",
    "transactions",
    "call_raise_validation",
    "human_player_option",
    "table",
    "game",
    "play",
    "fuzz_play",
]

local_run = isempty(get(ENV, "CI", ""))
n_fuzz = local_run ? 10 : 100
n_fuzz_10_players = local_run ? 1 : 10
n_fuzz_3_players = local_run ? 10 : 20

tests_to_debug = ["play", "fuzz_play"]
# tests_to_debug = submodules

for submodule in submodules
    println("Starting tests for $submodule")
    t = 0
    local logger
    if any(submodule in tests_to_debug)
        logger = ConsoleLogger(stderr,Logging.Info)
        # logger = ConsoleLogger(stderr,Logging.Debug;meta_formatter=metafmt)
    else
        logger = NullLogger()
    end
    with_logger(logger) do
        t = @elapsed include(submodule*".jl")
    end
    println("Completed tests for $submodule, $(round(Int, t)) seconds elapsed")
end


text:
using XDiff
using ReverseDiff: GradientTape, GradientConfig, gradient, gradient!, compile
using Base.Test
# using BenchmarkTools

function test_compare(f; inputs...)
    vals = ([val for (name, val) in inputs]...)
    df = xdiff(f; inputs...)
    # compare_test runs in an older world age than the generated function
    # need to invoke the latest one in the test
    dvals = Base.invokelatest(df, vals...)
    dvals_a = [dvals...]

    f_tape = GradientTape(f, vals)
    compiled_f_tape = compile(f_tape)
    cfg = GradientConfig(vals)
    results = map(similar, vals)
    gradient!(results, compiled_f_tape, vals)
    results_a = [results...]
    @test isapprox(results_a, dvals_a[2:end]; atol=0.1)
end

logistic(x) = 1 ./ (1 + exp.(-x))
@scalardiff logistic(x::Number) 1 (logistic(x) .* (1 .- logistic(x)))

include("vae.jl")
include("aggr.jl")
include("linreg.jl")
include("ann.jl")
include("autoencoder.jl")
include("others.jl")
include("destruct.jl")
include("special.jl")


text:
## Load data for testing
using Plots
using SparseArrays
using LinearAlgebra
using TextParse
using DelimitedFiles
using ColorSchemes

function read_graph(filename)
    src, dst, val = csvread(filename;
        spacedelim=true, colparsers=[Int,Int,Float64], header_exists=false)[1]
    src .+= 1
    dst .+= 1
    n = maximum(src)
    n = max(maximum(dst),n)
    A = sparse(src, dst, val, n, n)
    A = max.(A,A')
    return findnz(triu(A,1))[1:2]
end

src, dst = read_graph("../../dataset/lawlor-spectra-k32.edgelist.gz")
# x, y = csvread("$../../dataset/lawlor-spectra-k32.edgelist.xy";
#         spacedelim=true, colparsers=[Float64,Float64],
#         header_exists=false)[1]
V = readdlm("../../dataset/lawlor-spectra-k32.coords")
V[:,1] *= -10
V[:,2] *= -4

function plotperm(v)
  return invperm(sortperm(v))
end

x = Array{Float64}(plotperm(V[:,1]))
y = Array{Float64}(plotperm(V[:,2]))

## Form a random subset of edges
using Random
Random.seed!(0)
p = randperm(length(src))[1:10^5] # random set of 1M edges
srcsub = Int.(src[p])
dstsub = Int.(dst[p])

## Hexbin graph plots
include("hexbinplots_orig.jl")
function shapecoords(hx::HexBinPlots.HexHistogram)
  h,vh = HexBinPlots.make_shapes(hx)
  vmax = maximum(vh)
  xsh = Vector{Float64}()
  ysh = Vector{Float64}()
  for k in eachindex(h)
      append!(xsh,h[k].x)
      push!(xsh,h[k].x[1])
      push!(xsh,NaN)
      append!(ysh,h[k].y)
      push!(ysh,h[k].y[1])
      push!(ysh,NaN)
  end
  return (xsh, ysh, vh)
end
using Interpolations
function control_point(xi, xj, yi, yj, dist_from_mid)
    xmid = 0.5 * (xi+xj)
    ymid = 0.5 * (yi+yj)
    # get the angle of y relative to x
    theta = atan((yj-yi) / (xj-xi)) + 0.5pi
    # dist = sqrt((xj-xi)^2 + (yj-yi)^2)
    # dist_from_mid = curvature_scalar * 0.5dist
    # now we have polar coords, we can compute the position, adding to the midpoint
    (xmid + dist_from_mid * cos(theta),
     ymid + dist_from_mid * sin(theta))
end
function _add_edge_points!(xlines,ylines,xsi,ysi,xdi,ydi,curve_amount,interplen)
  xpt, ypt = control_point(xsi,xdi,ysi,ydi,curve_amount)
  xpts = [xsi, xpt, xdi]
  ypts = [ysi, ypt, ydi]
  t = range(0, stop=1, length=3)
  A = hcat(xpts, ypts)
  itp = scale(interpolate(A, BSpline(Cubic(Natural(OnGrid())))), t, 1:2)
  tfine = range(0, 1, length=interplen)
  for t in tfine
    push!(xlines, itp(t, 1))
  end
  for t in tfine
    push!(ylines, itp(t, 2))
  end
end

# Mutating version from HexBinPlots
function _xy2counts!(counts::Dict{(Tuple{Int, Int}), Int},
                    x::AbstractArray,y::AbstractArray,
                    xsize::Real,ysize::Real,x0,y0)
  @inbounds for i in eachindex(x)
      h = convert(HexagonOffsetOddR,
                  cube_round(x[i] - x0,y[i] - y0,xsize, ysize))
      idx = (h.q, h.r)
      counts[idx] = 1 + get(counts,idx,0)
  end
  counts
end
using Statistics
using Hexagons
# Generate counts for a graph edge plot histogram
# hhparam = (xsize,ysize,x0,y0) are the hex-histogrma params
function graphedgehist(src,dst,x,y,ptscale,hhparams,curvature = 0.05)
  xlines = similar(x,0)
  ylines = similar(y,0)
  sizehint!(xlines,max(length(src)+1000,10000000)) # we'll process them in batches this big
  sizehint!(ylines,max(length(src)+1000,10000000))

  # this is the counts of the edge vectors...
  edgecounts = Dict{(Tuple{Int, Int}), Int}()

  # now, run through each edges...
  for ei in eachindex(src)
    si = src[ei]
    di = dst[ei]
    xsi = x[si]
    ysi = y[si]
    xdi = x[di]
    ydi = y[di]
    curve_amount = curvature*sign((di-si))

    edgelen = sqrt((xdi-xsi)^2 + (ydi-ysi)^2)
    interplen = max(ceil(Int,2*edgelen/ptscale),2)

    _add_edge_points!(xlines, ylines, xsi, ysi, xdi, ydi, curve_amount, interplen)

    if length(xlines) >= length(src)
      # compress all to edgecounts
      println("Compressing counts... $(ei) / $(length(src))")
      _xy2counts!(edgecounts, xlines, ylines, hhparams...)
      resize!(xlines,0)
      resize!(ylines,0)
    end
  end
  # handle final update
  return _xy2counts!(edgecounts, xlines, ylines, hhparams...)
end

function hexbingraphplot(src,dst,x,y;nbins=500)
  xbins = nbins
  ybins = nbins

  # compute params for hex histogram bins
  xmin,xmax = extrema(x)
  ymin,ymax = extrema(y)
  xspan, yspan = xmax - xmin, ymax - ymin
  xsize, ysize = xspan / xbins, yspan / ybins
  x0, y0 = xmin - xspan / 2,ymin - yspan / 2
  hhparams = (xsize, ysize, x0, y0) # save them as a tuple

  edgescale = 6
  hsize = min((xmax-xmin)/(xbins*edgescale),(ymax-ymin)/(ybins*edgescale))

  # need to compute hhparams

  # get coords for edges
  counts = graphedgehist(src,dst,x,y,hsize,hhparams)
  # do this for the nodes...
  #counts = HexBinPlots.xy2counts_(x,y,hhparams...)

  xh,yh,vh = HexBinPlots.counts2xy_(counts,hhparams...)
  hedges = HexBinPlots.HexHistogram(xh,yh,vh,xsize,ysize,false)

  #hx = fit(HexBinPlots.HexHistogram,x,y,nbins)
  xsh,ysh,vh = shapecoords(hedges)

  #plot(xsh, ysh, seriestype=:shape, alpha=log10.(vh)/maximum(log10.(vh)), fill_z=vh, legend=false,linealpha=0)
  return hedges
end
@time hedges = hexbingraphplot(src,dst,x,y,nbins=1000)
##
xsh,ysh,vh = shapecoords(hedges)

val_range = (minimum(V[:,1]),0.008)
colors = [ColorSchemes.get(ColorSchemes.CMRmap,i,val_range) for i in V[:,1]]

axform = (log10.(vh).^3)/maximum(log10.(vh).^3)
avals = axform/maximum(axform)
plot(xsh, ysh, seriestype=:shape, alpha=avals, fillcolor=:darkblue, legend=false,linealpha=0, framestyle=:none,background=nothing)
scatter!(x,y,markercolor=colors,alpha=0.1,markersize=1,legend=false,colorbar=false,markerstrokewidth=0,background=nothing)
##
plot!(dpi=400,size=(1200,1200)) # size is inches times 100, so this should be about 3600x3600
png("hexbin_plot_global_high_res.png")

axform = (log10.(vh).^3)/maximum(log10.(vh).^3)
avals = axform/maximum(axform)
plot(xsh, ysh, seriestype=:shape, alpha=avals, fillcolor=:darkblue, legend=false,linealpha=0, framestyle=:none,background=nothing)
scatter!(x,y,markercolor=colors,alpha=0.1,markersize=1,legend=false,colorbar=false,markerstrokewidth=0,background=nothing)
plot!(dpi=150,size=(1200,1200)) # size is inches times 100, so this should be about 3600x3600
png("hexbin_plot_global_low_res.png")

max_stars_repo_name: dgleich/flowpaper-code

max_issues_repo_name: dgleich/flowpaper-code

max_forks_repo_name: dgleich/flowpaper-code


text:
module TestQuadsum

using QuantEcon
using Base.Test
using FactCheck

rough_kwargs = Dict(:atol => 1e-10, :rtol => 1e-10)

facts("Testing quadsums.jl") do
    context("test simple var sum") do
        beta = .95
        A = 1.
        C = 0.
        H = 1.
        x0 = 1.

        val = var_quadratic_sum(A, C, H, beta, x0)

        @fact val --> roughly(20.0; rough_kwargs...)
    end

    context("test identity var sum") do
        beta = .95
        A = eye(3)
        C = zeros(3, 3)
        H = eye(3)
        x0 = ones(3)

        val = var_quadratic_sum(A, C, H, beta, x0)

        @fact val --> roughly(60.0; rough_kwargs...)
    end

end  # facts
end  # module