text:
struct ServerAction
    id::String
    desc::String
    kind::Union{CodeActionKind,Missing}
    preferred::Union{Bool,Missing}
    when::Function
    handler::Function
end

function client_support_action_kind(s::LanguageServerInstance, _::CodeActionKind=CodeActionKinds.Empty)
    if s.clientCapabilities !== missing &&
       s.clientCapabilities.textDocument !== missing &&
       s.clientCapabilities.textDocument.codeAction !== missing &&
       s.clientCapabilities.textDocument.codeAction.codeActionLiteralSupport !== missing
       s.clientCapabilities.textDocument.codeAction.codeActionLiteralSupport.codeActionKind !== missing
        # From the spec of CodeActionKind (https://microsoft.github.io/language-server-protocol/specifications/lsp/3.17/specification/#codeActionClientCapabilities):
        #
        #     The code action kind values the client supports. When this
        #     property exists the client also guarantees that it will
        #     handle values outside its set gracefully and falls back
        #     to a default value when unknown.
        #
        # so we can always return true here(?).
        return true
    else
        return false
    end
end

function client_preferred_support(s::LanguageServerInstance)::Bool
    if s.clientCapabilities !== missing &&
       s.clientCapabilities.textDocument !== missing &&
       s.clientCapabilities.textDocument.codeAction !== missing &&
       s.clientCapabilities.textDocument.codeAction.isPreferredSupport !== missing
       return s.clientCapabilities.textDocument.codeAction.isPreferredSupport
   else
       return false
    end
end

# TODO: All currently supported CodeActions in LS.jl can be converted "losslessly" to
#       Commands but this might not be true in the future so unless the client support
#       literal code actions those need to be filtered out.
function convert_to_command(ca::CodeAction)
    return ca.command
end

function textDocument_codeAction_request(params::CodeActionParams, server::LanguageServerInstance, conn)
    actions = CodeAction[]
    doc = getdocument(server, params.textDocument.uri)
    offset = index_at(doc, params.range.start)
    x = get_expr(getcst(doc), offset)
    arguments = Any[params.textDocument.uri, offset] # use the same arguments for all commands
    if x isa EXPR
        for (_, sa) in LSActions
            if sa.when(x, params)
                kind = sa.kind
                if sa.kind !== missing && sa.kind == CodeActionKinds.SourceOrganizeImports &&
                    server.clientInfo !== missing && occursin("code", lowercase(server.clientInfo.name))
                    # SourceOrganizeImports doesn't show up in the VS Code UI, so make this a
                    # RefactorRewrite instead
                    kind = CodeActionKinds.RefactorRewrite
                end
                action = CodeAction(
                    sa.desc, # title
                    kind,    # kind
                    missing, # diagnostics
                    client_preferred_support(server) ? sa.preferred : missing, # isPreferred
                    missing, # edit
                    Command(sa.desc, sa.id, arguments), # command
                )
                push!(actions, action)
            end
        end
    end
    if client_support_action_kind(server)
        return actions
    else
        # TODO: Future CodeActions might have to be filtered here.
        return convert_to_command.(actions)
    end
end

function workspace_executeCommand_request(params::ExecuteCommandParams, server::LanguageServerInstance, conn)
    if haskey(LSActions, params.command)
        uri = URI(params.arguments[1])
        offset = params.arguments[2]
        doc = getdocument(server, uri)
        x = get_expr(getcst(doc), offset)
        LSActions[params.command].handler(x, server, conn)
    end
end

function find_using_statement(x::EXPR)
    for ref in refof(x).refs
        if StaticLint.is_in_fexpr(ref, x -> headof(x) === :using || headof(x) === :import)
            return parentof(ref)
        end
    end
    return nothing
end

function explicitly_import_used_variables(x::EXPR, server, conn)
    !(refof(x) isa StaticLint.Binding && refof(x).val isa SymbolServer.ModuleStore) && return
    using_stmt = find_using_statement(x)
    using_stmt isa Nothing && return

    tdes = Dict{URI,TextDocumentEdit}()
    vars = Set{String}() # names that need to be imported
    # Find uses of `x` and mark edits
    for ref in refof(x).refs
        if parentof(ref) isa EXPR && CSTParser.is_getfield_w_quotenode(parentof(ref)) && parentof(ref).args[1] == ref
            childname = parentof(ref).args[2].args[1]
            StaticLint.hasref(childname) && refof(childname) isa StaticLint.Binding && continue # check this isn't the name of something being explictly overwritten
            !haskey(refof(x).val.vals, Symbol(valof(childname))) && continue # skip, perhaps mark as missing ref ?

            file, offset = get_file_loc(ref)
            if !haskey(tdes, get_uri(file))
                tdes[get_uri(file)] = TextDocumentEdit(VersionedTextDocumentIdentifier(get_uri(file), get_version(file)), TextEdit[])
            end
            push!(tdes[get_uri(file)].edits, TextEdit(Range(file, offset .+ (0:parentof(ref).span)), valof(childname)))
            push!(vars, valof(childname))
        end
    end
    isempty(tdes) && return

    # Add `using x: vars...` statement
    if parentof(using_stmt) isa EXPR && (headof(parentof(using_stmt)) === :block || headof(parentof(using_stmt)) === :file)
        # this should cover all cases
        i1 = 0
        for i = 1:length(parentof(using_stmt).args)
            if using_stmt === parentof(using_stmt).args[i]
                i1 = i
                break
            end
        end
        i1 == 0 && return WorkspaceEdit(missing, missing)

        file, offset = get_file_loc(using_stmt)
        insertpos = get_next_line_offset(using_stmt)
        insertpos == -1 && return

        if !haskey(tdes, get_uri(file))
            tdes[get_uri(file)] = TextDocumentEdit(VersionedTextDocumentIdentifier(get_uri(file), get_version(file)), TextEdit[])
        end
        push!(tdes[get_uri(file)].edits, TextEdit(Range(file, insertpos .+ (0:0)), string("using ", valof(x), ": ", join(vars, ", "), "\n")))
    else
        return
    end

    JSONRPC.send(conn, workspace_applyEdit_request_type, ApplyWorkspaceEditParams(missing, WorkspaceEdit(missing, collect(values(tdes)))))
end

is_single_line_func(x) = CSTParser.defines_function(x) && headof(x) !== :function

function expand_inline_func(x, server, conn)
    func = _get_parent_fexpr(x, is_single_line_func)
    length(func) < 3 && return
    sig = func.args[1]
    op = func.head
    body = func.args[2]
    if headof(body) == :block && length(body) == 1
        file, offset = get_file_loc(func)
        tde = TextDocumentEdit(VersionedTextDocumentIdentifier(get_uri(file), get_version(file)), TextEdit[
            TextEdit(Range(file, offset .+ (0:func.fullspan)), string("function ", get_text(file)[offset .+ (1:sig.span)], "\n    ", get_text(file)[offset + sig.fullspan + op.fullspan .+ (1:body.span)], "\nend\n"))
        ])
        JSONRPC.send(conn, workspace_applyEdit_request_type, ApplyWorkspaceEditParams(missing, WorkspaceEdit(missing, TextDocumentEdit[tde])))
    elseif (headof(body) === :begin || CSTParser.isbracketed(body)) &&
        headof(body.args[1]) === :block && length(body.args[1]) > 0
        file, offset = get_file_loc(func)
        newtext = string("function ", get_text(file)[offset .+ (1:sig.span)])
        blockoffset = offset + sig.fullspan + op.fullspan + body.trivia[1].fullspan
        for i = 1:length(body.args[1].args)
            newtext = string(newtext, "\n    ", get_text(file)[blockoffset .+ (1:body.args[1].args[i].span)])
            blockoffset += body.args[1].args[i].fullspan
        end
        newtext = string(newtext, "\nend\n")
        tde = TextDocumentEdit(VersionedTextDocumentIdentifier(get_uri(file), get_version(file)), TextEdit[TextEdit(Range(file, offset .+ (0:func.fullspan)), newtext)])
        JSONRPC.send(conn, workspace_applyEdit_request_type, ApplyWorkspaceEditParams(missing, WorkspaceEdit(missing, TextDocumentEdit[tde])))
    end
end

function is_in_fexpr(x::EXPR, f)
    if f(x)
        return true
    elseif parentof(x) isa EXPR
        return is_in_fexpr(parentof(x), f)
    else
        return false
    end
end

function _get_parent_fexpr(x::EXPR, f)
    if f(x)
        return x
    elseif parentof(x) isa EXPR
        return _get_parent_fexpr(parentof(x), f)
    end
end
function get_next_line_offset(x)
    file, offset = get_file_loc(x)
    # get next line after using_stmt
    insertpos = -1
    line_offsets = get_line_offsets(get_text_document(file))
    for i = 1:length(line_offsets) - 1
        if line_offsets[i] < offset + x.span <= line_offsets[i + 1]
            insertpos = line_offsets[i + 1]
        end
    end
    return insertpos
end

function reexport_package(x::EXPR, server, conn)
    (refof(x) isa SymbolServer.ModuleStore || refof(x).type === StaticLint.CoreTypes.Module || (refof(x).val isa StaticLint.Binding && refof(x).val.type === StaticLint.CoreTypes.Module)) || (refof(x).val isa SymbolServer.ModuleStore) || return
    mod = if refof(x) isa SymbolServer.ModuleStore
        refof(x)
    elseif refof(x).val isa SymbolServer.ModuleStore
        refof(x).val
    else
        return
    end
    using_stmt = parentof(x)
    file, _ = get_file_loc(x)
    insertpos = get_next_line_offset(using_stmt)
    insertpos == -1 && return

    text_document = get_text_document(file)

    tde = TextDocumentEdit(VersionedTextDocumentIdentifier(get_uri(text_document), get_version(text_document)), TextEdit[
        TextEdit(Range(file, insertpos .+ (0:0)), string("export ", join(sort([string(n) for (n, v) in mod.vals if StaticLint.isexportedby(n, mod)]), ", "), "\n"))
    ])

    JSONRPC.send(conn, workspace_applyEdit_request_type, ApplyWorkspaceEditParams(missing, WorkspaceEdit(missing, TextDocumentEdit[tde])))
end

# TODO move to StaticL  int
# to be called where typof(x) === CSTParser.ModuleH/BareModule
function find_exported_names(x::EXPR)
    exported_vars = EXPR[]
    for i in 1:length(x.args[3].args)
        expr = x.args[3].args[i]
        if headof(expr) === :export
            for j = 2:length(expr)
                if CSTParser.isidentifier(expr.args[j]) && StaticLint.hasref(expr.args[j])
                    push!(exported_vars, expr.args[j])
                end
            end
        end
    end
    return exported_vars
end

function reexport_module(x::EXPR, server, conn)
    using_stmt = parentof(x)
    mod_expr = refof(x).val isa StaticLint.Binding ? refof(x).val.val : refof(x).val
    (mod_expr.args isa Nothing || length(mod_expr.args) < 3 || headof(mod_expr.args[3]) !== :block || mod_expr.args[3].args isa Nothing) && return # module expr without block
    # find export EXPR
    exported_names = find_exported_names(mod_expr)

    isempty(exported_names) && return
    file, _ = get_file_loc(x)
    insertpos = get_next_line_offset(using_stmt)
    insertpos == -1 && return
    names = filter!(s -> !isempty(s), collect(CSTParser.str_value.(exported_names)))
    tde = TextDocumentEdit(VersionedTextDocumentIdentifier(get_uri(file), get_version(file)), TextEdit[
        TextEdit(Range(file, insertpos .+ (0:0)), string("export ", join(sort(names), ", "), "\n"))
    ])

    JSONRPC.send(conn, workspace_applyEdit_request_type, ApplyWorkspaceEditParams(missing, WorkspaceEdit(missing, TextDocumentEdit[tde])))
end

function wrap_block(x, server, type, conn) end
function wrap_block(x::EXPR, server, type, conn)
    file, offset = get_file_loc(x) # rese
    if type == :if
        tde = TextDocumentEdit(VersionedTextDocumentIdentifier(get_uri(file), get_version(file)), TextEdit[
            TextEdit(Range(file, offset .+ (0:0)), "if CONDITION\n"),
            TextEdit(Range(file, offset + x.span .+ (0:0)), "\nend")
        ])
    end

    JSONRPC.send(conn, workspace_applyEdit_request_type, ApplyWorkspaceEditParams(missing, WorkspaceEdit(missing, TextDocumentEdit[tde])))
end

function is_fixable_missing_ref(x::EXPR, cac::CodeActionContext)
    if !isempty(cac.diagnostics) && any(startswith(d.message, "Missing reference") for d::Diagnostic in cac.diagnostics) && CSTParser.isidentifier(x)
        xname = StaticLint.valofid(x)
        tls = StaticLint.retrieve_toplevel_scope(x)
        if tls isa StaticLint.Scope && tls.modules !== nothing
            for m in values(tls.modules)
                if (m isa SymbolServer.ModuleStore && haskey(m, Symbol(xname))) || (m isa StaticLint.Scope && StaticLint.scopehasbinding(m, xname))
                    return true
                end
            end
        end
    end
    return false
end

function applymissingreffix(x, server, conn)
    xname = StaticLint.valofid(x)
    file, offset = get_file_loc(x)
    tls = StaticLint.retrieve_toplevel_scope(x)
    if tls.modules !== nothing
        for (n, m) in tls.modules
            if (m isa SymbolServer.ModuleStore && haskey(m, Symbol(xname))) || (m isa StaticLint.Scope && StaticLint.scopehasbinding(m, xname))
                tde = TextDocumentEdit(VersionedTextDocumentIdentifier(get_uri(file), get_version(file)), TextEdit[
                    TextEdit(Range(file, offset .+ (0:0)), string(n, "."))
                ])
                JSONRPC.send(conn, workspace_applyEdit_request_type, ApplyWorkspaceEditParams(missing, WorkspaceEdit(missing, TextDocumentEdit[tde])))
            end
        end
    end
end

function remove_farg_name(x, server, conn)
    x1 = StaticLint.get_parent_fexpr(x, x -> StaticLint.haserror(x) && StaticLint.errorof(x) == StaticLint.UnusedFunctionArgument)
    file, offset = get_file_loc(x1)
    if CSTParser.isdeclaration(x1)
        tde = TextDocumentEdit(VersionedTextDocumentIdentifier(get_uri(file), get_version(file)), TextEdit[
                        TextEdit(Range(file, offset .+ (0:x1.args[1].fullspan)), "")
                    ])
    else
        tde = TextDocumentEdit(VersionedTextDocumentIdentifier(get_uri(file), get_version(file)), TextEdit[
                        TextEdit(Range(file, offset .+ (0:x1.fullspan)), "_")
                    ])
    end
    JSONRPC.send(conn, workspace_applyEdit_request_type, ApplyWorkspaceEditParams(missing, WorkspaceEdit(missing, TextDocumentEdit[tde])))
end

function remove_unused_assignment_name(x, _, conn)
    x1 = StaticLint.get_parent_fexpr(x, x -> StaticLint.haserror(x) && StaticLint.errorof(x) == StaticLint.UnusedBinding && x isa EXPR && x.head === :IDENTIFIER)
    file, offset = get_file_loc(x1)
    tde = TextDocumentEdit(VersionedTextDocumentIdentifier(get_uri(file), get_version(file)), TextEdit[
                    TextEdit(Range(file, offset .+ (0:x1.span)), "_")
                ])
    JSONRPC.send(conn, workspace_applyEdit_request_type, ApplyWorkspaceEditParams(missing, WorkspaceEdit(missing, TextDocumentEdit[tde])))
end

function double_to_triple_equal(x, _, conn)
    x1 = StaticLint.get_parent_fexpr(x, y -> StaticLint.haserror(y) && StaticLint.errorof(y) in (StaticLint.NothingEquality, StaticLint.NothingNotEq))
    file, offset = get_file_loc(x1)
    tde = TextDocumentEdit(VersionedTextDocumentIdentifier(get_uri(file), get_version(file)), TextEdit[
        TextEdit(Range(file, offset .+ (0:x1.span)), StaticLint.errorof(x1) == StaticLint.NothingEquality ? "===" : "!==")
    ])
    JSONRPC.send(conn, workspace_applyEdit_request_type, ApplyWorkspaceEditParams(missing, WorkspaceEdit(missing, TextDocumentEdit[tde])))
end

function get_spdx_header(doc::Document)
    m = match(r"(*ANYCRLF)^# SPDX-License-Identifier:\s+((?:[\w\.-]+)(?:\s+[\w\.-]+)*)\s*$"m, get_text(doc))
    return m === nothing ? m : String(m[1])
end

function in_same_workspace_folder(server::LanguageServerInstance, file1::URI, file2::URI)
    for ws in server.workspaceFolders
        if startswith(uri2filepath(file1), ws) &&
           startswith(uri2filepath(file2), ws)
           return true
       end
    end
    return false
end

function identify_short_identifier(server::LanguageServerInstance, file::Document)
    # First look in tracked files (in the same workspace folder) for existing headers
    candidate_identifiers = Set{String}()
    for doc in getdocuments_value(server)
        in_same_workspace_folder(server, get_uri(file), get_uri(doc)) || continue
        id = get_spdx_header(doc)
        id === nothing || push!(candidate_identifiers, id)
    end
    if length(candidate_identifiers) == 1
        return first(candidate_identifiers)
    end
    # Fallback to looking for a license file in the same workspace folder
    candidate_files = String[]
    for dir in server.workspaceFolders
        for f in joinpath.(dir, ["LICENSE", "LICENSE.md"])
            if in_same_workspace_folder(server, get_uri(file), filepath2uri(f)) && safe_isfile(f)
                push!(candidate_files, f)
            end
        end
    end
    length(candidate_files) == 1 || return nothing
    license = read(first(candidate_files), String)

    # This is just a heuristic, but should be OK since this is not something automated, and
    # the programmer will see directly if the wrong license is added.
    # TODO: Add more licenses...
    if contains(license, r"MIT\s+(\"?Expat\"?\s+)?License")
        return "MIT"
    end
    return nothing
end

function add_license_header(x, server::LanguageServerInstance, conn)
    file, _ = get_file_loc(x)
    get_spdx_header(file) === nothing || return # TODO: Would be nice to check this already before offering the action
    short_identifier = identify_short_identifier(server, file)
    short_identifier === nothing && return
    tde = TextDocumentEdit(VersionedTextDocumentIdentifier(get_uri(file), get_version(file)), TextEdit[
        TextEdit(Range(file, 0:0), "# SPDX-License-Identifier: $(short_identifier)\n\n")
    ])
    JSONRPC.send(conn, workspace_applyEdit_request_type, ApplyWorkspaceEditParams(missing, WorkspaceEdit(missing, TextDocumentEdit[tde])))
end

function organize_import_block(x, _, conn)
    if !StaticLint.is_in_fexpr(x, x -> headof(x) === :using || headof(x) === :import)
        return
    end

    # Collect all sibling blocks that are also using/import expressions
    siblings = EXPR[]
    using_stmt = StaticLint.get_parent_fexpr(x, y -> headof(y) === :using || headof(y) === :import)
    push!(siblings, using_stmt)
    block = using_stmt.parent
    if block !== nothing
        myidx = findfirst(x -> x === using_stmt, block.args)
        # Find older and younger siblings
        for direction in (-1, 1)
            i = direction
            while true
                s = get(block.args, myidx + i, nothing)
                if s isa EXPR && (s.head === :using || s.head === :import)
                    (direction == 1 ? push! : pushfirst!)(siblings, s)
                    i += direction
                else
                    break
                end
            end
        end
    end

    # Collect all modules and symbols
    using_mods = Set{String}()
    using_syms = Dict{String,Set{String}}()
    import_mods = Set{String}()
    import_syms = Dict{String,Set{String}}()

    for s in siblings
        isusing = s.head === :using
        for a in s.args
            if CSTParser.is_colon(a.head)
                mod = a.args[1].args[1].val
                set = get!(Set, isusing ? using_syms : import_syms, mod)
                for i in 2:length(a.args)
                    push!(set, join(y.val for y in a.args[i]))
                end
            elseif CSTParser.is_dot(a.head)
                push!(isusing ? using_mods : import_mods, join((y.val for y in a.args), "."))
            else
                error("unreachable?")
            end
        end
    end

    # Rejoin and sort
    # TODO: Currently regular string sorting is used, which roughly will correspond to
    #       BlueStyle (modules, types, ..., functions) since usually CamelCase is used for
    #       modules, types, etc, but possibly this can be improved by using information
    #       available from SymbolServer
    import_lines = String[]
    for m in import_mods
        push!(import_lines, "import " * m)
    end
    for (m, s) in import_syms
        push!(import_lines, "import " * m * ": " * join(sort!(collect(s)), ", "))
    end
    using_lines = String[]
    for m in using_mods
        push!(using_lines, "using " * m)
    end
    for (m, s) in using_syms
        push!(using_lines, "using " * m * ": " * join(sort!(collect(s)), ", "))
    end
    io = IOBuffer()
    join(io, sort!(import_lines), "\n")
    length(import_lines) > 0 && print(io, "\n\n")
    join(io, sort!(using_lines), "\n")
    str_to_fmt = String(take!(io))

    # Format the new string
    # TODO: Fetch user configuration?
    formatted = JuliaFormatter.format_text(str_to_fmt; join_lines_based_on_source=true)

    # Compute range of original blocks
    file, firstoffset = get_file_loc(first(siblings))
    _, lastoffset = get_file_loc(last(siblings))
    lastoffset += last(siblings).span

    tde = TextDocumentEdit(VersionedTextDocumentIdentifier(get_uri(file), get_version(file)), TextEdit[
        TextEdit(Range(file, firstoffset:lastoffset), formatted)
    ])
    JSONRPC.send(conn, workspace_applyEdit_request_type, ApplyWorkspaceEditParams(missing, WorkspaceEdit(missing, TextDocumentEdit[tde])))
end

# Adding a CodeAction requires defining:
# * a unique id
# * a description
# * an action kind (optionally)
# * a function (.when) called on the currently selected expression and parameters of the CodeAction call;
# * a function (.handler) called on three arguments (current expression, server and the jr connection) to implement the command.
const LSActions = Dict{String,ServerAction}()

LSActions["ExplicitPackageVarImport"] = ServerAction(
    "ExplicitPackageVarImport",
    "Explicitly import used package variables.",
    missing,
    missing,
    (x, params) -> refof(x) isa StaticLint.Binding && refof(x).val isa SymbolServer.ModuleStore,
    explicitly_import_used_variables
)

LSActions["ExpandFunction"] = ServerAction(
    "ExpandFunction",
    "Expand function definition.",
    CodeActionKinds.Refactor,
    missing,
    (x, params) -> is_in_fexpr(x, is_single_line_func),
    expand_inline_func,
)

LSActions["FixMissingRef"] = ServerAction(
    "FixMissingRef",
    "Fix missing reference",
    missing,
    missing,
    (x, params) -> is_fixable_missing_ref(x, params.context),
    applymissingreffix,
)

LSActions["ReexportModule"] = ServerAction(
    "ReexportModule",
    "Re-export package variables.",
    missing,
    missing,
    (x, params) -> StaticLint.is_in_fexpr(x, x -> headof(x) === :using || headof(x) === :import) && (refof(x) isa StaticLint.Binding && (refof(x).type === StaticLint.CoreTypes.Module || (refof(x).val isa StaticLint.Binding && refof(x).val.type === StaticLint.CoreTypes.Module) || refof(x).val isa SymbolServer.ModuleStore) || refof(x) isa SymbolServer.ModuleStore),
    reexport_package,
)

LSActions["DeleteUnusedFunctionArgumentName"] = ServerAction(
    "DeleteUnusedFunctionArgumentName",
    "Delete name of unused function argument.",
    CodeActionKinds.QuickFix,
    missing,
    (x, params) -> StaticLint.is_in_fexpr(x, x -> StaticLint.haserror(x) && StaticLint.errorof(x) == StaticLint.UnusedFunctionArgument),
    remove_farg_name,
)

LSActions["ReplaceUnusedAssignmentName"] = ServerAction(
    "ReplaceUnusedAssignmentName",
    "Replace unused assignment name with _.",
    CodeActionKinds.QuickFix,
    missing,
    (x, params) -> StaticLint.is_in_fexpr(x, x -> StaticLint.haserror(x) && StaticLint.errorof(x) == StaticLint.UnusedBinding && x isa EXPR && x.head === :IDENTIFIER),
    remove_unused_assignment_name,
)

LSActions["CompareNothingWithTripleEqual"] = ServerAction(
    "CompareNothingWithTripleEqual",
    "Change ==/!= to ===/!==.",
    CodeActionKinds.QuickFix,
    true,
    (x, _) -> StaticLint.is_in_fexpr(x, y -> StaticLint.haserror(y) && (StaticLint.errorof(y) in (StaticLint.NothingEquality, StaticLint.NothingNotEq))),
    double_to_triple_equal,
)

LSActions["AddLicenseIdentifier"] = ServerAction(
    "AddLicenseIdentifier",
    "Add SPDX license identifier.",
    missing,
    missing,
    (_, params) -> params.range.start.line == 0,
    add_license_header,
)

LSActions["OrganizeImports"] = ServerAction(
    "OrganizeImports",
    "Organize `using` and `import` statements.",
    CodeActionKinds.SourceOrganizeImports,
    true,
    (x, _) -> StaticLint.is_in_fexpr(x, x -> headof(x) === :using || headof(x) === :import),
    organize_import_block,
)


text:
fade(x::String) = HTML("$x")

icon(x::String) = HTML("")

function interleave(xs, j)
  ys = []
  for x in xs
    push!(ys, x, j)
  end
  isempty(xs) || pop!(ys)
  return ys
end

dims(xs...) = Row(interleave(xs, fade("×"))...)

const UNDEF = fade("#undef")

function undefs(xs)
  xs′ = similar(xs, Any)
  for i in eachindex(xs)
    xs′[i] = isassigned(xs, i) ? xs[i] : UNDEF
  end
  return xs′
end

errtrace(e, trace) = trace

errmsg(e) = sprint(io -> showerror(IOContext(io, :limit => true), e))

view(x) =
  Dict(:type    => :html,
       :content => stringmime(MIME"text/html"(), x))

ext: jl


text:
"""
    NonCommonVertices()

    Given the indices of vertices of a source simplex inside the circumsphere of a target simplex,
    find the indices of the vertices that are not shared between simplices, but still lie
    inside the circumference of the target simplex.

Input arguments
---------------
IndexIn::Array{Int, 1}              Row vector. Indices of the vertices of the source simplex inside the circumsphere of the target simplex.
InternalComIndex::Array{Int, 1}     Row vector. Indices of the vertices that both simplices share, but in the 'IndexIn' row vector of indices instead of the 'IndexComVert' vector.
Nin::Int                            Integer. Number of vertices inside the circumsphere of the target simplex. REDUNDANT? Essentially, length(IndexIn) - but remember to check for zeros
Ncomm::Int                          Integer. Number of common vertices between between source and target simplices.

Returns
-------
Res                                 Row vector. Indices of the vertices of the source simplex that are inside the circumsphere of the target simplex, but is not shared between the simplices.
"""

function NonCommonVertices(InternalComIndex,IndexIn,
                          n_vertices_in_circumsphere,
                          n_commonvertices)
  # IndexComVert: indices of the common vertices
  # IndexIn: indices of the vertices contained in the circumsphere
  # InternalComIndex: Indices of IndexComVert in IndexIn
  # Nin: number of vertices in the circumsphere

  # OUTCOME
  @show n_vertices_in_circumsphere
  @show n_commonvertices
  @show IndexIn
  @show InternalComIndex
  # Res: indices of the non common vertices still inside the circumsphere
  # if Res is empty, it returns 0
  if n_commonvertices == 0
    #println("n_commonvertices == 0")
    return IndexIn
  else
      if n_commonvertices == n_vertices_in_circumsphere
        #println("Number of vertices in circumsphere equals number of common vertices")
        return [0]
      else
        #println("Number of vertices in circumsphere does not equal number of common vertices")
        return IndexIn[complementary(InternalComIndex, n_vertices_in_circumsphere)]
     end
 end

end


text:
"""
Code for estimating the radial velocity based on projecting the flux onto the derivative of a template spectrum.

Author: Eric Ford
"""

""" Return mean flux (averaging over observations at different times, variance weighted) based on a common set of wavelengths.
   Inputs: flux & var (2d: pixel, time)
"""
function calc_mean_spectrum(flux::AbstractArray{T1,2}, var::AbstractArray{T2,2} ) where { T1<:Real, T2<:Real }
    flux_mean = vec(sum(flux./var,dims=2)./sum(1.0./var,dims=2))
end

""" Estimate numerical derivative of fluxes given wavelengths. """
function calc_dfluxdlnlambda(flux::AbstractArray{T1,1}, λ::AbstractArray{T2,1}) where { T1<:Real, T2<:Real }
    @assert size(flux) == size(λ)
    @assert length(flux) >= 3
    dfdlogλ = Array{T1,1}(undef,length(flux))
    calc_dfluxdlnlambda!(dfdlogλ,flux,λ)
    #dfdlogλ[1] = 0.5 * (flux[2]-flux[1])/(λ[2]-λ[1])*(λ[2]+λ[1])
    #dfdlogλ[2:end-1] .= 0.5 .* (flux[3:end].-flux[1:end-2])./(λ[3:end].-λ[1:end-2]).*(λ[3:end].+λ[1:end-2])
    #dfdlogλ[end] = 0.5 * (flux[end]-flux[end-1])/(λ[end]-λ[end-1])*(λ[end]+λ[end-1])
    return dfdlogλ
end

""" Estimate numerical derivative of fluxes given wavelengths. """
function calc_dfluxdlnlambda!(dfdlogλ::AbstractArray{T1,1}, flux::AbstractArray{T2,1}, λ::AbstractArray{T3,1}) where { T1<:Real, T2<:Real, T3<:Real }
    @assert size(flux) == size(λ)
    @assert size(dfdlogλ) == size(flux)
    @assert length(flux) >= 3
    #dfdlogλ = Array{T1,1}(undef,length(flux))
    dfdlogλ[1] = 0.5 * (flux[2]-flux[1])/(λ[2]-λ[1])*(λ[2]+λ[1])
    dfdlogλ[2:end-1] .= 0.5 .* (flux[3:end].-flux[1:end-2])./(λ[3:end].-λ[1:end-2]).*(λ[3:end].+λ[1:end-2])
    dfdlogλ[end] = 0.5 * (flux[end]-flux[end-1])/(λ[end]-λ[end-1])*(λ[end]+λ[end-1])
    return dfdlogλ
end

""" Estimate numerical second derivative of fluxes given wavelengths. """
function calc_d2fluxdlnlambda2(flux::AbstractArray{T1,1}, λ::AbstractArray{T2,1}) where { T1<:Real, T2<:Real }
    @assert size(flux) == size(λ)
    @assert length(flux) >= 3
    logλ = log.(λ)
    d2fdlogλ2 = Array{T1,1}(undef,length(flux))
    #d2fdlogλ2[2:end-1] .= 0.25*(flux[3:end].+flux[1:end-2].-2.0.*flux[2:end-1])./(λ[3:end].+λ[1:end-2].-2.0*λ[2:end-1]).*(λ[3:end].+λ[end-2]).^2
    #d2fdlogλ2[2:end-1] .= 0.5 * (flux[3:end].+flux[1:end-2].-2.0.*flux[2:end-1]).* ((λ[3:end].+λ[1:end-2])./(logλ[3:end].-logλ[1:end-2])).^2
    #d2fdlogλ2[end] = d2fdlogλ2[end-1]
    #d2fdlogλ2[1] = d2fdlogλ2[2]
    calc_d2fluxdlnlambda2!(d2fdlogλ2,flux,λ)
    return d2fdlogλ2
end

""" Estimate numerical second derivative of fluxes given wavelengths. """
function calc_d2fluxdlnlambda2!(d2fdlogλ2::AbstractArray{T1,1}, flux::AbstractArray{T1,1}, λ::AbstractArray{T2,1}) where { T1<:Real, T2<:Real, T3<:Real }
    @assert size(flux) == size(λ)
    @assert size(d2fdlogλ2) == size(flux)
    @assert length(flux) >= 3
    #logλ = log.(λ)
    #d2fdlogλ2 = Array{T1,1}(undef,length(flux))
    #d2fdlogλ2[2:end-1] .= 0.25*(flux[3:end].+flux[1:end-2].-2.0.*flux[2:end-1])./(λ[3:end].+λ[1:end-2].-2.0*λ[2:end-1]).*(λ[3:end].+λ[end-2]).^2
    d2fdlogλ2[2:end-1] .= 0.5 * (flux[3:end].+flux[1:end-2].-2.0.*flux[2:end-1]).* ((λ[3:end].+λ[1:end-2]).^2 ./ (2.0.*log.(λ[3:end]./λ[1:end-2])))
    d2fdlogλ2[end] = d2fdlogλ2[end-1]
    d2fdlogλ2[1] = d2fdlogλ2[2]
    return d2fdlogλ2
end

""" Return mean numerical derivative of fluxes based on a common set of wavelengths.
    Inputs: flux & var (2d) and λ (1d)
 """
function calc_mean_dfluxdlnlambda(flux::AbstractArray{T1,2}, var::AbstractArray{T1,2}, λ::AbstractArray{T3,1},
        chunk_map::AbstractArray{URT,1}) where
    { T1<:Real, T2<:Real, T3<:Real, URT<:AbstractUnitRange} #, RT<:AbstractRange }
    flux_mean = calc_mean_spectrum(flux,var)
    deriv = Array{T1,1}(undef,length(flux_mean))
    #map(c->deriv[c] .= calc_dfluxdlnlambda(flux_mean[c],λ[c]),chunk_map )
    for c in chunk_map
         calc_dfluxdlnlambda!(deriv[c],flux_mean[c],λ[c])
    end
    return deriv
end

function calc_mean_d2fluxdlnlambda2(flux::AbstractArray{T1,2}, var::AbstractArray{T1,2}, λ::AbstractArray{T3,1},
        chunk_map::AbstractArray{URT,1}) where
    { T1<:Real, T2<:Real, T3<:Real, URT<:AbstractUnitRange} #, RT<:AbstractRange }
    flux_mean = calc_mean_spectrum(flux,var)
    deriv2 = Array{T1,1}(undef,length(flux_mean))
    #map(c->deriv2[c] .= calc_d2fluxdlnlambda2(flux_mean[c],λ[c]),chunk_map )
    for c in chunk_map
         calc_d2fluxdlnlambda2!(deriv2[c],flux_mean[c],λ[c])
    end
    return deriv2
end

function calc_rvs_from_taylor_expansion(spectra::STS; mean::MT = calc_mean_spectrum(spectra),
                deriv::DT = calc_mean_dfluxdlnlambda(spectra), idx::RT = 1:length(mean),
                equal_weight::Bool = true ) where {
                    STS<:AbstractSpectralTimeSeriesCommonWavelengths, T1<:Real, MT<:AbstractVector{T1},
                    T2<:Real, DT<:AbstractVector{T2}, RT<:AbstractUnitRange }
   @assert length(mean) == length(deriv)
   @assert size(spectra.flux,1) == length(mean)
   @assert minimum(idx) >=1
   @assert maximum(idx) <= length(mean)

   if equal_weight # Pixels equally-weighted
      norm = sum(abs2.(deriv[idx]))
      rv = sum((spectra.flux[idx,:].-mean[idx]).*deriv[idx],dims=1).*(speed_of_light_mps/norm)
      # TODO: WARN: Uncertinaties ignores correlations between pixels, particularly problematic when oversample pixels
      σ_rv = sqrt.(sum(spectra.var[idx,:].*abs2.(deriv[idx]),dims=1)) .*(speed_of_light_mps/norm)
  else # Pixels inverse variance weighted
      # TODO:  CHECK/FIX? math for inverse variance weighting.
      # Or maybe the math is right, but it's just a bad idea to have different weighting within one line/chunk
      @info "Warning: I think this is either wrong or a bad idea.  Need to check."
      norm = sum(abs2.(deriv[idx])./spectra.var[idx,:])
      rv = sum((spectra.flux[idx,:].-mean[idx]).*deriv[idx]./spectra.var[idx,:],dims=1).*(speed_of_light_mps/norm)
      σ_rv = sqrt.(sum(abs2.(deriv)./spectra.var[idx,:],dims=1)).*(speed_of_light_mps/norm)
   end
   return (rv=-vec(rv), σ_rv=vec(σ_rv))
end

function calc_rvs_from_taylor_expansion_alt(spectra::STS; mean::MT = calc_mean_spectrum(spectra),
                deriv::DT = calc_mean_dfluxdlnlambda(spectra),
                deriv2::DT = calc_d2fluxdlnlambda2(spectra), idx::RT = 1:length(mean),
                equal_weight::Bool = true ) where {
                    STS<:AbstractSpectralTimeSeriesCommonWavelengths, T1<:Real, MT<:AbstractVector{T1},
                    T2<:Real, DT<:AbstractVector{T2}, RT<:AbstractUnitRange }
   @assert length(mean) == length(deriv)
   @assert size(spectra.flux,1) == length(mean)
   @assert minimum(idx) >=1
   @assert maximum(idx) <= length(mean)

   if equal_weight # Pixels equally-weighted
      #norm = sum(abs2.(deriv[idx]))
      norm = sum(deriv[idx].^2 .- (spectra.flux[idx,:].-mean[idx]).*deriv2[idx],dims=1)
      rv = sum(((spectra.flux[idx,:].-mean[idx]).*deriv[idx]),dims=1)./norm
      rv *= speed_of_light_mps
      # TODO: WARN: Uncertinaties ignores correlations between pixels, particularly problematic when oversample pixels
      σ_rv = sqrt.(sum(spectra.var[idx,:].*abs2.(deriv[idx]),dims=1)./norm)
      σ_rv *= speed_of_light_mps
  else # Pixels inverse variance weighted
      # TODO:  CHECK/FIX? math for inverse variance weighting.
      # Or maybe the math is right, but it's just a bad idea to have different weighting within one line/chunk
      @info "Warning: I think this is either wrong or a bad idea.  Need to check."
      norm = sum(abs2.(deriv[idx]./spectra.var[idx,:]))  # TODO comparable to norm_old, if new norm works, then check updating this also helps
      rv = sum((spectra.flux[idx,:].-mean[idx]).*deriv[idx]./spectra.var[idx,:],dims=1).*(speed_of_light_mps/norm)
      σ_rv = sqrt.(sum(abs2.(deriv)./spectra.var[idx,:],dims=1)).*(speed_of_light_mps/norm)
   end
   return (rv=-vec(rv), σ_rv=vec(σ_rv))
end

function calc_chunk_rvs_from_taylor_expansion(spectra::STS; mean::MT = calc_mean_spectrum(spectra),
                deriv::DT = calc_mean_dfluxdlnlambda(spectra),
                equal_weight::Bool = false ) where {
                    STS<:AbstractSpectralTimeSeriesCommonWavelengths, T1<:Real, MT<:AbstractVector{T1},
                    T2<:Real, DT<:AbstractVector{T2}, RT<:AbstractUnitRange }
   @assert length(mean) == length(deriv)
   @assert size(spectra.flux,1) == length(mean)

   map(idx->calc_rvs_from_taylor_expansion(spectra,mean=mean,deriv=deriv,idx=idx),spectra.chunk_map)
end


text:
using CSV, Tables, Statistics, BenchmarkTools

"""
    read_data(csv_file::String, T::Type{<:AbstractFloat})::Matrix{T}

Read the voltage-current data in `csv_file` into a matrix, where each column is a (V, I) point.
`T` is the data type we want.
"""
function read_data(csv_file::String, T::Type{<:AbstractFloat})::Matrix{T}
    data = CSV.File(csv_file) |> Tables.matrix
    # we want each column as a (V, I) point
    return convert.(T, data')
end

"""
    measure_RMSE_stats(f::Function, n_points; n_trials=30)

Measure the runtime of function `f`. `n_points` is the number of data points (used in computing RMSE).
"""
function measure_RMSE_stats(f::Function, n_points; n_trials=30)
    rmse_list = zeros(n_trials)
    for i = 1:n_trials
        res = f()
        rmse_list[i] = sqrt(res.sse / n_points)
    end
    return Dict(
        :min => minimum(rmse_list),
        :mean => mean(rmse_list),
        :max => maximum(rmse_list),
        :std => std(rmse_list)
    )
end

function benchmark(f::Function; n_trials=30)
    @benchmark f() evals = 1 samples = n_trials seconds = 1000
end

max_stars_repo_name: ShuhuaGao/PV-PPSO

max_issues_repo_name: ShuhuaGao/PV-PPSO

max_forks_repo_name: ShuhuaGao/PV-PPSO


text:
# standard complex number functions

import Base: conj, abs, sqrt, exp, angle, real, imag, log

function real(a::RC)
    if isnan(a)
        return NaN
    end
    if isinf(a)
        return Inf
    end
    return real(a.val)
end

function imag(a::RC)
    if isnan(a)
        return NaN
    end
    if isinf(a)
        return Inf
    end
    return imag(a.val)
end

function conj(a::RC)::RC
    if isinf(a) || isnan(a)
        return a
    end
    return RC(conj(a.val))
end

function abs(a::RC)::Real
    if isinf(a)
        return Inf
    end

    if isnan(a)
        return NaN
    end

    return abs(a.val)
end

function sqrt(a::RC)::RC
    if isnan(a) || isinf(a)
        return a
    end
    return RC(sqrt(a.val))
end

function exp(a::RC)::RC
    if isnan(a) || isinf(a)
        return ComplexNaN
    end
    return RC(exp(a.val))
end

function log(a::RC)::RC
    if isnan(a) || iszero(a) || isinf(a)
        return ComplexNaN
    end
    return RC(log(a.val))
end

function angle(a::RC)::Real
    if isnan(a) || isinf(a)
        return NaN
    end
    return angle(a.val)
end


text:
using Distributed
addprocs()

@everywhere begin
  using StochasticDiffEq, DiffEqProblemLibrary, ParallelDataTransfer, Random, Distributed
  using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
  prob = DiffEqProblemLibrary.SDEProblemLibrary.oval2ModelExample(largeFluctuations=true,useBigs=false)
  Random.seed!(99 + myid())
  prob = remake(prob,tspan=(0.0,1.0))
  println("Solve once to compile.")
  sol = solve(prob,EM(),dt=1/2^(18),adaptive=false,save_everystep=false)
  sol = solve(prob,RKMil(),dt=1/2^(18),adaptive=false,save_everystep=false)
  sol = solve(prob,SRIW1(),dt=1/2^(18),adaptive=false,save_everystep=false)
  sol = solve(prob,SRI(),dt=1/2^(18),adaptive=false,save_everystep=false)
  sol = solve(prob,SOSRI(),dt=1/2^(18),adaptive=false,save_everystep=false)
  sol = solve(prob,SOSRI2(),dt=1/2^(18),adaptive=false,save_everystep=false)
  Int(sol.u[1]!=NaN)
  println("Compilation complete.")
  js = 16:21
  dts = 1.0 ./ 2.0 .^ (js)
  fails = Array{Int}(undef,length(dts),3)
  times = Array{Float64}(undef,length(dts),3)
  numRuns = 10000
end
println("Setup Complete")

## Timing Runs

@everywhere function runAdaptiveSRIW1(i)
  sol = solve(prob,SRIW1(),abstol=2.0^(-13),reltol=2.0^(-7),maxIters=Int(1e11),qmax=1.125,save_everystep=false)
  Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
@everywhere Random.seed!(99 + myid())
adaptiveTime = @elapsed numFails = sum(pmap(runAdaptiveSRIW1,1:numRuns))
println("The number of Adaptive Fails is $numFails. Elapsed time was $adaptiveTime")

## Timing Runs

@everywhere function runAdaptiveSRI(i)
    sol = solve(prob,SRI(error_terms=2),abstol=2.0^(-13),reltol=2.0^(-7),maxIters=Int(1e11),qmax=1.125,save_everystep=false)
  Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
@everywhere Random.seed!(99 + myid())
adaptiveTime = @elapsed numFails = sum(pmap(runAdaptiveSRI,1:numRuns))
println("The number of Adaptive Fails is $numFails. Elapsed time was $adaptiveTime")

## Timing Runs

@everywhere function runAdaptiveSRI(i)
  sol = solve(prob,SRI(),abstol=2.0^(-14),reltol=2.0^(-18),maxIters=Int(1e11),qmax=1.125,save_everystep=false)
  Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
@everywhere Random.seed!(99 + myid())
adaptiveTime = @elapsed numFails = sum(pmap(runAdaptiveSRI,1:numRuns))
println("The number of Adaptive Fails is $numFails. Elapsed time was $adaptiveTime")

## Timing Runs

@everywhere function runAdaptiveSRIOpt1(i)
  sol = solve(prob,SRI(tableau=StochasticDiffEq.constructSRIOpt1()),abstol=2.0^(-7),reltol=2.0^(-4),maxIters=Int(1e11),qmax=1.125,save_everystep=false)
  Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
@everywhere Random.seed!(99 + myid())
adaptiveTime = @elapsed numFails = sum(pmap(runAdaptiveSRIOpt1,1:numRuns))
println("The number of Adaptive Fails is $numFails. Elapsed time was $adaptiveTime")

## Timing Runs

@everywhere function runAdaptiveSRIOpt1(i)
  sol = solve(prob,SOSRI(),abstol=2.0^(-7),reltol=2.0^(-4),maxIters=Int(1e11),qmax=1.125,save_everystep=false)
  Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
@everywhere Random.seed!(99 + myid())
adaptiveTime = @elapsed numFails = sum(pmap(runAdaptiveSRIOpt1,1:numRuns))
println("The number of Adaptive Fails is $numFails. Elapsed time was $adaptiveTime")

## Timing Runs

@everywhere function runAdaptiveSRIOpt1(i)
  sol = solve(prob,SOSRI(),abstol=2.0^(-7),reltol=2.0^(-6),maxIters=Int(1e11),qmax=1.125,save_everystep=false)
  Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
@everywhere Random.seed!(99 + myid())
adaptiveTime = @elapsed numFails = sum(pmap(runAdaptiveSRIOpt1,1:numRuns))
println("The number of Adaptive Fails is $numFails. Elapsed time was $adaptiveTime")

## Timing Runs

@everywhere function runAdaptiveSRIOpt1(i)
  sol = solve(prob,SOSRI(),abstol=2.0^(-12),reltol=2.0^(-15),maxIters=Int(1e11),qmax=1.125,save_everystep=false)
  Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
@everywhere Random.seed!(99 + myid())
adaptiveTime = @elapsed numFails = sum(pmap(runAdaptiveSRIOpt1,1:numRuns))
println("The number of Adaptive Fails is $numFails. Elapsed time was $adaptiveTime")

## Timing Runs

@everywhere function runAdaptiveSRIOpt1(i)
  sol = solve(prob,SOSRI(),abstol=2.0^(-13),reltol=2.0^(-7),maxIters=Int(1e11),qmax=1.125,save_everystep=false)
  Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
@everywhere Random.seed!(99 + myid())
adaptiveTime = @elapsed numFails = sum(pmap(runAdaptiveSRIOpt1,1:numRuns))
println("The number of Adaptive Fails is $numFails. Elapsed time was $adaptiveTime")

## Timing Runs

@everywhere function runAdaptiveSRIOpt1(i)
  sol = solve(prob,SOSRI(),abstol=2.0^(-12),reltol=2.0^(-15),maxIters=Int(1e11),qmax=1.125,save_everystep=false)
  Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
@everywhere Random.seed!(99 + myid())
adaptiveTime = @elapsed numFails = sum(pmap(runAdaptiveSRIOpt1,1:numRuns))
println("The number of Adaptive Fails is $numFails. Elapsed time was $adaptiveTime")

## Timing Runs

@everywhere function runAdaptiveSRIOpt2(i)
    sol = solve(prob,SOSRI2(),abstol=2.0^(-12),reltol=2.0^(-15),maxIters=Int(1e11),qmax=1.125,save_everystep=false)
    Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
@everywhere Random.seed!(99 + myid())
adaptiveTime = @elapsed numFails = sum(pmap(runAdaptiveSRIOpt2,1:numRuns))
println("The number of Adaptive Fails is $numFails. Elapsed time was $adaptiveTime")

## Timing Runs

@everywhere function runAdaptiveSRIOpt2(i)
    sol = solve(prob,SOSRI2(),abstol=2.0^(-13),reltol=2.0^(-11),maxIters=Int(1e11),qmax=1.125,save_everystep=false)
    Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
@everywhere Random.seed!(99 + myid())
adaptiveTime = @elapsed numFails = sum(pmap(runAdaptiveSRIOpt2,1:numRuns))
println("The number of Adaptive Fails is $numFails. Elapsed time was $adaptiveTime")

## Timing Runs

@everywhere function runAdaptiveSRIOpt2(i)
    sol = solve(prob,SOSRI2(),abstol=2.0^(-16),reltol=2.0^(-9),maxIters=Int(1e11),qmax=1.125,save_everystep=false)
    Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
@everywhere Random.seed!(99 + myid())
adaptiveTime = @elapsed numFails = sum(pmap(runAdaptiveSRIOpt2,1:numRuns))
println("The number of Adaptive Fails is $numFails. Elapsed time was $adaptiveTime")

@everywhere function runEM(i,j)
  sol =solve(prob,EM(),dt=dts[j],maxIters=Int(1e11),save_everystep=false,verbose=false)
  Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
for j in eachindex(js)
  println("j = $j")
  sendto(workers(), j=j)
  @everywhere Random.seed!(99 + myid())
  t1 = @elapsed numFails = sum(pmap((i)->runEM(i,j),1:numRuns))
  println("The number of Euler-Maruyama Fails is $numFails. Elapsed time was $t1")
  fails[j,1] = numFails
  times[j,1] = t1
end

@everywhere function runSRI(i,j)
  sol =solve(prob,SRIW1(),dt=dts[j],maxIters=Int(1e11),adaptive=false,save_everystep=false,verbose=false)
  Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
for j in 1:4
  println("j = $j")
  sendto(workers(), j=j)
  @everywhere Random.seed!(99 + myid())
  t2 = @elapsed numFails = sum(pmap((i)->runSRI(i,j),1:numRuns))
  println("The number of Rossler-SRI Fails is $numFails. Elapsed time was $t2")
  fails[j,2] = numFails
  times[j,2] = t2
end

@everywhere js = 17:21
@everywhere dts = 1.0 ./2.0 .^ (js)
@everywhere function runIEM(i,j)
  sol =solve(prob,ImplicitEM(),dt=dts[j],maxIters=Int(1e11),save_everystep=false,verbose=false)
  Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
for j in 1:6
  println("j = $j")
  sendto(workers(), j=j)
  @everywhere Random.seed!(99 + myid())
  t2 = @elapsed numFails = sum(pmap((i)->runIEM(i,j),1:numRuns))
    println("The number of Implicit-EM Fails is $numFails. Elapsed time was $t2")
  fails[j,2] = numFails
  times[j,2] = t2
end

@everywhere js = 17:21
@everywhere dts = 1.0 ./ 2.0 .^(js)
@everywhere function runIRM(i,j)
    sol =solve(prob,ImplicitRKMil(),dt=dts[j],maxIters=Int(1e11),save_everystep=false,verbose=false)
  Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
for j in 1:4
  println("j = $j")
  sendto(workers(), j=j)
  @everywhere Random.seed!(99 + myid())
  t2 = @elapsed numFails = sum(pmap((i)->runIRM(i,j),1:numRuns))
    println("The number of Implicit-RKMil Fails is $numFails. Elapsed time was $t2")
  fails[j,2] = numFails
  times[j,2] = t2
end

@everywhere function runMil(i,j)
  sol =solve(prob,RKMil(),dt=dts[j],maxIters=Int(1e11),save_everystep=false,verbose=false)
  Int(any(isnan,sol[end]) || sol.t[end] != 1)
end
for j in eachindex(js)
  println("j = $j")
  sendto(workers(), j=j)
  @everywhere Random.seed!(99 + myid())
  t3 = @elapsed numFails = sum(pmap((i)->runMil(i,j),1:numRuns))
  println("The number of RK-Milstein Fails is $numFails. Elapsed time was $t3")
  fails[j,3] = numFails
  times[j,3] = t3
end

using Plots
lw = 3
p2 = plot(dts,times,xscale=:log2,yscale=:log2,guidefont=font(16),tickfont=font(14),yguide="Elapsed Time (s)",xguide=L"Chosen $\Delta t$",top_margin=50px,linewidth=lw,lab=["Euler-Maruyama" "RK-Mil" "RosslerSRI"],legendfont=font(14))
plot!(dts,repmat([adaptiveTime],11),linewidth=lw,line=:dash,lab="ESRK+RSwM3",left_margin=75px)
scatter!([2.0^(-20);2.0^(-20);2.0^(-18)],[times[5,1];times[5,2];times[3,3]],markersize=20,c=:red,lab="")
plot(p2,size=(800,800))

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
# Derived from the following cool example by @nixterrimus: (MIT licensed)
# https://github.com/JensRavens/Interstellar/blob/master/Interstellar.playground/Pages/6%20Using%20flatMap%20on%20signals%20of%20signals.xcplaygroundpage/Contents.swift

"""
It's your first day at Space Inc, and you've been given the job to write some code
to refuel an orbiting space station.  There's already been some code written, and
at the center of this code are Signals.

Talking to the space station, takes a lot of time.  Commands need to be send via a
radio tower and then results returned to the radio tower.  These kind of operations
are perfect to be Signals- they return values... eventually.
"""
module SpaceStation

using ReactiveBasics

"""
Get fuel levels for the space station
Returns an Int for the fuel level

In this stub implementation we return a value right away
"""
get_fuel_levels_from_space_station() = Signal(255)

"""
Sends the command to open the fuel port
Returns a Bool representing if the fuel port is open (true) or closed (false)

In this stub implementation we return a value right away
"""
open_fuel_port() = Signal(true)

"""
Sends the command to refuel from a fuel pack
WARNING: Make sure the fuel port is open before you do this!
WARNING: The space station can only contain 1000 fuel units
Returns an Int representing the fuel level for the space station

In this stub implementation we return a value right away
"""
refuel_space_station_from_fuel_pack(fuel_pack_size::Int) = Signal(255 + fuel_pack_size)

"""
We've been tasked with taking this library code and writing a refueling function.
This is a perfect case for `flatmap`.

`flatmap` is a tool that takes a function that returns a Signal and returns its inner
Signal's values.

Let's start by getting the fuel level and opening the opening the fuel port.

The result of running this code is:
The space station now has 755 fuel units available

Notice how the inner blocks return a Signal but the next blocks reference the value that
is returned by the inner Signal.

`flatmap` is a wonderful tool when you are dealing with multiple asynchrous operations that
need to be chained.
"""
function example()
    fuel_pack_size = 500
    s = get_fuel_levels_from_space_station()
    s = flatmap(s) do fuel_level
        if fuel_level < 500
            return open_fuel_port()
        else
            return Signal(false)
        end
    end
    s = flatmap(s) do fuel_port_is_open
        if fuel_port_is_open
            return refuel_space_station_from_fuel_pack(fuel_pack_size)
        else
            return Signal(0)
        end
    end
    map(space_station_fuel ->
        println("The space station now has $space_station_fuel fuel units available"), s)
end

example()

end # module


text:
# Selectors for nodes in `XExpr`s.

#=
select(doc, "pre[lang='julia']")

SelectTag(:pre) & SelectAttributeEq(:lang, "julia")

select(doc, "h1, h2, h3, h4")

SelectTag(:h1) | SelectTag(:h2) | SelectTag(:h3) | SelectTag(:h4)

=#

# Selectors

abstract type Selector end

matches(sel::Selector, x) = false

struct SelectCond <: Selector
    f
end
matches(sel::SelectCond, x) = sel.f(x)

struct SelectAll <: Selector end
matches(::SelectAll, x) = true

struct SelectNode <: Selector end
matches(::SelectNode, ::XExpr) = true
matches(::SelectNode, _) = false

struct SelectTag <: Selector
    tag::Symbol
end

matches(sel::SelectTag, x::XExpr) = sel.tag == x.tag

struct SelectOr{T<:Tuple} <: Selector
    selectors::T
end

Base.:(|)(sel1::Selector, sel2::Selector) = SelectOr((sel1, sel2))
Base.:(|)(selor::SelectOr, sel::Selector) = SelectOr((selor.seletors..., sel))

matches(sel::SelectOr, x) = any(matches(s, x) for s in sel.selectors)

struct SelectAnd{T<:Tuple} <: Selector
    selectors::T
end

Base.:(&)(sel1::Selector, sel2::Selector) = SelectAnd((sel1, sel2))
Base.:(&)(seland::SelectAnd, sel::Selector) = SelectAnd((seland.selectors..., sel))

matches(sel::SelectAnd, x) = all(matches(s, x) for s in sel.selectors)

struct SelectNot <: Selector
    sel::Selector
end
matches(sel::SelectNot, x) = !matches(sel.sel, x)

Base.:(!)(sel::Selector) = SelectNot(sel)

"""
    SelectAttrEq(attr, val)

Select an x-expression with attribute `attr == val`
"""
struct SelectAttrEq{T} <: Selector
    attr::Symbol
    val::T
end

matches(sel::SelectAttrEq, x::XExpr) = get(x.attributes, sel.attr, nothing) == sel.val

struct SelectHasAttr <: Selector
    attr::Symbol
end

matches(sel::SelectHasAttr, x::XExpr) = haskey(attributes(x), sel.attr)

# Traversal

function traverse(doc::XExpr; ordering = PostOrderDFS)
    return (x for x in ordering(doc) if x isa XExpr)
end

function dfstraversal(doc; childfn = xchildren)

end

xchildren(x::XExpr) = children(x)
xchildren(_) = ()
#=

down(tree, i, )
up(tree, state)
=#

# Select

function select(doc::XExpr, sel::Selector; ordering = PreOrderDFS)
    return (x for x in traverse(doc; ordering = ordering) if matches(sel, x))
end

function selectfirst(doc::XExpr, sel::Selector; ordering = PreOrderDFS)
    for x in traverse(doc; ordering = ordering)
        if matches(sel, x)
            return x
        end
    end
    return nothing
end

function Base.map(f, doc::XExpr, sel::Selector)
    return fold(doc, nothing) do x, s
        return matches(sel, x) ? (f(x), nothing) : (x, nothing)
    end |> first
end

#= Tree operations

`selectfirst(doc, sel)`
`replacefirst(doc, sel)`

- traverse tree

`insertbefore(doc, subdoc, sel)`
`insertafter(doc, subdoc, sel)`
`groupchildren(doc, sel)`

=#
function fold(f, doc::XExpr, state)
    children = []
    for child in doc.children
        child, state = fold(f, child, state)
        push!(children, child)
    end
    return f(xexpr(doc.tag, doc.attributes, children), state)
end

function Base.foreach(f, doc::XExpr, sel::Selector; kwargs...)
    for x in traverse(doc; kwargs...)
        if matches(sel, x)
            f(x)
        end
    end
end

function mapfirst(f, doc::XExpr, sel::Selector)
    return fold(doc, false) do x, hasreplaced
        hasreplaced && return (x, true)
        matches(sel, x) ? (f(x), true) : (x, false)
    end |> first
end

replacefirst(doc, x, sel) = mapfirst(_ -> x, doc, sel)

function foldchildren(f, doc::XExpr, state)
    children = []
    for child in doc.children
        if child isa XExpr
            child, state = foldchildren(f, child, state)
        end
        push!(children, child)
    end
    children, state = f(children, state)

    return xexpr(doc.tag, doc.attributes, children...), state
end

fold(f, val, state) = f(val, state)

"""
    insertafter(sel, doc::XExpr, child)

Insert `child` into `doc` after the first match of `sel`.
"""
function insertafter(sel, doc::XExpr, child)
    return foldchildren(doc, false) do children, inserted
        # only insert once
        if inserted
            return children, true
        # insert into children after match
        else
            matchidx = findfirst([matches(sel, c) for c in children])
            if isnothing(matchidx)
                return children, false
            else
                return [children[1:matchidx]..., child, children[matchidx+1:end]...], true
            end
        end
    end |> first
end

function insertfirstchild(sel, doc::XExpr, child)
    return fold(doc, false) do x, inserted
        if !inserted && matches(sel, x)
            return xexpr(x.tag, x.attributes, child, x.children...), true
        else
            return x, inserted
        end
    end |> first
end

function Base.filter(sel::Selector, doc::XExpr)
    return foldchildren(doc, nothing) do children, _
        return [child for child in children if matches(sel, child)], nothing
    end |> first
end

#=
Positions:

- Before(sel)
- After(sel)
- Child(sel, n)

=#

function replacemany(doc, blocks, sel)
    return fold(doc, 1) do x, i
        if matches(sel, x)
            return blocks[i], i + 1
        else
            return x, i
        end
    end |> first
end


text:
function ComplexODEStateAndDerivative(arg0::jdouble, arg1::Vector{JComplex}, arg2::Vector{JComplex})
    return ComplexODEStateAndDerivative((jdouble, Vector{JComplex}, Vector{JComplex}), arg0, arg1, arg2)
end

function ComplexODEStateAndDerivative(arg0::jdouble, arg1::Vector{JComplex}, arg2::Vector{JComplex}, arg3::Vector{Vector{JComplex}}, arg4::Vector{Vector{JComplex}})
    return ComplexODEStateAndDerivative((jdouble, Vector{JComplex}, Vector{JComplex}, Vector{Vector{JComplex}}, Vector{Vector{JComplex}}), arg0, arg1, arg2, arg3, arg4)
end

function get_complete_derivative(obj::ComplexODEStateAndDerivative)
    return jcall(obj, "getCompleteDerivative", Vector{JComplex}, ())
end

function get_primary_derivative(obj::ComplexODEStateAndDerivative)
    return jcall(obj, "getPrimaryDerivative", Vector{JComplex}, ())
end

function get_secondary_derivative(obj::ComplexODEStateAndDerivative, arg0::jint)
    return jcall(obj, "getSecondaryDerivative", Vector{JComplex}, (jint,), arg0)
end


text:
"""
Pi
"""
const DPI = 3.141592653589793

"""
2Pi
"""
const D2PI = 6.283185307179586

"""
Radians to degrees
"""
const DR2D = 57.29577951308232

"""
Degrees to radians
"""
const DD2R = 0.017453292519943295

"""
Radians to arcseconds
"""
const DR2AS = 206264.80624709636

"""
Arcseconds to radians
"""
const DAS2R = 4.84813681109536e-6

"""
Seconds of time to radians
"""
const DS2R = 7.27220521664304e-5

"""
Arcseconds in a full circle
"""
const TURNAS = 1.296e6

"""
Milliarcseconds to radians
"""
const DMAS2R = DAS2R / 1000.0

"""
Length of tropical year B1900 (days)
"""
const DTY = 365.242198781

"""
Seconds per day
"""
const DAYSEC = 86400.0

"""
Days per Julian year
"""
const DJY = 365.25

"""
Days per Julian century
"""
const DJC = 36525.0

"""
Days per Julian millennium
"""
const DJM = 365250.0

"""
Reference epoch (J2000.0), Julian Date
"""
const DJ00 = 2.451545e6

"""
Julian Date of Modified Julian Date zero
"""
const DJM0 = 2.4000005e6

"""
Reference epoch (J2000.0), Modified Julian Date
"""
const DJM00 = 51544.5

"""
1977 Jan 1.0 as MJD
"""
const DJM77 = 43144.0

"""
TT minus TAI (s)
"""
const TTMTAI = 32.184

"""
Astronomical unit (m, IAU 2012)
"""
const DAU = 1.4959787e11

"""
Speed of light (m/s)
"""
const CMPS = 2.99792458e8

"""
Light time for 1 au (s)
"""
const AULT = 499.004782

"""
Speed of light (au per day)
"""
const DC = DAYSEC / AULT

"""
L_G = 1 - d(TT)/d(TCG)
"""
const ELG = 6.969290134e-10

"""
L_B = 1 - d(TDB)/d(TCB)
"""
const ELB = 1.550519768e-8

"""
TDB (s) at TAI 1977/1/1.0
"""
const TDB0 = -6.55e-5

"""
Schwarzschild radius of the Sun (au) =
2 * 1.32712440041e20 / (2.99792458e8)^2 / 1.49597870700e11
"""
const SRS = 1.97412574336e-8

@enum Ellipsoid WGS84 = 1 GRS80 = 2 WGS72 = 3

mutable struct ASTROM
    pmt::Cdouble
    eb::NTuple{3,Cdouble}
    eh::NTuple{3,Cdouble}
    em::Cdouble
    v::NTuple{3,Cdouble}
    bm1::Cdouble
    bpn::NTuple{9,Cdouble}
    along::Cdouble
    phi::Cdouble
    xpl::Cdouble
    ypl::Cdouble
    sphi::Cdouble
    cphi::Cdouble
    diurab::Cdouble
    eral::Cdouble
    refa::Cdouble
    refb::Cdouble
end

struct LDBODY
    bm::Cdouble
    dl::Cdouble
    pv::NTuple{6,Cdouble}
end

struct ERFAException <: Exception
    msg::String
end

Base.showerror(io::IO, ex::ERFAException) = print(io, ex.msg)

export ERFAException


text:
struct DenseConvRBM{Vd,Vc,H,Wd,Wc,K,K2}
    visible_dense::Vd # dense visible part
    visible_conv::Vc # convolutional visible part
    hidden::H # hidden layer
    w_dense::Wd # dense weights
    w_conv::Wc # convolutional weights
    stride::NTuple{K,Int}
    pad::NTuple{K2,Int}
    dilation::NTuple{K,Int}
    groups::Int
    pool::Bool
    """
        DenseConvRBM(visible_dense, visible_conv, hidden, weights_dense, weights_conv)

    Constructs an RBM where a portion of the visible layer is dense and another portion is
    convolutional.
    """
    function DenseConvRBM(
        visible_dense::AbstractLayer,
        visible_conv::AbstractLayer,
        hidden::AbstractLayer,
        w_dense::AbstractArray,
        w_conv::AbstractArray,
        stride::NTuple{K,Int}, pad::NTuple{K2,Int}, dilation::NTuple{K,Int}, groups::Int,
        pool::Bool = false
    ) where {K,K2}
        # dense size checks
        @assert size(w_dense) == (size(visible_dense)..., size(hidden)...)
        # conv size checks
        @assert K == ndims(w_conv) - ndims(visible_conv) - ndims(hidden)
        @assert K2 == 2K
        @assert size(w_conv)[begin:ndims(visible_conv)] == size(visible_conv)
        @assert size(w_conv)[(ndims(visible_conv) + K + 1):end] == size(hidden)
        # types
        V_dens = typeof(visible_dense)
        V_conv = typeof(visible_conv)
        H = typeof(hidden)
        W_dens = typeof(w_dense)
        W_conv = typeof(w_conv)
        # construct
        return new{V_dens, V_conv, H, W_dens, W_conv, K, K2}(
            visible_dense, visible_conv, hidden, w_dense, w_conv,
            stride, pad, dilation, groups, pool
        )
    end
end

function DenseConvRBM(
    visible_dense::AbstractLayer,
    visible_conv::AbstractLayer,
    hidden::AbstractLayer,
    w_dense::AbstractArray,
    w_conv::AbstractArray;
    stride = 1, pad = 0, dilation = 1, groups::Int = 1, pool::Bool = false
)
    K = ndims(w_conv) - ndims(visible_conv) - ndims(hidden)
    kernel_size = size(w_conv)[(ndims(visible_conv) + 1):(ndims(visible_conv) + K)]
    return DenseConvRBM(
        visible_dense, visible_conv, hidden, w_dense, w_conv,
        expand_tuple(Val(K), stride),
        parsepad(kernel_size, dilation, pad),
        expand_tuple(Val(K), dilation),
        groups, pool
    )
end

"""
    parts(denseconvrbm)

Extracts the dense and convolutional parts of a `DenseConvRBM`.
"""
parts(rbm::DenseConvRBM) = (
    dense = RBM(rbm.visible_dense, rbm.hidden, rbm.w_dense),
    conv = ConvRBM(
        rbm.visible_conv, rbm.hidden, rbm.w_conv;
        rbm.stride, rbm.pad, rbm.dilation, rbm.groups, rbm.pool
    )
)

RBMs.hidden(rbm::DenseConvRBM) = rbm.hidden

channel_ndims(rbm::DenseConvRBM) = channel_ndims(parts(rbm).conv)
kernel_ndims(rbm::DenseConvRBM) = kernel_ndims(parts(rbm).conv)
kernel_size(rbm::DenseConvRBM) = kernel_size(parts(rbm).conv)
hsizes(rbm::DenseConvRBM, h::AbstractArray) = hsizes(parts(rbm).conv, h)

input_dims(rbm::DenseConvRBM) = (channel_ndims(rbm) + 1):(ndims(rbm.w_conv) - ndims(hidden(rbm)))
output_dims(rbm::DenseConvRBM) = (ndims(hidden(rbm)) + 1):(ndims(hidden(rbm)) + kernel_ndims(rbm))

output_size(rbm::DenseConvRBM, v_conv::AbstractArray) = output_size(parts(rbm).conv, v_conv)

function RBMs.batch_size(
    rbm::DenseConvRBM, v_dense::AbstractArray, v_conv::AbstractArray, h::AbstractArray
)
    bsz_dense = batch_size(rbm.visible_dense, v_dense)
    bsz_conv = batch_size(parts(rbm).conv, v_conv, h)
    if isempty(bsz_dense)
        return bsz_conv
    elseif isempty(bsz_conv)
        return bsz_dense
    else
        @assert bsz_dense == bsz_conv
        return bsz_dense
    end
end

function RBMs.batch_size(rbm::DenseConvRBM, v_dense::AbstractArray, v_conv::AbstractArray)
    bsz_dense = batch_size(rbm.visible_dense, v_dense)
    bsz_conv = batch_size(parts(rbm).conv, v_conv)
    if isempty(bsz_dense)
        return bsz_conv
    elseif isempty(bsz_conv)
        return bsz_dense
    else
        @assert bsz_dense == bsz_conv
        return bsz_dense
    end
end

function RBMs.inputs_v_to_h(rbm::DenseConvRBM, v_dense::AbstractArray, v_conv::AbstractArray)
    I_dense = inputs_v_dense_to_h(rbm, v_dense)
    I_conv = inputs_v_conv_to_h(rbm, v_conv)
    return I_dense .+ I_conv
end

function inputs_v_conv_to_h(rbm::DenseConvRBM, v_conv::AbstractArray)
    return inputs_v_to_h(parts(rbm).conv, v_conv)
end

function inputs_v_dense_to_h(rbm::DenseConvRBM, v_dense::AbstractArray)
    I_dense = inputs_v_to_h(parts(rbm).dense, v_dense)
    bsz = batch_size(rbm.visible_dense, v_dense)
    @assert size(I_dense) == (size(rbm.hidden)..., bsz...)
    return reshape(I_dense, size(rbm.hidden)..., one.(kernel_size(rbm))..., bsz...)
end

function inputs_h_to_v_dense(rbm::DenseConvRBM, h::AbstractArray)
    I_dense = inputs_h_to_v(parts(rbm).dense, h)
    hsz = hsizes(rbm, h)
    @assert size(I_dense) == (size(rbm.visible_dense)..., hsz.output_size..., hsz.batch_size...)
    kdims = ndims(rbm.visible_dense) .+ (1:kernel_ndims(rbm))
    return reshape(sum(I_dense; dims=kdims), size(rbm.visible_dense)..., hsz.batch_size...)
end

inputs_h_to_v_conv(rbm::DenseConvRBM, h::AbstractArray) = inputs_h_to_v(parts(rbm).conv, h)

function RBMs.energy(
    rbm::DenseConvRBM, v_dense::AbstractArray, v_conv::AbstractArray, h::AbstractArray
)
    E_v_dense = energy(rbm.visible_dense, v_dense)
    E_v_conv = energy(rbm.visible_conv, v_conv)
    E_h = energy(rbm.hidden, h)
    E_w = interaction_energy(rbm, v_dense, v_conv, h)

    E_v_conv_reduced = reshape_maybe(
        sum(E_v_conv; dims=1:kernel_ndims(rbm)),
        vsizes(parts(rbm).conv, v_conv).batch_size
    )
    E_h_reduced = reshape_maybe(
        sum(E_h; dims=1:kernel_ndims(rbm)),
        hsizes(rbm, h).batch_size
    )

    return E_v_dense .+ E_v_conv_reduced .+ E_h_reduced .+ E_w
end

function RBMs.interaction_energy(
    rbm::DenseConvRBM, v_dense::AbstractArray, v_conv::AbstractArray, h::AbstractArray
)
    I = RBMs.inputs_v_to_h(rbm, v_dense, v_conv)
    @assert size(I) == (
        size(hidden(rbm))..., hsizes(rbm, h).output_size..., batch_size(rbm, v_dense, v_conv)...
    )
    E = -sum(I .* h; dims=1:(ndims(hidden(rbm)) + kernel_ndims(rbm)))
    return reshape_maybe(E, batch_size(rbm, v_dense, v_conv, h))
end

function RBMs.free_energy(rbm::DenseConvRBM, v_dense::AbstractArray, v_conv::AbstractArray; β::Real = 1)
    if rbm.pool
        return free_energy_pooled(rbm, v_dense, v_conv; β)
    else
        return free_energy_nopool(rbm, v_dense, v_conv; β)
    end
end

function RBMs.sample_h_from_v(
    rbm::DenseConvRBM, v_dense::AbstractArray, v_conv::AbstractArray; β::Real = true
)
    if rbm.pool
        return sample_h_from_v_pooled(rbm, v_dense, v_conv; β)
    else
        return sample_h_from_v_nopool(rbm, v_dense, v_conv; β)
    end
end

function sample_v_dense_from_h(rbm::DenseConvRBM, h::AbstractArray; β::Real = true)
    inputs = inputs_h_to_v_dense(rbm, h)
    return transfer_sample(rbm.visible_dense, inputs; β)
end

function sample_v_conv_from_h(rbm::DenseConvRBM, h::AbstractArray; β::Real = true)
    inputs = inputs_h_to_v_conv(rbm, h)
    return transfer_sample(rbm.visible_conv, inputs; β)
end

function RBMs.sample_v_from_v(
    rbm::DenseConvRBM, v_dense::AbstractArray, v_conv::AbstractArray;
    β::Real = true, steps::Int = 1
)
    for _ in 1:steps
        h = sample_h_from_v(rbm, v_dense, v_conv; β)
        v_dense = oftype(v_dense, sample_v_dense_from_h(rbm, h; β))
        v_conv = oftype(v_conv, sample_v_conv_from_h(rbm, h; β))
    end
    return (v_dense = v_dense, v_conv = v_conv)
end

function RBMs.sample_h_from_h(
    rbm::DenseConvRBM, h::AbstractArray; β::Real = true, steps::Int = 1
)
    for _ in 1:steps
        h = oftype(h, RBMs.sample_h_from_h_once(rbm, h; β))
    end
    return h
end

function RBMs.sample_h_from_h_once(rbm::DenseConvRBM, h::AbstractArray; β::Real = true)
    v_dense = sample_v_dense_from_h(rbm, h; β)
    v_conv = sample_v_conv_from_h(rbm, h; β)
    return RBMs.sample_h_from_v(rbm, v_dense, v_conv; β)
end


text:
"""
    (f_q::GP)(X::AbstractVector)

A GP evaluated at `X` is a finite-dimensional GP (i.e. a multivariate Normal).
"""
(f_q::GP)(X::AbstractVector) = GP(f_q, X)

"""
    (f_q::BlockGP)(X::BlockData)

Index the `c`th component `AbstractGP` of `f_q` at `X[c]`.

    (f_q::BlockGP)(X::AbstractVector)

Index each `c`th component `AbstractGP` of `f_q` at `X`.
"""
(f_q::BlockGP)(X::BlockData) = BlockGP(map((f, x)->f(x), f_q.fs, blocks(X)))
(f_q::BlockGP)(X::AbstractVector) = BlockGP([f(X) for f in f_q.fs])

μ_p′(f_q::GP, X::AVM) = finite(mean(f_q), X)
k_p′(f_q::GP, X::AVM) = finite(kernel(f_q), X)
k_p′p(f_q::GP, X::AVM, f_p::GP) = lhsfinite(kernel(f_q, f_p), X)
k_pp′(f_p::GP, f_q::GP, X′::AVM) = rhsfinite(kernel(f_p, f_q), X′)


text:
using BinaryBuilder

name = "dSFMT"
version = v"2.2.4"

# Collection of sources required to build dSFMT
sources = [
    GitSource("https://github.com/MersenneTwister-Lab/dSFMT.git",
              "5a02974a257ae74dbab12bde3ef6a5ffc0cfbbc2"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/dSFMT

FLAGS=(
    -O3 -finline-functions -fomit-frame-pointer -fno-strict-aliasing -Wmissing-prototypes -Wall -std=c99
    -DNDEBUG -DDSFMT_MEXP=19937
    -fPIC -shared -DDSFMT_SHLIB -DDSFMT_DO_NOT_USE_OLD_NAMES
)

if [[ ${target} == x86_64* ]]; then
    FLAGS+=(-msse2 -DHAVE_SSE2)
fi

mkdir -p "${libdir}"
${CC} ${FLAGS[@]} ${CFLAGS} ${CPPFLAGS} ${LDFLAGS} -o "${libdir}/libdSFMT.${dlext}" dSFMT.c
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(;experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libdSFMT", :libdSFMT),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
#Include KG_Module
include("../src/KG_Module.jl")
include("../src/KG_Plot_Module.jl")

#Run slices
slices_observer("./data/animation_centra/binary_bh_rotating_radius_1.5_alpha_10.0_Rorbit_3.0.h5")

#generate_animation_heatmap_makie2("./data/animation_centra/binary_bh_rotating_radius_1.5_alpha_10.0_Rorbit_3.0.h5" , "my_test_2.mp4")


text:
precompile(rmvn, (Int, Vector{Float64}, Matrix{Float64}))
precompile(rmvn, (Float64, Vector{Float64}, Matrix{Float64}))
precompile(rmvn, (Int, Matrix{Float64}))
precompile(rmvn, (Float64, Matrix{Float64}))

precompile(cor, (Matrix{Float64}, Pearson))
precompile(cor, (Matrix{Float64}, Spearman))
precompile(cor, (Matrix{Float64}, Kendall))

precompile(cor_fast, (Matrix{Float64}, Pearson))
precompile(cor_fast, (Matrix{Float64}, Spearman))
precompile(cor_fast, (Matrix{Float64}, Kendall))

precompile(cor_randPD, (Float64, Int, Int))
precompile(cor_randPSD, (Float64, Int, Int))

precompile(cor_fastPD!, (Matrix{Float64}, Float64))
precompile(cor_fastPD, (Matrix{Float64}, Float64))

precompile(cor_nearPD, (Matrix{Float64},))
precompile(cor_nearPD, (Matrix{Float64}, Float64, Float64))

precompile(pearson_bounds, (DUD, DUD))
precompile(pearson_bounds, (CUD, DUD))
precompile(pearson_bounds, (DUD, CUD))
precompile(pearson_bounds, (CUD, CUD))

precompile(pearson_match, (Float64, CUD, CUD))
precompile(pearson_match, (Float64, CUD, DUD))
precompile(pearson_match, (Float64, DUD, CUD))
precompile(pearson_match, (Float64, DUD, DUD))

precompile(GSDist, (DUD,))
precompile(GSDist, (DUD, Float64))
precompile(GSDist, (DUD, Float64, Int64, Function, Float64))


text:
# we write a function that caculates the product of two matrices A and B
function mult(A,B)

	# we find the sizes of A and B
	aRows, aCols = size(A)
	bRows, bCols = size(B)

	# we check if the matrix sizes are consistent with multiplication
	if aCols != bRows
		error("matrix sizes don't agree")
	end

	# we initialize C to a zero matrix of the appropriate size
	C = zeros(aRows, bCols)

	# we loop over all pairs of rows and columns
	for i = 1:aRows
		for j = 1:bCols
			# and for each pair we store their dot pruduct in C
			for k = 1:aCols
				C[i,j] += A[i,k] * B[k,j];
			end
		end
	end

	# we return C as the answer
	return C
end

max_stars_repo_path: 1a/Matrix Multiplication/mult.jl

max_issues_repo_path: 1a/Matrix Multiplication/mult.jl

max_forks_repo_path: 1a/Matrix Multiplication/mult.jl


text:
#=
this file includes functionality for estimating linear scaling regions and
defines the `generalized_dim` function.
=#

export linear_region, linear_regions, estimate_boxsizes, linreg
#####################################################################################
# Functions and methods to deduce linear scaling regions
#####################################################################################
using Statistics
using Statistics: covm, varm
# The following function comes from a version in StatsBase that is now deleted
# StatsBase is copyrighted under the MIT License with
# Copyright (c) 2012-2016: Dahua Lin, Simon Byrne, Andreas Noack, Douglas Bates,
# John Myles White, Simon Kornblith, and other contributors.
"""
    linreg(x, y) -> a, b
Perform a linear regression to find the best coefficients so that the curve:
`z = a + b*x` has the least squared error with `y`.
"""
function linreg(x::AbstractVector, y::AbstractVector)
    # Least squares given
    # Y = a + b*X
    # where
    # b = cov(X, Y)/var(X)
    # a = mean(Y) - b*mean(X)
    if size(x) != size(y)
        throw(DimensionMismatch("x has size $(size(x)) and y has size $(size(y)), " *
            "but these must be the same size"))
    end
    mx = Statistics.mean(x)
    my = Statistics.mean(y)
    # don't need to worry about the scaling (n vs n - 1)
    # since they cancel in the ratio
    b = covm(x, mx, y, my)/varm(x, mx)
    a = my - b*mx
    return a, b
end

slope(x, y) = linreg(x, y)[2]

"""
    linear_regions(x, y; dxi::Int = 1, tol = 0.25) -> (lrs, tangents)
Identify regions where the curve `y(x)` is linear, by scanning the
`x`-axis every `dxi` indices sequentially
(e.g. at `x[1] to x[5], x[5] to x[10], x[10] to x[15]` and so on if `dxi=5`).

If the slope (calculated via linear regression) of a region of width `dxi` is
approximatelly equal to that of the previous region,
within tolerance `tol`,
then these two regions belong to the same linear region.

Return the indices of `x` that correspond to linear regions, `lrs`,
and the _correct_ `tangents` at each region
(obtained via a second linear regression at each accumulated region).
"""
function linear_regions(
        x::AbstractVector, y::AbstractVector;
        method = :sequential, dxi::Int = method == :overlap ? 3 : 1, tol = 0.25,
    )
    @assert length(x) == length(y)
    return if method == :overlap
        linear_regions_overlap(x, y, dxi, tol)
    elseif method == :sequential
        linear_regions_sequential(x, y, dxi, tol)
    end
end

function linear_regions_sequential(x, y, dxi, tol)
    maxit = length(x) ÷ dxi

    tangents = Float64[slope(view(x, 1:max(dxi, 2)), view(y, 1:max(dxi, 2)))]

    prevtang = tangents[1]
    lrs = Int[1] #start of first linear region is always 1
    lastk = 1

    # Start loop over all partitions of `x` into `dxi` intervals:
    for k in 1:maxit-1
        tang = slope(view(x, k*dxi:(k+1)*dxi), view(y, k*dxi:(k+1)*dxi))
        if isapprox(tang, prevtang, rtol=tol, atol = 0)
            # Tanget is similar with initial previous one (based on tolerance)
            continue
        else
            # Tangent is not similar.
            # Push new tangent for a new linear region
            push!(tangents, tang)

            # Set the START of a new linear region
            # which is also the END of the previous linear region
            push!(lrs, k*dxi)
            lastk = k
        end

        # Set new previous tangent (only if it was not the same as current)
        prevtang = tang
    end
    push!(lrs, length(x))
    # create new tangents that do have linear regression weighted
    tangents = Float64[]
    for i in 1:length(lrs)-1
        push!(tangents, linreg(view(x, lrs[i]:lrs[i+1]), view(y ,lrs[i]:lrs[i+1]))[2])
    end
    return lrs, tangents
end

"""
    linear_region(x, y; kwargs...) -> ((ind1, ind2), slope)
Call [`linear_regions`](@ref) and identify and return the largest linear region
and its slope. The region starts and stops at `x[ind1:ind2]`.

The keywords `dxi, tol` are propagated as-is to [`linear_regions`](@ref).
The keyword `ignore_saturation = true` ignores saturation that (sometimes) happens
at the start and end of the curve `y(x)`, where the curve flattens.
The keyword `sat = 0.01` decides what saturation is (while `abs(y[i]-y[i+1])<sat` we
are in a saturation regime).

The keyword `warning = true` prints a warning if the linear region is less than 1/3
of the available x-axis.
"""
function linear_region(x::AbstractVector, y::AbstractVector;
    dxi::Int = 1, tol::Real = 0.2, ignore_saturation = true, warning = true, sat = 0.01)

    isat = 0
    if ignore_saturation
        j = findfirst(i -> abs(y[i] - y[i-1]) > sat, length(y):-1:2)
        if !isnothing(j)
            i = (length(y):-1:2)[j]
            x, y = x[1:i], y[1:i]
        end
        k = findfirst(i -> abs(y[i+1] - y[i]) > sat, 1:length(y)-1)
        if !isnothing(k)
            x, y = x[k:end], y[k:end]
            isat = k-1
        end
    end

    lrs, tangents = linear_regions(x,y; dxi, tol)
    # Find biggest linear region:
    j = findmax(diff(lrs))[2]
    if lrs[j+1] - lrs[j] ≤ length(x)÷3 && warning
        @warn "Found linear region spans less than a 3rd of the available x-axis "*
              "and might imply inaccurate slope or insufficient data. "*
              "Recommended: plot `x` vs `y`."
    end
    return (lrs[j] + isat, lrs[j+1] + isat), tangents[j]
end

#####################################################################################
# Autotomatic estimation for proper `ε` from a Dataset
#####################################################################################
"""
    estimate_boxsizes(A::Dataset; kwargs...) → εs
Return `k` exponentially spaced values: `εs = base .^ range(lower + w, upper + z; length = k)`,
that are a good estimate for sizes ε that are used in calculating a [Fractal Dimension](@ref).
It is strongly recommended to [`standardize`](@ref) input dataset `A` before using this
function.

Let `d₋` be the minimum pair-wise distance in `A` and `d₊` the average total length of `A`
along each of the dimensions of `A`.
Then `lower = log(base, d₋)` and `upper = log(base, d₊)`.
Because by default `w=1, z=-1`, the returned sizes are an order of mangitude
larger than the minimum distance, and an order of magnitude smaller than the maximum
distance.

## Keywords
* `w = 1, z = -1, k = 20` : as explained above.
* `base = MathConstants.e` : the base used in the `log` function.
* `warning = true`: Print some warnings for bad estimates.
* `autoexpand = true`: If the final estimated range does not cover at least 2 orders of
  magnitude, it is automatically expanded by setting `w -= we` and `z -= ze`.
  You can set different default values to the keywords `we = w, ze = z`.
"""
function estimate_boxsizes(
        A::AbstractDataset;
        k::Int = 20, z = -1, w = 1, base = MathConstants.e,
        warning = true, autoexpand = true, ze = z, we = w
    )

    mi, ma = minmaxima(A)
    max_d = mean(ma - mi)
    min_d, _ = minimum_pairwise_distance(A)
    if min_d == 0 && warning
        @warn(
        "Minimum distance in the dataset is zero! Probably because of having data "*
        "with low resolution, or duplicate data points. Setting to `d₊/base^4` for now.")
        min_d = max_d/(base^4)
    end

    lower = log(base, min_d)
    upper = log(base, max_d)

    if lower ≥ upper
        error("`lower ≥ upper`. There must be something fundamentally wrong with dataset.")
    elseif lower+w ≥ upper+z && warning
        @warn(
        "Automatic boxsize determination was inappropriate: `lower+w` was found ≥ than "*
        "`upper+z`. Returning `base .^ range(lower, upper; length = k)`. "*
        "Please adjust keywords or provide a bigger dataset.")
        εs = float(base) .^ range(lower, upper; length = k)
    elseif abs(upper+z - (lower+w)) < 2 && autoexpand
        if warning
            @warn(
            "Boxsize limits do not differ by at least 2 orders of magnitude. "*
            "Setting `w-=$(we)` and `z+=$(ze)`, please adjust keywords `w, z` otherwise.")
        end
        εs = float(base) .^ range(lower+w-we, upper+z-ze; length = k)
    else
        εs = float(base) .^ range(lower+w, upper+z; length = k)
    end
    return εs
end

"""
    minimum_pairwise_distance(A::Dataset, metric = Euclidean())
Return `min_d, min_pair`: the minimum pairwise distance
of all points in the dataset, and the corresponding point pair.
"""
function minimum_pairwise_distance(A::AbstractDataset, metric = Euclidean())
    tree = KDTree(A, metric)
    min_d = eltype(A[1])(Inf)
    max_d = -min_d
    min_pair = max_pair = (0, 0)
    theiler = Theiler(0)
    for i in 1:length(A)
        inds, dists = Neighborhood.knn(tree, A[i], 1, theiler(i); sortds=false)
        ind, dist = inds[1], dists[1]
        if dist < min_d
            min_d = dist
            min_pair = (i, ind)
        end
    end
    return min_d, min_pair
end


text:
#Function to deal with positioners
struct POSRobotBody
    #positioner body is theta arm and head is for phi arm
    #Any thing stored in two element array will always asume body first and head second
    max_seg::Integer #maximum number of segements this object can store

    xy::Array{Float64,1} #To store the xy position of the positioner
    arms::Array{Float64,1} #To store the theta.phi arm of the positioner
    arms2_lim::Array{Float64,1} #To store the square of minimum and maximum limit of the robot arm

    nseg::Array{Integer,1} #Stores the number of segment in body and head
    restricted_angels::Array{Bool,1} #Store if any angular restriction is needed

    body_minmax::Array{Float64,1} #Stores the minimum and maximum length of the verted for body
    body_orig::Array{Float64, 2} #Stores the orginal polygon of body
    head_minmax::Array{Float64,1} #Stores the minimum and maximum length of the verted for head
    head_orig::Array{Float64, 2} #Stores the original polygon of head

    axis_orig::Array{Float64,1} #To store the orignal axis

    #To get the theta limit and phi limit: Angles are stored in radian
    body_theta_limit::Array{Float64,1} #Store the limit of angles for body MIN_T,MAX_T,OFFSET_T
    head_theta_limit::Array{Float64,1} #Store the limit of angles for head MIN_P,MAX_P,OFFSET_P

    body_on_target::Array{Float64, 2} #Stores the body polygon on target
    head_on_target::Array{Float64, 2} #Stores the head polygon on target
    axis_on_target::Array{Float64,1} #This goes throgh the same transformation as head to keep tract of its axis

    function  POSRobotBody(maximum_seg::Integer)
        new(maximum_seg,zeros(2),zeros(2),zeros(2),zeros(2),
            zeros(Bool,2),zeros(2),zeros(maximum_seg,2),zeros(2),
            zeros(maximum_seg,2),zeros(2),zeros(3),zeros(3),
            zeros(maximum_seg,2),zeros(maximum_seg,2),zeros(2))
    end
end

"""
Initialize the pos:POSRobotBody with the position at index in fp_dic
"""
function load_positioner!(pos::POSRobotBody,index,fp_dic,exc_dic)

    #get the exclusion names
    excl_name=fp_dic["EXCLUSION"][index]
    excl_body_name="$(excl_name)-theta-segments"
    excl_head_name="$(excl_name)-phi-segments"

    #Assign number of segments
    pos.nseg[1]=size(exc_dic[excl_body_name],1)
    pos.nseg[2]=size(exc_dic[excl_head_name],1)

    #Assign polygon for body
    Assign_polygon_with_maxmin_vertex!(pos.max_seg,pos.nseg[1],
        exc_dic[excl_body_name],pos.body_orig,pos.body_minmax)

    #Assign polygon for head
    Assign_polygon_with_maxmin_vertex!(pos.max_seg,pos.nseg[2],
        exc_dic[excl_head_name],pos.head_orig,pos.head_minmax)

    #get the angular limit of the robots body
    pos.body_theta_limit[1]=deg2rad(fp_dic["MIN_T"][index])
    pos.body_theta_limit[2]=deg2rad(fp_dic["MAX_T"][index])
    pos.body_theta_limit[3]=deg2rad(fp_dic["OFFSET_T"][index])

    #get the angular limit of the robots head
    pos.head_theta_limit[1]=deg2rad(fp_dic["MIN_P"][index])
    pos.head_theta_limit[2]=deg2rad(fp_dic["MAX_P"][index])
    pos.head_theta_limit[3]=deg2rad(fp_dic["OFFSET_P"][index])

    #Allocate xy
    pos.xy[1]=fp_dic["OFFSET_X"][index]
    pos.xy[2]=fp_dic["OFFSET_Y"][index]

    #Assign arms
    pos.arms[1]=fp_dic["LENGTH_R1"][index]
    pos.arms[2]=fp_dic["LENGTH_R2"][index]

    #arms limits
    pos.arms2_lim[1]=(pos.arms[1]-pos.arms[2])^2 #minimum
    pos.arms2_lim[2]=(pos.arms[1]+pos.arms[2])^2 #maximum

end

function Assign_polygon_with_maxmin_vertex!(max_seg,nseg,polygon,pos_poly,minmax)
    cent=[mean(polygon[:,1]),mean(polygon[:,2])]
    min_vert2=100000
    max_vert2=0

    for ii in 1:nseg
        pos_poly[ii,1]=polygon[ii,1]
        pos_poly[ii,2]=polygon[ii,2]
        rad2=(polygon[ii,1]-cent[1])^2 + (polygon[ii,2]-cent[2])^2
        if(max_vert2<rad2)
            max_vert2=rad2
        end
        if(min_vert2>rad2)
            min_vert2=rad2
        end
    end

    for ii in nseg+1:max_seg
        pos_poly[ii,1]= 0
        pos_poly[ii,2]= 0
    end

    minmax[1]=sqrt(min_vert2)
    minmax[2]=sqrt(max_vert2)
end

"""
Applies the tanslation on poly_in by x,y and return result in poly_out
"""
function poly_translate!(poly_in,poly_out,x,y)
    for ii in 1:size(poly_in,1)
        poly_out[ii,1]=poly_in[ii,1]+x
        poly_out[ii,2]=poly_in[ii,2]+y
    end
end

function point_translate!(pt_in,pt_out,x,y)
    pt_out[1]= pt_in[1] + x
    pt_out[2]= pt_in[2] + y
end

"""Rotates the points in poly around the axis
by cos(theta),sin(theta)

diff_v: optional to use pre-allocated array for intermediate step for speed
cos_sin: optional to use pre-allocated array for intermediate step for speed
cos_sin_sum: optional to use pre-allocated array for intermediate step for speed
"""
function poly_rotate_axis!(poly_in,poly_out,cstheta,axis;
    diff_v=zeros(Float64,2),cos_sin=zeros(Float64,2),cos_sin_sum=zeros(Float64,2))
    for ii in 1:size(poly_in,1)
        point_rotate_axis!(view(poly_in,ii,:),cstheta,axis,view(poly_out,ii,:);
         diff_v=diff_v,cos_sin=cos_sin,cos_sin_sum=cos_sin_sum)
    end
end

"""
Apply a rotation to a point along the axis
point_in: input point x,y coordinate
point_out: x,y co-ordinate after rotation
cstheta: [cos(theta, sin(theta)] for rotation by theta
axis: [0,0] : center about which rotation is needed
diff_v: optional to use pre-allocated array for intermediate step for speed
cos_sin: optional to use pre-allocated array for intermediate step for speed
cos_sin_sum: optional to use pre-allocated array for intermediate step for speed
"""
function point_rotate_axis!(point_in,cstheta,axis,point_out;
        diff_v=zeros(Float64,2),cos_sin=zeros(Float64,2),cos_sin_sum=zeros(Float64,2))
    diff_v .= point_in .- axis
    if((diff_v[1]==0) & (diff_v[2]==0))
        return point_in[1],point_in[2]
    end

    rad=sqrt(sum(diff_v.^2))

    #cos_sin of difference vector
    cos_sin .= diff_v ./rad

    #sum of angles cos(A+B) and sin (A+B) rules
    cos_sin_sum .= [cos_sin[1]*cstheta[1] - cos_sin[2]*cstheta[2],
                 cos_sin[1]*cstheta[2] + cos_sin[2]*cstheta[1]  ]

    #rotation and translation
    point_out[1]=axis[1]+(rad*cos_sin_sum[1])
    point_out[2]=axis[2]+(rad*cos_sin_sum[2])
end

"""moves the positioner to a theta phi position
Adapted from desihub/fiberassign
diff_v: optional to use pre-allocated array for intermediate step for speed
cos_sin: optional to use pre-allocated array for intermediate step for speed
cos_sin_sum: optional to use pre-allocated array for intermediate step for speed
csthetha: optional to use pre-allocated array for intermediate step for speed
csphi: optional to use pre-allocated array for intermediate step for speed
"""
function move_positioner_to_theta_phi!(pos::POSRobotBody,theta,phi;
    diff_v=zeros(Float64,2),cos_sin=zeros(Float64,2),cos_sin_sum=zeros(Float64,2),
    cstheta=zeros(Float64,2),csphi=zeros(Float64,2))

    cstheta .= [cos(theta),sin(theta)]
    csphi .= [cos(phi),sin(phi)]

    #// Move the phi polygon into the fully extended position along the X axis.
    #moving to LENGHT_R1,0
    nT=pos.nseg[1]
    nP=pos.nseg[2]

    poly_translate!(view(pos.head_orig,1:nP,:),view(pos.head_on_target,1:nP,:),pos.arms[1],0)
    #translate the axis
    point_translate!(pos.axis_orig,pos.axis_on_target,pos.arms[1],0)

    #@show pos.head_on_target[1:pos.nseg[2],2] .== (pos.head_orig[1:pos.nseg[2],2] .+ 0)
    #@show "1: ",pos.axis_on_target,pos.arms[1]

    #// Rotate fully extended positioner an angle of theta about the origin.
    poly_rotate_axis!(view(pos.body_orig,1:nT,:),view(pos.body_on_target,1:nT,:),cstheta,[0.0,0.0],
        diff_v=diff_v,cos_sin=cos_sin,cos_sin_sum=cos_sin_sum)
    poly_rotate_axis!(view(pos.head_on_target,1:nP,:),view(pos.head_on_target,1:nP,:),cstheta,[0.0,0.0],
        diff_v=diff_v,cos_sin=cos_sin,cos_sin_sum=cos_sin_sum)
    #rotate the axis
    point_rotate_axis!(pos.axis_on_target,cstheta,[0.0,0.0],pos.axis_on_target,
        diff_v=diff_v,cos_sin=cos_sin,cos_sin_sum=cos_sin_sum)

    #@show "2: ",pos.axis_on_target,cstheta

    #// Rotate just the phi arm an angle phi about the theta arm center.
    #shpphi.rotation(csphi);
    poly_rotate_axis!(view(pos.head_on_target,1:nP,:),view(pos.head_on_target,1:nP,:),csphi,pos.axis_on_target,
        diff_v=diff_v,cos_sin=cos_sin,cos_sin_sum=cos_sin_sum)

    #// Translate the whole positioner to the center.
    poly_translate!(view(pos.head_on_target,1:nP,:),view(pos.head_on_target,1:nP,:),pos.xy[1],pos.xy[2])
    poly_translate!(view(pos.body_on_target,1:nT,:),view(pos.body_on_target,1:nT,:),pos.xy[1],pos.xy[2])

    #Not needed just for plotting
    #point_translate!(pos.axis_on_target,pos.axis_on_target,pos.xy[1],pos.xy[2])
end

"""Determines if a positioner is valid of not
fp_dic : dictionary contaning the status of focal plane
pp: index of positioner considered
pos_type: purpose of validity testing
           assignment: whether this positioner can be assigned to a target
           neighbour_collision: whether this positioner should be considered for neighbour collision
"""
function valid_positionter(fp_dic,pp;pos_type="assignment")
    if(pos_type=="assignment")
        return ((fp_dic["DEVICE_TYPE"][pp]=="POS") & (fp_dic["LENGTH_R1"][pp]!=0.0))
    elseif(pos_type=="neighbour_collision")
        return ((fp_dic["DEVICE_TYPE"][pp]=="POS") & (fp_dic["LENGTH_R1"][pp]!=0.0))
    end
end

"""
loads the positioner negibours given pindex and zones of the central positioner
pos_ngb_dic: a dictionary of structure to hold the negibours
"""
function load_positioner_neighbours!(pos_cur,pindex,fp_bound,fp_dic,exc_dic,pos_ngb_dic)
    #This zone is already loaded so nothing to be done
    if(pos_ngb_dic["index_central_positioner"]==pindex)
        #already this positioner neighbours are loaded
        return
    else
        zone_x=fp_dic["zone_x"][pindex]
        zone_y=fp_dic["zone_y"][pindex]
        zone_xy=get_zonexy(zone_x,zone_y,fp_bound)
        pos_ngb_dic["index_central_positioner"]=pindex
    end

    #@show zone_x,zone_y,zone_xy
    count_pos=1
    for zx in maximum([zone_x-1,1]):minimum([zone_x+1,fp_bound.nx])
        for zy in maximum([zone_y-1,1]):minimum([zone_y+1,fp_bound.ny])
            zone_xy=get_zonexy(zx,zy,fp_bound)

            #@show zx,zy,zone_xy

            if(!(zone_xy in fp_dic["zone_map_keys"]))
                continue
            end
            imin=fp_dic["zone_map"][zone_xy][1]
            imax=fp_dic["zone_map"][zone_xy][2]
            #@show zx,zy,zone_xy,targets_dic["zone_map"][zone_xy]
            #@show targets_dic["zone_x"][imin:imax]
            count_this=imax-imin
            pos_ngb_dic["index"][count_pos:count_pos+count_this] .= collect(imin:imax)
            count_pos=count_pos+count_this+1
        end
    end

    #set everything else to 0
    if(count_pos<=size(pos_ngb_dic["index"],1))
        pos_ngb_dic["index"][count_pos:end] .= 0
    end
    #Due to Julia indexing convetion 1 need to be subtracte
    count_pos=count_pos-1

    #Determine the maximum distance to look for within positioner
    #This is defined as the maximum arams length + the maximum vertex of head
    #The above is conservative and specific choice will not matter
    lookup_rad2=(pos_cur.arms[1]+pos_cur.arms[2]+pos_cur.head_minmax[1])^2
    lookup_rad2=4*lookup_rad2 #factor of 4 is for twice quare the distance
    #@show lookup_rad2,pos_cur.arms[1],pos_cur.arms[2],pos_cur.head_minmax[2]

    #scan for positioners within look_up radius
    count_ngb=0
    for ii in pos_ngb_dic["index"][1:count_pos]
        if((ii==pindex) || (!valid_positionter(fp_dic,ii;pos_type="neighbour_collision")))#To not count self positioner and invalid positioner
            continue
        end
        rad2=(fp_dic["OFFSET_X"][ii]-fp_dic["OFFSET_X"][pindex])^2
        rad2 += (fp_dic["OFFSET_Y"][ii]-fp_dic["OFFSET_Y"][pindex])^2
        #@show ii,rad2,lookup_rad2,rad2<lookup_rad2

        if(rad2<lookup_rad2)
            count_ngb+=1
            #load this positioner info
            load_positioner!(pos_ngb_dic["pos"][count_ngb],ii,fp_dic,exc_dic)
            pos_ngb_dic["pos_index"][count_ngb]=ii
        end
    end
    pos_ngb_dic["count_ngb"]=count_ngb
    #@show count_ngb

    if(false)#To make plots
        tind=pos_ngb_dic["index"][1:count_pos]
        plot([fp_dic["OFFSET_X"][pindex]],[fp_dic["OFFSET_Y"][pindex]],
            seriestype=:scatter,color=:green,markersize=10)
        for ii in 1:count_ngb
            xy=pos_ngb_dic["pos"][ii].xy
            plot!([xy[1]],[xy[2]],seriestype=:scatter,color=:red,markersize=10)
        end

        plot!(fp_dic["OFFSET_X"][tind],fp_dic["OFFSET_Y"][tind],seriestype=:scatter,color=:black) |> display
    end

end

"""Compare the moved positioner positioner to check if either the body or head-boy is collided
body-body collision is not considered as it is not possible in the focal plane
pos1: first positioner of type POSRobotBody
pos2: second positioner of type POSRobotBody
"""
function positioner_collided(pos1,pos2)

    head1=view(pos1.head_on_target,1:pos1.nseg[2],:)
    head2=view(pos2.head_on_target,1:pos2.nseg[2],:)

    #check for the head1-head2 collision
    if(PolyIntersect_boundingbox(head1,head2))
        if(PolygonsIntersect(head1,head2))
            return true
        end
    end

    body1=view(pos1.body_on_target,1:pos1.nseg[1],:)
    body2=view(pos2.body_on_target,1:pos2.nseg[1],:)

    #head1-body2 collision
    if(PolyIntersect_boundingbox(head1,body2))
        if(PolygonsIntersect(head1,body2))
            return true
        end
    end

    #head2-body1 collision
    if(PolyIntersect_boundingbox(head2,body1))
        if(PolygonsIntersect(head2,body1))
            return true
        end
    end

    return false
end

max_stars_repo_name: shadaba/MockFBA

max_issues_repo_name: shadaba/MockFBA

max_forks_repo_name: shadaba/MockFBA


text:
# ---
# title: 606. Construct String from Binary Tree
# id: problem606
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String, Tree
# link: <https://leetcode.com/problems/construct-string-from-binary-tree/description/>
# hidden: true
# ---
#
# You need to construct a string consists of parenthesis and integers from a
# binary tree with the preorder traversing way.
#
# The null node needs to be represented by empty parenthesis pair "()". And you
# need to omit all the empty parenthesis pairs that don't affect the one-to-one
# mapping relationship between the string and the original binary tree.
#
# **Example 1:**
#
#
#
#     Input: Binary tree: [1,2,3,4]
#            1
#          /   \
#         2     3
#        /
#       4
#
#     Output: "1(2(4))(3)"
#
#     Explanation: Originallay it needs to be "1(2(4)())(3()())",
#     but you need to omit all the unnecessary empty parenthesis pairs.
#     And it will be "1(2(4))(3)".
#
#
# **Example 2:**
#
#
#
#     Input: Binary tree: [1,2,3,null,4]
#            1
#          /   \
#         2     3
#          \
#           4
#
#     Output: "1(2()(4))(3)"
#
#     Explanation: Almost the same as the first example,
#     except we can't omit the first parenthesis pair to break the one-to-one mapping relationship between the input and the output.
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    map_impact_parameters(m::AbstractMetricParams{T}, u, α, β)

Map impact parameters `α` and `β` to a velocity vector at some position `u` in the given metric `m`.
"""
function map_impact_parameters(
    m::AbstractMetricParams{T},
    u::AbstractVector{T},
    α::N,
    β::N
) where {T,N<:Number}
    [0.0, alpha_beta_to_vel(m, u, α, β)...]
end

function map_impact_parameters(
    m::AbstractMetricParams{T},
    u::SVector{S,T},
    α::N,
    β::N
) where {S,T,N<:Number}
    SVector{S,T}(0.0, alpha_beta_to_vel(m, u, α, β)...)
end

function map_impact_parameters(
    m::AbstractMetricParams{T},
    u,
    α::AbstractVector{P},
    β::AbstractVector{P}
) where {T,P}
    curried(_α, _β) = map_impact_parameters(m, u, _α, _β)
    curried.(α, β)
end

function alpha_beta_to_vel(m::AbstractMetricParams{T}, u, α, β) where {T}
    reg = u[2]^2
    -1.0, -β / reg, -α / reg
end

export map_impact_parameters


text:
# This script contains some code for symbolic manipulation of linear least-squares problems. It can be used to derive or verify matrix expressions associated with problems on the form Ax = y
using SymPy
using SymPy: expand, coeff
SymPy.Sym(name, size::Tuple;kwargs...) = Sym[symbols("$name$i$j";kwargs...) for i = 1:size[1], j=1:size[2]]
skewcoords(R) = [R[3,2];R[1,3];R[2,1]]
twistcoords(xi) = [xi[1:3, 4]; skewcoords(xi[1:3, 1:3])]
skew(s) = [0 -s[3] s[2];s[3] 0 -s[1]; -s[2] s[1] 0]
trinv(T) = [T[1:3,1:3]' -T[1:3,1:3]'*T[1:3,4];0 0 0 1]

## Calib NAXP ==============================================================================
# N' * H*Pf = norm(N)^2

F = Sym("F",(3,3),real=true)
d = Sym("d",(3,1),real=true)
N =  Sym("N",(3,1),real=true)
P =  Sym("P_",(3,1),real=true)

RA = Sym("RA",(3,3), real=true)
TA = Sym("TA",(3,1), real=true)
RX = Sym("RX",(3,3), real=true)
TX = Sym("TX",(3,1), real=true)
X = [RX TX;0 0 0 1]
A = [RA TA;0 0 0 1]
Ps = [P;1]
N0 = [N;0]

xh = symbols("xh", real=true)
yh = symbols("yh", real=true)

# H = F'[[1 0;0 1;0 0], -d];

w = X[1:3,[1, 2, 4]]
w = w[:]
Prb = A*(X*Ps)
eq1 = N'*Prb[1:3]
eq = eq1[:] # Vector of interesting equations,
# terms of order zeros should be placed on RHS

AA     = Array{Sym}(length(eq), length(w))
BB     = Array{Sym}(length(eq), 1)
AAz    = falses(size(AA))

for i = 1:length(eq)
    for j = 1:length(w)
        AA[i,j] = coeff(expand(eq[i]),w[j]) # Find all coefficients of parameter j in eq i
    end
end

println("Trying to find expressions common in all terms")
change,res  = cse(AA[:]);
res         = reshape(res,size(AA));

byhand = (reshape(repmat((N'A[1:3,1:3]),3,1)',9,1) .*[reshape(repmat(P[1:2],1,3)',6,1);1;1;1])
byhand2 = [N'A[1:3,1:3]*Ps[1] N'A[1:3,1:3]*Ps[2] N'A[1:3,1:3]]

@assert all(simplify.(byhand) .== simplify.(byhand2[:]))
@assert all(simplify.(byhand) .== simplify.(AA[:]))


text:
module Revise

# We use a code structure where all `using` and `import`
# statements in the package that load anything other than
# a Julia base or stdlib package are located in this file here.
# Nothing else should appear in this file here, apart from
# the `include("packagedef.jl")` statement, which loads what
# we would normally consider the bulk of the package code.
# This somewhat unusual structure is in place to support
# the VS Code extension integration.

using OrderedCollections, CodeTracking, JuliaInterpreter, LoweredCodeUtils

using CodeTracking: PkgFiles, basedir, srcfiles, line_is_decl, basepath
using JuliaInterpreter: whichtt, is_doc_expr, step_expr!, finish_and_return!, get_return,
                        @lookup, moduleof, scopeof, pc_expr, is_quotenode_egal,
                        linetable, codelocs, LineTypes, is_GotoIfNot, isassign, isidentical
using LoweredCodeUtils: next_or_nothing!, trackedheads, structheads, callee_matches

include("packagedef.jl")

end # module