text:
module JuliaInterpreter

using Base.Meta
import Base: +, -, convert, isless
using Core: CodeInfo, SSAValue, SlotNumber, TypeMapEntry, SimpleVector, LineInfoNode, GotoNode, Slot,
            GeneratedFunctionStub, MethodInstance, NewvarNode, TypeName

using UUIDs
# The following are for circumventing #28, memcpy invalid instruction error,
# in Base and stdlib
using Random.DSFMT
using InteractiveUtils
using CodeTracking

export @interpret, Compiled, JuliaStackFrame,
       Breakpoints, breakpoint, @breakpoint, breakpoints, enable, disable, remove

module CompiledCalls
# This module is for handling intrinsics that must be compiled (llvmcall)
end

"""
`Compiled` is a trait indicating that any `:call` expressions should be evaluated
using Julia's normal compiled-code evaluation. The alternative is to pass `stack=JuliaStackFrame[]`,
which will cause all calls to be evaluated via the interpreter.
"""
struct Compiled end
Base.similar(::Compiled, sz) = Compiled()  # to support similar(stack, 0)

"""
    JuliaProgramCounter(next_stmt::Int)

A wrapper specifying the index of the next statement in the lowered code to be executed.
"""
struct JuliaProgramCounter
    next_stmt::Int
end
+(x::JuliaProgramCounter, y::Integer) = JuliaProgramCounter(x.next_stmt+y)
-(x::JuliaProgramCounter, y::Integer) = JuliaProgramCounter(x.next_stmt-y)
convert(::Type{Int}, pc::JuliaProgramCounter) = pc.next_stmt
isless(x::JuliaProgramCounter, y::Integer) = isless(x.next_stmt, y)

Base.show(io::IO, pc::JuliaProgramCounter) = print(io, "JuliaProgramCounter(", pc.next_stmt, ')')

# Breakpoint support
truecondition(frame) = true
falsecondition(frame) = false
const break_on_error = Ref(false)

"""
    BreakpointState(isactive=true, condition=JuliaInterpreter.truecondition)

`BreakpointState` represents a breakpoint at a particular statement in
a `JuliaFrameCode`. `isactive` indicates whether the breakpoint is currently
[`enable`](@ref)d or [`disable`](@ref)d. `condition` is a function that accepts
a single `JuliaStackFrame`, and `condition(frame)` must return either
`true` or `false`. Execution will stop at a breakpoint only if `isactive`
and `condition(frame)` both evaluate as `true`. The default `condition` always
returns `true`.

To create these objects, see [`breakpoint`](@ref).
"""
struct BreakpointState
    isactive::Bool
    condition::Function
end
BreakpointState(isactive::Bool) = BreakpointState(isactive, truecondition)
BreakpointState() = BreakpointState(true)

# A type used transiently in renumbering CodeInfo SSAValues (to distinguish a new SSAValue from an old one)
struct NewSSAValue
    id::Int
end

"""
`JuliaFrameCode` holds static information about a method or toplevel code.
One `JuliaFrameCode` can be shared by many `JuliaStackFrame` calling frames.

Important fields:
- `scope`: the `Method` or `Module` in which this frame is to be evaluated
- `code`: the `CodeInfo` object storing (optimized) lowered code
- `methodtables`: a vector, each entry potentially stores a "local method table" for the corresponding
  `:call` expression in `code` (undefined entries correspond to statements that do not
  contain `:call` expressions)
- `used`: a `BitSet` storing the list of SSAValues that get referenced by later statements.
"""
struct JuliaFrameCode
    scope::Union{Method,Module}
    code::CodeInfo
    methodtables::Vector{Union{Compiled,TypeMapEntry}} # line-by-line method tables for generic-function :call Exprs
    breakpoints::Vector{BreakpointState}
    used::BitSet
    wrapper::Bool
    generator::Bool
    # Display options
    fullpath::Bool
end

function JuliaFrameCode(frame::JuliaFrameCode; wrapper = frame.wrapper, generator=frame.generator, fullpath=frame.fullpath)
    JuliaFrameCode(frame.scope, frame.code, frame.methodtables, frame.breakpoints, frame.used,
                   wrapper, generator, fullpath)
end

function JuliaFrameCode(scope, code::CodeInfo; wrapper=false, generator=false, fullpath=true, optimize=true)
    if optimize
        code, methodtables = optimize!(copy_codeinfo(code), moduleof(scope))
    else
        code = copy_codeinfo(code)
        methodtables = Vector{Union{Compiled,TypeMapEntry}}(undef, length(code.code))
    end
    breakpoints = Vector{BreakpointState}(undef, length(code.code))
    used = find_used(code)
    return JuliaFrameCode(scope, code, methodtables, breakpoints, used, wrapper, generator, fullpath)
end

"""
`JuliaStackFrame` represents the current execution state in a particular call frame.

Important fields:
- `code`: the [`JuliaFrameCode`](@ref) for this frame
- `locals`: a vector containing the input arguments and named local variables for this frame.
  The indexing corresponds to the names in `frame.code.code.slotnames`.
- `ssavalues`: a vector containing the
  [Static Single Assignment](https://en.wikipedia.org/wiki/Static_single_assignment_form)
  values produced at the current state of execution
- `sparams`: the static type parameters, e.g., for `f(x::Vector{T}) where T` this would store
  the value of `T` given the particular input `x`.
- `pc`: the [`JuliaProgramCounter`](@ref) that typically represents the current position
  during execution. However, note that some internal functions instead maintain the `pc`
  as a local variable, and only update the frame's `pc` when pushing a frame on the stack.
"""
struct JuliaStackFrame
    code::JuliaFrameCode
    locals::Vector{Union{Nothing,Some{Any}}}
    ssavalues::Vector{Any}
    sparams::Vector{Any}
    exception_frames::Vector{Int}
    last_exception::Base.RefValue{Any}
    pc::Base.RefValue{JuliaProgramCounter}
    # A vector from names to the slotnumber of that name
    # for which a reference was last encountered.
    last_reference::Dict{Symbol,Int}
    callargs::Vector{Any}  # a temporary for processing arguments of :call exprs
end

function JuliaStackFrame(framecode::JuliaFrameCode, frame::JuliaStackFrame, pc::JuliaProgramCounter; kwargs...)
    pcref = frame.pc
    pcref[] = pc
    if !isempty(kwargs)
        framecode = JuliaFrameCode(framecode; kwargs...)
    end
    JuliaStackFrame(framecode, frame.locals,
                    frame.ssavalues, frame.sparams,
                    frame.exception_frames, frame.last_exception,
                    pcref, frame.last_reference, frame.callargs)
end

JuliaStackFrame(frame::JuliaStackFrame, pc::JuliaProgramCounter; kwargs...) =
    JuliaStackFrame(frame.code, frame, pc; kwargs...)

"""
`Variable` is a struct representing a variable with an asigned value.
By calling the function `locals`[@ref] on a `JuliaStackFrame`[@ref] a
`Vector` of `Variable`'s is returned.

Important fields:
- `value::Any`: the value of the local variable
- `name::Symbol`: the name of the variable as given in the source code
- `isparam::Bool`: if the variable is a type parameter, for example `T` in `f(x::T) where {T} = x` .
"""
struct Variable
    value::Any
    name::Symbol
    isparam::Bool
end
Base.show(io::IO, var::Variable) = (print(io, var.name, " = "); show(io,var.value))
Base.isequal(var1::Variable, var2::Variable) =
    var1.value == var2.value && var1.name == var2.name && var1.isparam == var2.isparam

"""
    local_variables = locals(frame::JuliaStackFrame)::Vector{Variable}

Return the local variables as a vector of `Variable`[@ref].
"""
function locals(frame::JuliaStackFrame)
    vars = Variable[]
    for (sym, idx) in frame.last_reference
        push!(vars, Variable(something(frame.locals[idx]), sym, false))
    end
    reverse!(vars)
    if frame.code.scope isa Method
        syms = sparam_syms(frame.code.scope)
        for i in 1:length(syms)
            push!(vars, Variable(frame.sparams[i], syms[i], true))
        end
    end
    return vars
end

"""
`framedict[method]` returns the `JuliaFrameCode` for `method`. For `@generated` methods,
see [`genframedict`](@ref).
"""
const framedict = Dict{Method,JuliaFrameCode}()                # essentially a method table for lowered code

"""
`genframedict[(method,argtypes)]` returns the `JuliaFrameCode` for a `@generated` method `method`,
for the particular argument types `argtypes`.

The framecodes stored in `genframedict` are for the code returned by the generator
(i.e, what will run when you call the method on particular argument types);
for the generator itself, its framecode would be stored in [`framedict`](@ref).
"""
const genframedict = Dict{Tuple{Method,Type},JuliaFrameCode}() # the same for @generated functions

"""
`meth ∈ compiled_methods` indicates that `meth` should be run using [`Compiled`](@ref)
rather than recursed into via the interpreter.
"""
const compiled_methods = Set{Method}()

const junk = JuliaStackFrame[]      # to allow re-use of allocated memory (this is otherwise a bottleneck)

const empty_svec = Core.svec()

include("localmethtable.jl")
include("interpret.jl")
include("builtins-julia$(Int(VERSION.major)).$(Int(VERSION.minor)).jl")

function show_stackloc(io::IO, stack, frame, pc=frame.pc[])
    indent = ""
    for f in stack
        println(io, indent, f.code.scope)
        indent *= "  "
    end
    println(io, indent, frame.code.scope, ", pc = ", convert(Int, pc))
end
function show_stackloc(io::IO, ::Compiled, frame, pc)
    println(io, "No stack, ::Compiled")
    println(io, frame.code.scope, ", pc = ", convert(Int, pc))
end
show_stackloc(stack, frame, pc=frame.pc[]) = show_stackloc(stderr, stack, frame, pc)

function moduleof(x)
    if isa(x, JuliaStackFrame)
        x = x.code.scope
    end
    return _moduleof(x)
end
_moduleof(scope::Method) = scope.module
_moduleof(scope::Module) = scope
Base.nameof(frame) = isa(frame.code.scope, Method) ? frame.code.scope.name : nameof(frame.code.scope)

is_loc_meta(expr, kind) = isexpr(expr, :meta) && length(expr.args) >= 1 && expr.args[1] === kind

"""
    isglobalref(g, mod, name)

Tests whether `g` is equal to `GlobalRef(mod, name)`.
"""
isglobalref(g, mod, name) = isa(g, GlobalRef) && g.mod === mod && g.name == name

function to_function(x)
    if isa(x, GlobalRef)
        getfield(x.mod, x.name)
    else
        x
    end
end

_Typeof(x) = isa(x,Type) ? Type{x} : typeof(x)

function is_function_def(ex)
    (isexpr(ex, :(=)) && isexpr(ex.args[1], :call)) ||
    isexpr(ex,:function)
end

function is_generated(meth)
    isdefined(meth, :generator)
end

function sparam_syms(meth::Method)
    s = Symbol[]
    sig = meth.sig
    while sig isa UnionAll
        push!(s, Symbol(sig.var.name))
        sig = sig.body
    end
    return s
end

function namedtuple(kwargs)
    names, types, vals = Symbol[], [], []
    for pr in kwargs
        if isa(pr, Expr)
            push!(names, pr.args[1])
            val = pr.args[2]
            push!(types, typeof(val))
            push!(vals, val)
        elseif isa(pr, Pair)
            push!(names, pr.first)
            val = pr.second
            push!(types, typeof(val))
            push!(vals, val)
        else
            error("unhandled entry type ", typeof(pr))
        end
    end
    return NamedTuple{(names...,), Tuple{types...}}(vals)
end

"""
    frun, allargs = prepare_args(fcall, fargs, kwargs)

Prepare the complete argument sequence for a call to `fcall`. `fargs = [fcall, args...]` is a list
containing both `fcall` (the `#self#` slot in lowered code) and the positional
arguments supplied to `fcall`. `kwargs` is a list of keyword arguments, supplied either as
list of expressions `:(kwname=kwval)` or pairs `:kwname=>kwval`.

For non-keyword methods, `frun === fcall`, but for methods with keywords `frun` will be the
keyword-sorter function for `fcall`.

# Example

```jldoctest; setup=(using JuliaInterpreter; empty!(JuliaInterpreter.junk))
julia> mymethod(x) = 1
mymethod (generic function with 1 method)

julia> mymethod(x, y; verbose=false) = nothing
mymethod (generic function with 2 methods)

julia> JuliaInterpreter.prepare_args(mymethod, [mymethod, 15], ())
(mymethod, Any[mymethod, 15])

julia> JuliaInterpreter.prepare_args(mymethod, [mymethod, 1, 2], [:verbose=>true])
(getfield( Symbol("#kw##mymethod"))(), Any[#kw##mymethod(), (verbose = true,), mymethod, 1, 2])
```
"""
function prepare_args(@nospecialize(f), allargs, kwargs)
    if !isempty(kwargs)
        f = Core.kwfunc(f)
        allargs = [f,namedtuple(kwargs),allargs...]
    elseif f === Core._apply
        f = to_function(allargs[2])
        allargs = Base.append_any((allargs[2],), allargs[3:end]...)
    end
    return f, allargs
end

function prepare_framecode(method::Method, argtypes; enter_generated=false)
    sig = method.sig
    isa(method, TypeMapEntry) && (method = method.func)
    if method.module == Core.Compiler || method.module == Base.Threads || method ∈ compiled_methods
        return Compiled()
    end
    # Get static parameters
    (ti, lenv::SimpleVector) = ccall(:jl_type_intersection_with_env, Any, (Any, Any),
                        argtypes, sig)::SimpleVector
    enter_generated &= is_generated(method)
    if is_generated(method) && !enter_generated
        framecode = get(genframedict, (method, argtypes), nothing)
    else
        framecode = get(framedict, method, nothing)
    end
    if framecode === nothing
        if is_generated(method) && !enter_generated
            # If we're stepping into a staged function, we need to use
            # the specialization, rather than stepping through the
            # unspecialized method.
            code = Core.Compiler.get_staged(Core.Compiler.code_for_method(method, argtypes, lenv, typemax(UInt), false))
            code === nothing && return nothing
            generator = false
        else
            if is_generated(method)
                code = get_source(method.generator)
                generator = true
            else
                code = get_source(method)
                generator = false
            end
        end
        framecode = JuliaFrameCode(method, code; generator=generator)
        if is_generated(method) && !enter_generated
            genframedict[(method, argtypes)] = framecode
        else
            framedict[method] = framecode
        end
    end
    return framecode, lenv
end

function get_framecode(method)
    framecode = get(framedict, method, nothing)
    if framecode === nothing
        code = get_source(method)
        framecode = JuliaFrameCode(method, code; generator=false)
        framedict[method] = framecode
    end
    return framecode
end

function whichtt(tt)
    m = ccall(:jl_gf_invoke_lookup, Any, (Any, UInt), tt, typemax(UInt))
    m === nothing && return nothing
    return m.func::Method
end

"""
    framecode, frameargs, lenv, argtypes = prepare_call(f, allargs; enter_generated=false)

Prepare all the information needed to execute lowered code for `f` given arguments `allargs`.
`f` and `allargs` are the outputs of [`prepare_args`](@ref).
For `@generated` methods, set `enter_generated=true` if you want to extract the lowered code
of the generator itself.

On return `framecode` is the [`JuliaFrameCode`](@ref) of the method.
`frameargs` contains the actual arguments needed for executing this frame (for generators,
this will be the types of `allargs`);
`lenv` is the "environment", i.e., the static parameters for `f` given `allargs`.
`argtypes` is the `Tuple`-type for this specific call (equivalent to the signature of the `MethodInstance`).

# Example

```jldoctest; setup=(using JuliaInterpreter; empty!(JuliaInterpreter.junk))
julia> mymethod(x::Vector{T}) where T = 1
mymethod (generic function with 1 method)

julia> framecode, frameargs, lenv, argtypes = JuliaInterpreter.prepare_call(mymethod, [mymethod, [1.0,2.0]]);

julia> framecode
JuliaInterpreter.JuliaFrameCode(mymethod(x::Array{T,1}) where T in Main at none:1, CodeInfo(
1 ─     return 1
), Union{Compiled, TypeMapEntry}[#undef], JuliaInterpreter.BreakpointState[#undef], BitSet([]), false, false, true)

julia> frameargs
2-element Array{Any,1}:
 mymethod
 [1.0, 2.0]

julia> lenv
svec(Float64)

julia> argtypes
Tuple{typeof(mymethod),Array{Float64,1}}
```
"""
function prepare_call(@nospecialize(f), allargs; enter_generated = false)
    # Can happen for thunks created by generated functions
    if isa(f, Core.Builtin) || isa(f, Core.IntrinsicFunction)
        return nothing
    elseif any(x->isa(x, Type) && x <: Vararg, allargs)
        return nothing  # https://github.com/JuliaLang/julia/issues/30995
    end
    argtypes = Tuple{map(_Typeof,allargs)...}
    method = whichtt(argtypes)
    if method === nothing
        # Call it to generate the exact error
        f(allargs[2:end]...)
    end
    args = allargs
    ret = prepare_framecode(method, argtypes; enter_generated=enter_generated)
    # Exceptional returns
    if ret === nothing
        # The generator threw an error. Let's generate the same error by calling it.
        f(allargs[2:end]...)
    end
    isa(ret, Compiled) && return ret
    # Typical return
    framecode, lenv = ret
    if is_generated(method) && enter_generated
        args = Any[_Typeof(a) for a in args]
    end
    return framecode, args, lenv, argtypes
end

"""
    frame = prepare_thunk(mod::Module, expr::Expr)

Prepare `expr` for evaluation in `mod`. `expr` should be a "straightforward" expression,
one that does not require special top-level handling (see [`JuliaInterpreter.split_expressions`](@ref)).
"""
function prepare_thunk(mod::Module, thunk::Expr, recursive=false)
    if isexpr(thunk, :thunk)
        framecode = JuliaFrameCode(mod, thunk.args[1])
    elseif isexpr(thunk, :error) || isexpr(thunk, :incomplete)
        error("lowering returned an error, ", thunk)
    elseif recursive
        thunk = Meta.lower(mod, Expr(:block, nothing, thunk))
        framecode = JuliaFrameCode(mod, thunk.args[1])
    else
        lwr = Meta.lower(mod, thunk)
        isa(lwr, Expr) && return prepare_thunk(mod, lwr, true)
        return nothing
    end
    return prepare_locals(framecode, [])
end
prepare_thunk((mod, ex)::Tuple{Module,Expr}) = prepare_thunk(mod, ex)

"""
    modexs, docexprs = split_expressions(mod::Module, expr::Expr; extract_docexprs=false)

Break `expr` into a list `modexs` of sequential blocks. This is often needed when `expr`
needs to be evaluated at top level.

`modexs[i]` is a `(mod::Module, ex::Expr)` tuple, where `ex` is to be evaluated in `mod`.

# Toplevel evaluation

For code that defines new structs, new methods, or new macros, it can be important to evaluate
these expressions carefully:

    stack = JuliaStackFrame[]
    for modex in modexs    # or use `for (mod, ex) in modexs` to split the tuple
        frame = JuliaInterpreter.prepare_thunk(modex)
        while true
            JuliaInterpreter.through_methoddef_or_done!(stack, frame) === nothing && break
        end
    end

The `while` loop here deserves some explanation. Occasionally, a frame may define new methods
(e.g., anonymous or local functions) and then call those methods. In such cases, running
the entire frame as a single block (e.g., with [`JuliaInterpreter.finish_and_return!`](@ref)
can trigger "method is too new..." errors. Instead, the approach above runs each frame,
but returns to the caller after any new method is defined. When this loop is running at
top level (e.g., in the REPL), this allows the world age to update and thus avoid
"method is too new..." errors.

Putting the above nested loop inside a function defeats its purpose, because inside a
compiled function the world age will not update. If necessary, use the following strategy:

    Core.eval(somemodule, Expr(:toplevel, quote
        body
    ))

where `body` contains the nested loop, plus any preparatory statements required to make the
necessary variables available at top level in `somemodule`.
"""
function split_expressions(mod::Module, expr::Expr; filename=nothing, kwargs...)
    modexs = Tuple{Module,Expr}[]
    docexprs = Dict{Module,Vector{Expr}}()
    if filename === nothing
        # On Julia 1.2+, the first line of a :toplevel expr may contain line info
        if length(expr.args) >= 1 && isa(expr.args[1], LineNumberNode)
            filename = expr.args[1].file
        else
            filename="toplevel"
        end
    end
    return split_expressions!(modexs, docexprs, mod, expr; filename=filename, kwargs...)
end

"""
    isdocexpr(ex)

Test whether expression `ex` is a `@doc` expression.
"""
function isdocexpr(ex)
    docsym = Symbol("@doc")
    if isexpr(ex, :macrocall)
        a = ex.args[1]
        isglobalref(a, Core, docsym) && return true
        isa(a, Symbol) && a == docsym && return true
        if isexpr(a, :.)
            mod, name = a.args[1], a.args[2]
            return mod === :Core && isa(name, QuoteNode) && name.value == docsym
        end
    end
    return false
end

split_expressions!(modexs, docexprs, mod::Module, ex::Expr; kwargs...) =
    split_expressions!(modexs, docexprs, Expr(:block), mod, ex; kwargs...)

function split_expressions!(modexs, docexprs, lex::Expr, mod::Module, ex::Expr; extract_docexprs=false, filename="toplevel")
    # lex is the expression we'll lower; it will accumulate LineNumberNodes and a
    # single top-level expression. We split blocks, module defs, etc.
    if ex.head == :toplevel || ex.head == :block
        split_expressions!(modexs, docexprs, lex, mod, ex.args; extract_docexprs=extract_docexprs, filename=filename)
    elseif ex.head == :module
        newname = ex.args[2]::Symbol
        if isdefined(mod, newname)
            newmod = getfield(mod, newname)
        else
            if (id = Base.identify_package(mod, String(newname))) !== nothing
                newmod = Base.root_module(id)
            else
                newmod = Core.eval(mod, :(module $newname end))
            end
        end
        split_expressions!(modexs, docexprs, lex, newmod, ex.args[3]; extract_docexprs=extract_docexprs, filename=filename)
    elseif extract_docexprs && isdocexpr(ex)
        docexs = get(docexprs, mod, nothing)
        if docexs === nothing
            docexs = docexprs[mod] = Expr[]
        end
        push!(docexs, ex)
        body = ex.args[4]
        if isa(body, Expr) && body.head != :call
            split_expressions!(modexs, docexprs, lex, mod, body; extract_docexprs=extract_docexprs, filename=filename)
        end
    else
        if isempty(lex.args)
            push!(modexs, (mod, copy(ex)))
        else
            push!(lex.args, ex)
            push!(modexs, (mod, copy(lex)))
            empty!(lex.args)
        end
    end
    return modexs, docexprs
end

function split_expressions!(frames, docexprs, lex, mod::Module, args::Vector{Any}; filename="toplevel", kwargs...)
    for a in args
        if isa(a, Expr)
            split_expressions!(frames, docexprs, lex, mod, a; filename=filename, kwargs...)
        elseif isa(a, LineNumberNode)
            if a.file === nothing  # happens with toplevel expressions on Julia 1.2
                push!(lex.args, LineNumberNode(a.line, Symbol(filename)))
            else
                push!(lex.args, a)
            end
        else
            push!(lex.args, a)
        end
    end
end

"""
    framecode, frameargs, lenv, argtypes = determine_method_for_expr(expr; enter_generated = false)

Prepare all the information needed to execute a particular `:call` expression `expr`.
For example, try `JuliaInterpreter.determine_method_for_expr(:($(sum)([1,2])))`.
See [`JuliaInterpreter.prepare_call`](@ref) for information about the outputs.
"""
function determine_method_for_expr(expr; enter_generated = false)
    f = to_function(expr.args[1])
    allargs = expr.args
    # Extract keyword args
    local kwargs = Expr(:parameters)
    if length(allargs) > 1 && isexpr(allargs[2], :parameters)
        kwargs = splice!(allargs, 2)
    end
    f, allargs = prepare_args(f, allargs, kwargs.args)
    return prepare_call(f, allargs; enter_generated=enter_generated)
end

function get_source(meth)
    if isa(meth.source, Array{UInt8,1})
        return ccall(:jl_uncompress_ast, Any, (Any, Any), meth, meth.source)
    else
        return meth.source
    end
end

function get_source(g::GeneratedFunctionStub)
    b = g(g.argnames...)
    b isa CodeInfo && return b
    return eval(b)
end

function copy_codeinfo(code::CodeInfo)
    @static if VERSION < v"1.1.0-DEV.762"
        newcode = ccall(:jl_new_struct_uninit, Any, (Any,), CodeInfo)::CodeInfo
        for (i, name) in enumerate(fieldnames(CodeInfo))
            if isdefined(code, name)
                val = getfield(code, name)
                ccall(:jl_set_nth_field, Cvoid, (Any, Csize_t, Any), newcode, i-1, val===nothing || isa(val, Union{Type, Method}) ? val : copy(val))
            end
        end
        return newcode
    else
        # Inline this when support for VERSION above is dropped
        return copy(code)
    end
end

const calllike = Set([:call, :foreigncall])

function extract_inner_call!(stmt, idx, once::Bool=false)
    isa(stmt, Expr) || return nothing
    (stmt.head == :toplevel || stmt.head == :thunk) && return nothing
    once |= stmt.head ∈ calllike
    for (i, a) in enumerate(stmt.args)
        isa(a, Expr) || continue
        # Make sure we don't "damage" special syntax that requires literals
        if i == 1 && stmt.head == :foreigncall
            continue
        end
        if i == 2 && stmt.head == :call && stmt.args[1] == :cglobal
            continue
        end
        ret = extract_inner_call!(a, idx, once) # doing this first extracts innermost calls
        ret !== nothing && return ret
        iscalllike = a.head ∈ calllike
        if once && iscalllike
            stmt.args[i] = NewSSAValue(idx)
            return a
        end
    end
    return nothing
end

function replace_ssa!(stmt, ssalookup)
    isa(stmt, Expr) || return nothing
    for (i, a) in enumerate(stmt.args)
        if isa(a, SSAValue)
            stmt.args[i] = SSAValue(ssalookup[a.id])
        elseif isa(a, NewSSAValue)
            stmt.args[i] = SSAValue(a.id)
        else
            replace_ssa!(a, ssalookup)
        end
    end
    return nothing
end

function renumber_ssa!(stmts::Vector{Any}, ssalookup)
    for (i, stmt) in enumerate(stmts)
        if isa(stmt, GotoNode)
            stmts[i] = GotoNode(ssalookup[stmt.label])
        elseif isa(stmt, SSAValue)
            stmts[i] = SSAValue(ssalookup[stmt.id])
        elseif isa(stmt, NewSSAValue)
            stmts[i] = SSAValue(stmt.id)
        elseif isa(stmt, Expr)
            replace_ssa!(stmt, ssalookup)
            if (stmt.head == :gotoifnot || stmt.head == :enter) && isa(stmt.args[end], Int)
                stmt.args[end] = ssalookup[stmt.args[end]]
            end
        end
    end
    return stmts
end

# Pre-frame-construction lookup
function lookup_stmt(stmts, arg)
    if isa(arg, SSAValue)
        arg = stmts[arg.id]
    end
    if isa(arg, QuoteNode)
        arg = arg.value
    end
    return arg
end

function smallest_ref(stmts, arg, idmin)
    if isa(arg, SSAValue)
        idmin = min(idmin, arg.id)
        return smallest_ref(stmts, stmts[arg.id], idmin)
    elseif isa(arg, Expr)
        for a in arg.args
            idmin = smallest_ref(stmts, a, idmin)
        end
    end
    return idmin
end

function lookup_global_refs!(ex::Expr)
    (ex.head == :isdefined || ex.head == :thunk || ex.head == :toplevel) && return nothing
    for (i, a) in enumerate(ex.args)
        if isa(a, GlobalRef)
            r = getfield(a.mod, a.name)
            ex.args[i] = QuoteNode(r)
        elseif isa(a, Expr)
            lookup_global_refs!(a)
        end
    end
    return nothing
end

"""
    optimize!(code::CodeInfo, mod::Module)

Perform minor optimizations on the lowered AST in `code` to reduce execution time
of the interpreter.
Currently it looks up `GlobalRef`s (for which it needs `mod` to know the scope in
which this will run) and ensures that no statement includes nested `:call` expressions
(splitting them out into multiple SSA-form statements if needed).
"""
function optimize!(code::CodeInfo, mod::Module)
    code.inferred && error("optimization of inferred code not implemented")
    # TODO: because of builtins.jl, for CodeInfos like
    #   %1 = Core.apply_type
    #   %2 = (%1)(args...)
    # it would be best to *not* resolve the GlobalRef at %1

    ## Replace GlobalRefs with QuoteNodes
    for (i, stmt) in enumerate(code.code)
        if isa(stmt, GlobalRef)
            code.code[i] = QuoteNode(getfield(stmt.mod, stmt.name))
        elseif isa(stmt, Expr)
            if stmt.head == :call && stmt.args[1] == :cglobal  # cglobal requires literals
                continue
            else
                lookup_global_refs!(stmt)
            end
        end
    end

    ## Un-nest :call expressions (so that there will be only one :call per line)
    # This will allow us to re-use args-buffers rather than having to allocate new ones each time.
    old_code, old_codelocs = code.code, code.codelocs
    code.code = new_code = eltype(old_code)[]
    code.codelocs = new_codelocs = Int32[]
    ssainc = fill(1, length(old_code))
    for (i, stmt) in enumerate(old_code)
        loc = old_codelocs[i]
        inner = extract_inner_call!(stmt, length(new_code)+1)
        while inner !== nothing
            push!(new_code, inner)
            push!(new_codelocs, loc)
            ssainc[i] += 1
            inner = extract_inner_call!(stmt, length(new_code)+1)
        end
        push!(new_code, stmt)
        push!(new_codelocs, loc)
    end
    # Fix all the SSAValues and GotoNodes
    ssalookup = cumsum(ssainc)
    renumber_ssa!(new_code, ssalookup)
    code.ssavaluetypes = length(new_code)

    # Replace :llvmcall and :foreigncall with compiled variants. See
    # https://github.com/JuliaDebug/JuliaInterpreter.jl/issues/13#issuecomment-464880123
    methodtables = Vector{Union{Compiled,TypeMapEntry}}(undef, length(code.code))
    for (idx, stmt) in enumerate(code.code)
        if isexpr(stmt, :call)
            # Check for :llvmcall
            arg1 = stmt.args[1]
            if arg1 == :llvmcall || lookup_stmt(code.code, arg1) == Base.llvmcall
                uuid = uuid4()
                ustr = replace(string(uuid), '-'=>'_')
                methname = Symbol("llvmcall_", ustr)
                nargs = length(stmt.args)-4
                argnames = [Symbol("arg", string(i)) for i = 1:nargs]
                # Run a mini-interpreter to extract the types
                framecode = JuliaFrameCode(CompiledCalls, code; optimize=false)
                frame = prepare_locals(framecode, [])
                idxstart = idx
                for i = 2:4
                    idxstart = smallest_ref(code.code, stmt.args[i], idxstart)
                end
                frame.pc[] = JuliaProgramCounter(idxstart)
                while true
                    pc = step_expr!(Compiled(), frame)
                    convert(Int, pc) == idx && break
                    pc === nothing && error("this should never happen")
                end
                str, RetType, ArgType = @lookup(frame, stmt.args[2]), @lookup(frame, stmt.args[3]), @lookup(frame, stmt.args[4])
                def = quote
                    function $methname($(argnames...))
                        return Base.llvmcall($str, $RetType, $ArgType, $(argnames...))
                    end
                end
                f = Core.eval(CompiledCalls, def)
                stmt.args[1] = QuoteNode(f)
                deleteat!(stmt.args, 2:4)
                methodtables[idx] = Compiled()
            end
        end
    end

    return code, methodtables
end

function prepare_locals(framecode, argvals::Vector{Any})
    if isa(framecode.scope, Method)
        meth, code = framecode.scope::Method, framecode.code
        ssavt = code.ssavaluetypes
        ng = isa(ssavt, Int) ? ssavt : length(ssavt::Vector{Any})
        nargs = length(argvals)
        if !isempty(junk)
            oldframe = pop!(junk)
            locals, ssavalues, sparams = oldframe.locals, oldframe.ssavalues, oldframe.sparams
            exception_frames, last_reference = oldframe.exception_frames, oldframe.last_reference
            callargs = oldframe.callargs
            last_exception, pc = oldframe.last_exception, oldframe.pc
            resize!(locals, length(code.slotflags))
            resize!(ssavalues, ng)
            # for check_isdefined to work properly, we need sparams to start out unassigned
            resize!(sparams, 0)
            empty!(exception_frames)
            empty!(last_reference)
            last_exception[] = nothing
            pc[] = JuliaProgramCounter(1)
        else
            locals = Vector{Union{Nothing,Some{Any}}}(undef, length(code.slotflags))
            ssavalues = Vector{Any}(undef, ng)
            sparams = Vector{Any}(undef, 0)
            exception_frames = Int[]
            last_reference = Dict{Symbol,Int}()
            callargs = Any[]
            last_exception = Ref{Any}(nothing)
            pc = Ref(JuliaProgramCounter(1))
        end
        for i = 1:meth.nargs
            last_reference[framecode.code.slotnames[i]] = i
            if meth.isva && i == meth.nargs
                locals[i] = nargs < i ? Some{Any}(()) : (let i=i; Some{Any}(ntuple(k->argvals[i+k-1], nargs-i+1)); end)
                break
            end
            locals[i] = nargs >= i ? Some{Any}(argvals[i]) : Some{Any}(())
        end
        # add local variables initially undefined
        for i = (meth.nargs+1):length(code.slotnames)
            locals[i] = nothing
        end
    else
        code = framecode.code
        locals = Vector{Union{Nothing,Some{Any}}}(undef, length(code.slotflags))
        fill!(locals, nothing)
        ssavalues = Vector{Any}(undef, length(code.code))
        sparams = Any[]
        exception_frames = Int[]
        last_reference = Dict{Symbol,Int}()
        callargs = Any[]
        last_exception = Ref{Any}(nothing)
        pc = Ref(JuliaProgramCounter(1))
    end
    JuliaStackFrame(framecode, locals, ssavalues, sparams, exception_frames, last_exception,
                    pc, last_reference, callargs)
end

"""
    frame = build_frame(framecode::JuliaFrameCode, frameargs, lenv)

Construct a new `JuliaStackFrame` for `framecode`, given lowered-code arguments `frameargs` and
static parameters `lenv`. See [`JuliaInterpreter.prepare_call`](@ref) for information about how to prepare the inputs.
"""
function build_frame(framecode, args, lenv)
    frame = prepare_locals(framecode, args)
    resize!(frame.sparams, length(lenv))
    # Add static parameters to environment
    for i = 1:length(lenv)
        T = lenv[i]
        isa(T, TypeVar) && continue  # only fill concrete types
        frame.sparams[i] = T
    end
    return frame
end

"""
    frame = enter_call_expr(expr; enter_generated=false)

Build a `JuliaStackFrame` ready to execute the expression `expr`. Set `enter_generated=true`
if you want to execute the generator of a `@generated` function, rather than the code that
would be created by the generator.

# Example

```jldoctest; setup=(using JuliaInterpreter; empty!(JuliaInterpreter.junk))
julia> mymethod(x) = x+1
mymethod (generic function with 1 method)

julia> JuliaInterpreter.enter_call_expr(:(\$mymethod(1)))
JuliaStackFrame(JuliaInterpreter.JuliaFrameCode(mymethod(x) in Main at none:1, CodeInfo(
1 ─ %1 = ($(QuoteNode(+)))(x, 1)
└──      return %1
), Union{Compiled, TypeMapEntry}[#undef, #undef], JuliaInterpreter.BreakpointState[#undef, #undef], BitSet([1]), false, false, true), Union{Nothing, Some{Any}}[Some(mymethod), Some(1)], Any[#undef, #undef], Any[], Int64[], Base.RefValue{Any}(nothing), Base.RefValue{JuliaInterpreter.JuliaProgramCounter}(JuliaProgramCounter(1)), Dict(Symbol("#self#")=>1,:x=>2), Any[])

julia> mymethod(x::Vector{T}) where T = 1
mymethod (generic function with 2 methods)

julia> a = [1.0, 2.0]
2-element Array{Float64,1}:
 1.0
 2.0

julia> JuliaInterpreter.enter_call_expr(:(\$mymethod(\$a)))
JuliaStackFrame(JuliaInterpreter.JuliaFrameCode(mymethod(x::Array{T,1}) where T in Main at none:1, CodeInfo(
1 ─     return 1
), Union{Compiled, TypeMapEntry}[#undef], JuliaInterpreter.BreakpointState[#undef], BitSet([]), false, false, true), Union{Nothing, Some{Any}}[Some(mymethod), Some([1.0, 2.0])], Any[#undef], Any[Float64], Int64[], Base.RefValue{Any}(nothing), Base.RefValue{JuliaInterpreter.JuliaProgramCounter}(JuliaProgramCounter(1)), Dict(Symbol("#self#")=>1,:x=>2), Any[])
```

See [`enter_call`](@ref) for a similar approach not based on expressions.
"""
function enter_call_expr(expr; enter_generated = false)
    r = determine_method_for_expr(expr; enter_generated = enter_generated)
    if isa(r, Tuple)
        return build_frame(r[1:end-1]...)
    end
    nothing
end

"""
    frame = enter_call(f, args...; kwargs...)

Build a `JuliaStackFrame` ready to execute `f` with the specified positional and keyword arguments.

# Example

```jldoctest; setup=(using JuliaInterpreter; empty!(JuliaInterpreter.junk))
julia> mymethod(x) = x+1
mymethod (generic function with 1 method)

julia> JuliaInterpreter.enter_call(mymethod, 1)
JuliaStackFrame(JuliaInterpreter.JuliaFrameCode(mymethod(x) in Main at none:1, CodeInfo(
1 ─ %1 = ($(QuoteNode(+)))(x, 1)
└──      return %1
), Union{Compiled, TypeMapEntry}[#undef, #undef], JuliaInterpreter.BreakpointState[#undef, #undef], BitSet([1]), false, false, true), Union{Nothing, Some{Any}}[Some(mymethod), Some(1)], Any[#undef, #undef], Any[], Int64[], Base.RefValue{Any}(nothing), Base.RefValue{JuliaInterpreter.JuliaProgramCounter}(JuliaProgramCounter(1)), Dict(Symbol("#self#")=>1,:x=>2), Any[])

julia> mymethod(x::Vector{T}) where T = 1
mymethod (generic function with 2 methods)

julia> JuliaInterpreter.enter_call(mymethod, [1.0, 2.0])
JuliaStackFrame(JuliaInterpreter.JuliaFrameCode(mymethod(x::Array{T,1}) where T in Main at none:1, CodeInfo(
1 ─     return 1
), Union{Compiled, TypeMapEntry}[#undef], JuliaInterpreter.BreakpointState[#undef], BitSet([]), false, false, true), Union{Nothing, Some{Any}}[Some(mymethod), Some([1.0, 2.0])], Any[#undef], Any[Float64], Int64[], Base.RefValue{Any}(nothing), Base.RefValue{JuliaInterpreter.JuliaProgramCounter}(JuliaProgramCounter(1)), Dict(Symbol("#self#")=>1,:x=>2), Any[])
```

For a `@generated` function you can use `enter_call((f, true), args...; kwargs...)`
to execute the generator of a `@generated` function, rather than the code that
would be created by the generator.

See [`enter_call_expr`](@ref) for a similar approach based on expressions.
"""
function enter_call(@nospecialize(finfo), @nospecialize(args...); kwargs...)
    if isa(finfo, Tuple)
        f = finfo[1]
        enter_generated = finfo[2]::Bool
    else
        f = finfo
        enter_generated = false
    end
    f, allargs = prepare_args(f, Any[f, args...], kwargs)
    # Can happen for thunks created by generated functions
    if isa(f, Core.Builtin) || isa(f, Core.IntrinsicFunction)
        error(f, " is a builtin or intrinsic")
    end
    r = prepare_call(f, allargs; enter_generated=enter_generated)
    if isa(r, Tuple)
        return build_frame(r[1:end-1]...)
    end
    return nothing
end

lower(mod, arg) = false ? expand(arg) : Meta.lower(mod, arg)

# This is a version of gen_call_with_extracted_types, except that is passes back the call expression
# for further processing.
function extract_args(__module__, ex0)
    if isa(ex0, Expr)
        kws = collect(filter(x->isexpr(x,:kw),ex0.args))
        if !isempty(kws)
            names = []
            values = Tuple(map(x-> begin
                push!(names,x.args[1])
                x.args[2]
            end,kws))
            names = Tuple(names)
            return Expr(:tuple,:(Core.kwfunc($(ex0.args[1]))),
                Expr(:call, NamedTuple{names,typeof(values)}, values),
                map(x->isexpr(x, :parameters) ? QuoteNode(x) : x,
                filter(x->!isexpr(x, :kw),ex0.args))...)
        elseif ex0.head == :.
            return Expr(:tuple, :getproperty, ex0.args...)
        elseif ex0.head == :(<:)
            return Expr(:tuple, :(<:), ex0.args...)
        else
            return Expr(:tuple,
                map(x->isexpr(x,:parameters) ? QuoteNode(x) : x, ex0.args)...)
        end
    end
    if isa(ex0, Expr) && ex0.head == :macrocall # Make @edit @time 1+2 edit the macro by using the types of the *expressions*
        return error("Macros are not supported in @enter")
    end
    ex = lower(__module__, ex0)
    exret = Expr(:none)
    if !isa(ex, Expr)
        return error("expression is not a function call or symbol")
    elseif ex.head == :call
        return Expr(:tuple,
                map(x->isexpr(x,:parameters) ? QuoteNode(x) : x, ex.args)...)
    elseif ex.head == :body
        a1 = ex.args[1]
        if isa(a1, Expr) && a1.head == :call
            a11 = a1.args[1]
            if a11 == :setindex!
                return Expr(:tuple,
                map(x->isexpr(x,:parameters) ? QuoteNode(x) : x, arg.args)...)
            end
        end
    end
    return error("expression is not a function call, "
               * "or is too complex for @enter to analyze; "
               * "break it down to simpler parts if possible")
end

"""
    @interpret f(args; kwargs...)

Evaluate `f` on the specified arguments using the interpreter.

# Example

```jldoctest; setup=(using JuliaInterpreter; empty!(JuliaInterpreter.junk))
julia> a = [1, 7]
2-element Array{Int64,1}:
 1
 7

julia> sum(a)
8

julia> @interpret sum(a)
8
```
"""
macro interpret(arg)
    args = try
        extract_args(__module__, arg)
    catch e
        return :(throw($e))
    end
    quote
        theargs = $(esc(args))
        stack = JuliaStackFrame[]
        frame = JuliaInterpreter.enter_call_expr(Expr(:call, theargs...))
        if frame === nothing
            return eval(Expr(:call, map(QuoteNode, theargs)...))
        end
        if shouldbreak(frame, 1)
            push!(stack, frame)
            return stack, BreakpointRef(frame.code, 1)
        end
        ret = finish_and_return!(stack, frame)
        isa(ret, BreakpointRef) ? (stack, ret) : ret
    end
end

function set_compiled_methods()
    # Work around #28 by preventing interpretation of all Base methods that have a ccall to memcpy
    push!(compiled_methods, which(vcat, (Vector,)))
    push!(compiled_methods, first(methods(Base._getindex_ra)))
    push!(compiled_methods, first(methods(Base._setindex_ra!)))
    push!(compiled_methods, which(Base.decompose, (BigFloat,)))
    push!(compiled_methods, which(DSFMT.dsfmt_jump, (DSFMT.DSFMT_state, DSFMT.GF2X)))
    if Sys.iswindows()
        push!(compiled_methods, which(InteractiveUtils.clipboard, (AbstractString,)))
    end
    # issue #76
    push!(compiled_methods, which(unsafe_store!, (Ptr{Any}, Any, Int)))
    push!(compiled_methods, which(unsafe_store!, (Ptr, Any, Int)))
    # issue #92
    push!(compiled_methods, which(objectid, Tuple{Any}))
end

function __init__()
    set_compiled_methods()
end

include("breakpoints.jl")
using .Breakpoints
using .Breakpoints: shouldbreak, BreakpointRef

include("precompile.jl")
_precompile_()

end # module


text:
include("testplotting.jl")
include("Para_Didi_test.jl")
println("Tests passed!")

ext: jl


text:
module MyMatching

export my_deferred_acceptance

#many to many

function my_deferred_acceptance(prop_prefs::Vector{Vector{Int}},
                             resp_prefs::Vector{Vector{Int}},
                             p_caps::Vector{Int},
                             r_caps::Vector{Int})
    m = length(prop_prefs)
    n = length(resp_prefs)
    prop_matched = zeros(Int64,sum(p_caps))
    resp_matched = zeros(Int64,sum(r_caps))
    count = ones(Int64,m) #提案側の提案回数
    accept = zeros(Int64,n) #受入側の受入人数
    p_count = zeros(Int64,n) #受入側が提案された回数（選好順位の考察に使用）

#iがマッチングしたリストを表すための変数
    p_indptr = Array{Int}(m+1)
    r_indptr = Array{Int}(n+1)
    p_indptr[1] = 1
    for i in 1:m
        p_indptr[i+1] = p_indptr[i] + p_caps[i]
    end
    r_indptr[1] = 1
    for i in 1:n
        r_indptr[i+1] = r_indptr[i] + r_caps[i]
    end

    for j in 1:sum(p_caps)*sum(r_caps)
        for i in 1:m
            p_list = prop_matched[p_indptr[i]:p_indptr[i+1]-1]
            if count[i] <= length(prop_prefs[i]) && findfirst(p_list, 0) != 0
                k = prop_prefs[i][count[i]] #このラウンドでiが提案する相手k
                if findfirst(resp_prefs[k], i) != 0
                    if accept[k] < r_caps[k]
                        resp_matched[r_indptr[k+1]-1 - accept[k]] = i
                        prop_matched[p_indptr[i] + findfirst(p_list, 0) - 1] = k
                        accept[k]+=1
                    else
                        list = resp_matched[r_indptr[k]:r_indptr[k+1]-1] #受入側kの保留相手のリスト
                        ranking = zeros(Int64,r_caps[k]) #リスト内での保留相手の順序とkのそれぞれの選好順位を対応させた
                        for l in 1:r_caps[k]
                            ranking[l] = findfirst(resp_prefs[k], list[l])
                        end
                        if 0 < findfirst(resp_prefs[k], i) < maximum(ranking)
                            d_list = prop_matched[p_indptr[list[indmax(ranking)]]:p_indptr[list[indmax(ranking)]+1]-1] #保留から外される提案側のマッチングリスト
                            resp_matched[r_indptr[k] + indmax(ranking) - 1] = i
                            prop_matched[p_indptr[list[indmax(ranking)]] + findfirst(d_list, k) - 1] = 0
                            prop_matched[p_indptr[i] + findfirst(p_list, 0) - 1] = k
                        end
                    end
                end
                count[i]+=1
                p_count[k]+=1
            end
        end
    end
    return prop_matched, resp_matched, p_indptr, r_indptr, p_count
end

#many to one
function my_deferred_acceptance(prop_prefs::Vector{Vector{Int}},
                                resp_prefs::Vector{Vector{Int}},
                                r_caps::Vector{Int})
    p_caps = ones(Int, length(prop_prefs))
    prop_matched, resp_matched, p_indptr, r_indptr =
        my_deferred_acceptance(prop_prefs, resp_prefs, p_caps, r_caps)
    return prop_matched, resp_matched, r_indptr
end

#one to one
function my_deferred_acceptance(prop_prefs::Vector{Vector{Int}},
                                resp_prefs::Vector{Vector{Int}})
    p_caps = ones(Int, length(prop_prefs))
    r_caps = ones(Int, length(resp_prefs))
    prop_matched, resp_matched, p_indptr, r_indptr, p_count =
        my_deferred_acceptance(prop_prefs, resp_prefs, p_caps, r_caps)
    return prop_matched, resp_matched, p_count
end

end


text:
"""
A direct translation of the Matlab function mlf by Igor Podlubny and Martin
Kacenak.  https://www.mathworks.com/matlabcentral/fileexchange/8738-mittag-leffler-function
"""
module MLF

using ..MittagLefflerFunctions: rΓ

export mlf

"""
    mlf(α, β, z, p)

Evaluate the Mittag-Leffler function E_αβ(z) with accuracy 10^(-p).
"""
function mlf(α::T, β::T, z::Complex{T}, p::Integer) where T <: AbstractFloat
    ten = parse(T, "10.0")
    r0 = radius(α, β, z, p)
    if α == one(T) && β == one(T)
        E = exp(z)
    else
        if (   ( α < 1 && abs(z) ≤ 1 )
            || ( ( 1 ≤ α < 2 ) && abs(z) ≤ floor(20/(2.1-α)^(5.5-2α)))
            || (  α ≥ 2 && abs(z) ≤ 50 ) )
            E = mlf1(α, β, z, p)
        else
            if  α<=1 && abs(z) ≤ floor(5*α+ten)
                if (   ( abs(angle(z)) > π*α )
                    && ( abs(abs(angle(z)) - (π*α)) > ten^(-p) ) )
                    if β ≤ 1
                        E = mlf2(α, β, z, p, r0)
                    else
                        E = mlf3(α, β, z, p, r0)
                    end
                elseif (   abs(angle(z)) < π*α
                        && abs(abs(angle(z))-(π*α)) > ten^(-p) )
                    if β <= 1
                        E = mlf4(α, β, z, p, r0)
                    else
                        E = mlf5(α, β, z, p, r0)
                    end
                else
                    E = mlf6(α, β, z, p, r0)
                end
            else
                if α ≤ 1
                    if abs(angle(z)) < (π*α/2+min(π,π*α))/2
                        E = mlf7(α, β, z, p)
                    else
                        E = mlf8(α, β, z, p)
                    end
                else
                    if α ≥ 2
                        E = mlf9(α, β, z, p)
                    else
                        E = mlf10(α, β, z, p)
                    end
                end
            end
        end
    end
    return E
end

function mlf(α::T, β::T, x::T, p::Integer) where T <: AbstractFloat
    z = Complex(x, zero(T))
    return real(mlf(α, β, z, p))
end

function radius(α::T, β::T, z::Complex{T}, p::Integer) where T <: AbstractFloat
    ten = parse(T, "10.0")
    if β < 0
        rc = ( -2*log( ten^(-p)*π / (6*(abs(β)+2)*(2*abs(β))^(abs(β))) ) )^α
    else
        rc = ( -2*log( ten^(-p)*π / 6 ) )^α
    end
    r0 = max(one(T), 2*abs(z), rc)
    return r0
end

function mlf1(α::T, β::T, z::Complex{T}, p::Integer) where T <: AbstractFloat
    oldsum=0
    k = 0
    powz = Complex(one(T),zero(T))
    while β + k*α ≤ 0
        powz *= z # = z^k
        k += 1
    end
    newsum = powz * rΓ(β+k*α)
    while newsum != oldsum
        oldsum = newsum
        k += 1
        powz *= z
        term = powz * rΓ(β+k*α)
        newsum += term
        k += 1
        powz *= z
        term = powz * rΓ(β+k*α)
        newsum += term
    end
    E = newsum
    return E
end

function mlf2(α::T, β::T, z::Complex{T}, p::Integer,
              r0::T) where T <: AbstractFloat
    E = Romberg(zero(T), r0, p) do r
        return K(r, α, β, z)
    end
    return E
end

function mlf3(α::T, β::T, z::Complex{T}, p::Integer,
              r0::T) where T <: AbstractFloat
    ϵ = one(T)
    E1 = Romberg(ϵ, r0, p) do r
        return K(r, α, β, z)
    end
    E2 = Romberg(-π*α, π*α, p) do r
        return P(r, α, β, z, ϵ)
    end
    E = E1 + E2
    return E
end

function mlf4(α::T, β::T, z::Complex{T}, p::Integer,
              r0::T) where T <: AbstractFloat
    E1 = Romberg(zero(T), r0, p) do r
        return K(r, α, β, z)
    end
    E = E1 + z^((1-β)/α) * exp(z^(1/α)) / α
    return E
end

function mlf5(α::T, β::T, z::Complex{T}, p::Integer,
              r0::T) where T <: AbstractFloat
    ϵ = abs(z)/2
    E1 = Romberg(ϵ, r0, p) do r
        return K(r, α, β, z)
    end
    E2 = Romberg(-π*α, π*α, p) do r
        return P(r, α, β, z, ϵ)
    end
    E = E1 + E2 + (z^((1-β)/α)) * (exp(z^(1/α))/α)
    return E
end

function mlf6(α::T, β::T, z::Complex{T}, p::Integer,
              r0::T) where T <: AbstractFloat
    ϵ = abs(z) + one(T)/2
    E1 = Romberg(ϵ, r0, p) do r
        return K(r, α, β, z)
    end
    E2 = Romberg(-π*α, π*α, p) do r
        return P(r, α, β, z, ϵ)
    end
    E = E1 + E2
    return E
end

function mlf7(α::T, β::T, z::Complex{T}, p::Integer) where T <: AbstractFloat
    newsum = z^((1-β)/α) * exp(z^(1/α))/α
    powz = Complex(one(T), zero(T))
    for k=1:floor(Int64, p/log10(abs(z)))
        powz /= z # = z^(-k)
        newsum -= powz * rΓ(β-k*α)
    end
    E = newsum
    return E
end

function mlf8(α::T, β::T, z::Complex{T}, p::Integer) where T <: AbstractFloat
    newsum = zero(Complex{T})
    powz = Complex(one(T), zero(T))
    for k=1:floor(Int64, p/log10(abs(z)))
        powz /= z # = z^(-k)
        newsum -= powz * rΓ(β-k*α)
    end
    E = newsum
    return E
end

function mlf9(α::T, β::T, z::Complex{T}, p::Integer) where T <: AbstractFloat
    m = floor(Int64, α/2)
    s = zero(Complex{T})
    for h = 0:m
        zn = z^(one(T)/(m+1)) * exp_i(2π*h/(m+1))
        s += mlf(α/(m+1), β, zn, p)
    end
    E = (one(T)/(m+1)) * s
    return E
end

function mlf10(α::T, β::T, z::Complex{T}, p::Integer) where T <: AbstractFloat
    E = ( mlf(α/2, β, sqrt(z), p)
        + mlf(α/2, β, -sqrt(z), p) ) / 2
    return E
end

function exp_i(t::T) where T <: AbstractFloat
    return exp(Complex(zero(T),t))
end

function K(r::T, α::T, β::T, z::Complex{T}) where T <: AbstractFloat
    return r^((1-β)/α) * exp(-r^(1/α)) * ( (r*sinpi(1-β)-z*sinpi(1-β+α))
                                          /(π*α*(r^2-2*r*z*cospi(α)+z^2)) )
end

function P(r::T, α::T, β::T, z::Complex{T}, ϵ::T) where T <: AbstractFloat
    w = ϵ^(1/α) * sin(r/α) + r * ( 1 + (1-β)/α )
    numer =  exp((ϵ^(1/α)) * cos(r/α)) * exp_i(w)
    return ( ϵ^(1+(1-β)/α) / (2π*α) ) * numer / ( ϵ*exp_i(r) - z )
end

function Romberg(f::Function, a::T, b::T, order=6) where T <: AbstractFloat
    R = zeros(Complex{T}, 2, order)
    h = b - a
    R[1,1] = ( h/2 ) * ( f(a) + f(b) )
    ipower = 1
    for i = 2:order
        s = zero(Complex{T})
        for j = 1:ipower
            s += f( a + (2j-1)*h/2 )
        end
        R[2,1] = ( R[1,1] + h * s ) / 2
        pow4 = one(T)
        for k = 1:i-1
            pow4 *= 4
            R[2,k+1] = ( pow4 * R[2,k] - R[1,k] ) / ( pow4 - 1 )
        end
        for j = 0:i-1
            R[1,j+1] = R[2,j+1]
        end
        ipower *= 2
        h /= 2
    end
    return R[1,order]
end

end # module


text:
const Space = SVector{2, Float64}

struct SpaceTime
    xy::Space
    t::Float64
end
SpaceTime(x, y, t) = SpaceTime(Space(x, y), t)

# const SpaceTime = typeof(StructArray((xy = Space.(rand(2), rand(2)), t = rand(2))))

struct POI
    xyt::StructVector{SpaceTime}
    video::SystemPath
end
POI(x, y, t, v) = POI(StructVector(SpaceTime.(x, y, t)), v)
space(x::POI) = x.xyt.xy
time(x::POI) = x.xyt.t
space(x::POI, i) = x.xyt[i].xy
time(x::POI, i) = x.xyt[i].t

correctedges(x, nframes) = x == nframes ? nframes - 1 :
                           x == 1 ? 2 :
                           x

function resfile2coords(resfile, videofile, poi_names)
    matopen(string(resfile)) do io
        xdata = read(io, "xdata")
        fr = read(io, "status")["FrameRate"]
        nframes = size(xdata, 1)
        rows = rowvals(xdata)
        xvals = nonzeros(xdata)
        yvals = nonzeros(read(io, "ydata"))
        pois = Dict{Symbol, POI}()
        for (j, name) in enumerate(poi_names)
            i = nzrange(xdata, j)
            if !isempty(i) && !all(isnumeric, String(name))
                r = map(x -> correctedges(x, nframes), rows[i])
                pois[name] = POI(xvals[i], yvals[i], r/fr, videofile)
            end
        end
        return pois
    end
end

cleanpoi(x::POI) = length(x.xyt) == 1 ? space(x, 1) : x.xyt


text:
push!(LOAD_PATH,"../src/")

using Documenter, CvxCompress

makedocs(
    sitename="CvxCompress.jl",
    modules=[CvxCompress],
    pages = [
        "Home" => "index.md",
        "User Guide" => "manual.md",
        "reference.md",
        "Performance" => "perf/perf.md"
        ]
)

deploydocs(
    repo = "github.com/ChevronETC/CvxCompress.jl"
)


text:
"""
    grdinfo(cmd0::String="", arg1=[]; kwargs...)

Reads a 2-D grid file and reports metadata and various statistics for the (x,y,z) data in the grid file

Full option list at [`grdinfo`](http://gmt.soest.hawaii.edu/doc/latest/grdinfo.html)

Parameters
----------

- **C** : **numeric** : -- Str or Number --

    Formats the report using tab-separated fields on a single line.
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/grdinfo.html#c)
- **D** : **tiles** : -- Number or Str --

    Divide a single grid’s domain (or the -R domain, if no grid given) into tiles of size
    dx times dy (set via -I).
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/grdinfo.html#d)
- **F** : -- Bool or [] --

    Report grid domain and x/y-increments in world mapping format.
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/grdinfo.html#f)
- **I** : **nearest_multiple** : -- Number or Str --

    Report the min/max of the region to the nearest multiple of dx and dy, and output
    this in the form -Rw/e/s/n
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/grdinfo.html#i)
- **L** : **force_scan** : -- Number or Str --

    Report stats after actually scanning the data.
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/grdinfo.html#l)
- **M** : **zmin_max** : -- Bool or [] --

    Find and report the location of min/max z-values.
    [`-M`](http://gmt.soest.hawaii.edu/doc/latest/grdinfo.html#m)
- $(GMT.opt_R)
- **T** : **nan_t** : -- Number or Str --
    Determine min and max z-value.
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/grdinfo.html#t)
- $(GMT.opt_V)
- $(GMT.opt_f)
"""
function grdinfo(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && !isa(arg1, GMTgrid) && return monolitic("grdinfo", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)
	cmd, = parse_R("", d)
	cmd = parse_V_params(cmd, d)
	cmd, = parse_f(cmd, d)

	cmd = add_opt(cmd, 'C', d, [:C :numeric])
	cmd = add_opt(cmd, 'D', d, [:D :tiles])
	cmd = add_opt(cmd, 'F', d, [:F])
	cmd = add_opt(cmd, 'I', d, [:I :nearest_multiple])
	cmd = add_opt(cmd, 'L', d, [:L :force_scan])
	cmd = add_opt(cmd, 'M', d, [:M :location])
	cmd = add_opt(cmd, 'T', d, [:T :zmin_max])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "grdinfo", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grdinfo(arg1, cmd0::String=""; kw...) = grdinfo(cmd0, arg1; kw...)


text:
module AgnosticBayesEnsemble
  using Optim
  export TDistPosteriorEstimation, GMatrix,
         dirichletPosteriorEstimation,
         dirichletPosteriorEstimationV2, dirichletPosteriorEstimation!,
         metaParamSearchValidationDirichlet,
         bootstrapPosteriorEstimation, bootstrapPosteriorEstimation!,
         bootstrapPosteriorCorEstimation,
         TDistPosteriorEstimation, TDistPosteriorEstimationReference,
         predictEnsemble,
         directOptimNaiveMSE, directOptimHinge,
         posteriorLinearBasis

  include( "argmaxProb.jl" )
  include( "argminProb.jl" )
  include( "bootstrapPosteriorCorEstimation.jl" )
  include( "bootstrapPosteriorEstimation.jl" )
  include( "directSolution.jl" )
  include( "dirichletPosteriorEstimation.jl" )
  include( "predictEnsemble.jl" )

end


text:
using DrWatson
@quickactivate :ConditionalNFs4Imaging

using JLD
using PyPlot
using ArgParse
using ProgressMeter

function create_OOD_shots_and_RTM(args::Dict)

    sim_name = "create_OOD_shots_and_RTM"
    data_args = Dict(
        "nsrc" => args["nsrc"],
        "nrec" => args["nrec"],
        "sigma" => args["sigma"],
        "sim_name" => sim_name,
    )

    dm, lin_data, rtm = load_shot_data(data_args)

    if (dm == lin_data == rtm == nothing)

        sigma = args["sigma"]
        idx = args["idx"]
        nsrc = args["nsrc"]
        nrec = args["nrec"]

        # Define raw data directory
        mkpath(datadir("training-data"))
        data_path = datadir("training-data", "seismic_samples_256_by_256_num_10k.jld")
        label_path =
            datadir("training-data", "seismic_samples_256_by_256_num_10k_labels.jld")

        # Download the dataset into the data directory if it does not exist
        if isfile(data_path) == false
            run(`wget https://www.dropbox.com/s/vapl62yhh8fxgwy/'
                'seismic_samples_256_by_256_num_10k.jld -q -O $data_path`)
        end
        if isfile(label_path) == false
            run(`wget https://www.dropbox.com/s/blxkh6tdszudlcq/'
                'seismic_samples_256_by_256_num_10k_labels.jld -q -O $label_path`)
        end

        # Load OOD seismic images
        X_OOD = JLD.jldopen(data_path, "r") do file
            X_OOD = read(file, "X")

            labels = JLD.jldopen(label_path, "r")["labels"][:]
            idx = findall(x -> x == -1, labels)[idx]
            return X_OOD[:, :, :, idx:idx]
        end

        dm = X_OOD[:, :, 1, 1]'
        dm[:, 1:10] .= 0.0
        dm = convert(Array{Float32,1}, vec(dm))

        # Create forward modeling operator
        J, noise_dist =
            create_operator(; nrec = nrec, nsrc = nsrc, sigma = sigma, sim_src = false)
        lin_data = zeros(Float32, nsrc, J.source.geometry.nt[1] * nrec)
        rtm = zeros(Float32, size(dm))

        p = Progress(nsrc)
        for j = 1:nsrc
            Base.flush(Base.stdout)

            lin_data[j, :] = J[j] * dm + rand(noise_dist)
            grad = adjoint(J[j]) * lin_data[j, :]
            rtm = (rtm * (j - 1) + grad) / j

            ProgressMeter.next!(p; showvalues = [(:source, j)])
        end

        dm = reshape(reshape(dm, size(X_OOD)[1:2])', size(X_OOD)[1], size(X_OOD)[2], 1, 1)
        rtm = reshape(reshape(rtm, size(X_OOD)[1:2])', size(X_OOD)[1], size(X_OOD)[2], 1, 1)

        save_dict = @strdict sigma idx nsrc nrec dm lin_data rtm sim_name
        @tagsave(
            datadir(sim_name, savename(save_dict, "jld2"; digits = 6)),
            save_dict;
            safe = true
        )

    end
    return dm, lin_data, rtm
end


text:
function simulate_fusion(seq, chr, start, profile_fusion, assembly)
    seqstr = seq.seq
    len = length(seq)
    changed = false
    l=0
    r=0
    for fusion in profile_fusion
        on_fusion = false
        if is_in_fusion(fusion["left"], chr, start, len)
            on_fusion = true
            if fusion["left"]["strand"]=="+"
                l = fusion["left"]["pos"] - start
            else
                l = start - fusion["left"]["pos"]
            end
            l = min(len-1, max(1, l))
            r = len - l
        elseif is_in_fusion(fusion["right"], chr, start, len)
            on_fusion = true
            if fusion["right"]["strand"]=="+"
                r = (start+len) - fusion["right"]["pos"]
            else
                r = fusion["right"]["pos"] - (start-len)
            end
            r = min(len-1, max(1, r))
            l = len - r
        end
        if on_fusion && rand() < fusion["rate"]
            return make_fusion_seq(l, r, fusion["left"], fusion["right"], assembly)
        end
    end
    return seq
end

function make_fusion_seq(l, r, fusion_left, fusion_right, assembly)
    leftseq = nothing
    rightseq = nothing
    if fusion_left["strand"]=="+"
        leftseq = assembly[fusion_left["chrom"]].sequence[fusion_left["pos"]-l-1:fusion_left["pos"]]
    else
        leftseq = assembly[fusion_left["chrom"]].sequence[fusion_left["pos"]:fusion_left["pos"]+l-1]
        leftseq = ~leftseq
    end
    if fusion_right["strand"]=="+"
        rightseq = assembly[fusion_right["chrom"]].sequence[fusion_right["pos"]:fusion_right["pos"]+r-1]
    else
        rightseq = assembly[fusion_right["chrom"]].sequence[fusion_right["pos"]-r-1:fusion_right["pos"]]
        rightseq = ~rightseq
    end
    # connect left/right seq to a fusion seq
    fusion_seq = dna(uppercase(leftseq.seq * rightseq.seq))
    # println(">fusion_left:", l, "_right:", r)
    # println(fusion_seq.seq)
    return fusion_seq
end

function is_in_fusion(fusion, chr, start, len)
    if fusion["strand"] == "+"
        return fusion["chrom"] == chr && start < fusion["pos"] && fusion["pos"] < start + len
    else
        return fusion["chrom"] == chr && start-len < fusion["pos"] && fusion["pos"] < start
    end
end

function reverse(fusion)
    strand = "+"
    if fusion["strand"] == "+"
        strand = "-"
    end
    return Dict("chrom"=>fusion["chrom"], "pos"=>fusion["pos"], "strand"=>strand)
end


text:
############################################################################################
### Types
############################################################################################
"""
Empty placeholder for models that don't need an imply part.

# Constructor

    ImplyEmpty(;specification, kwargs...)

# Arguments
- `specification`: either a `RAMMatrices` or `ParameterTable` object

# Extended help

## Interfaces
- `identifier(::RAMSymbolic) `-> Dict containing the parameter labels and their position
- `n_par(::RAMSymbolic)` -> Number of parameters

## Implementation
Subtype of `SemImply`.
"""
struct ImplyEmpty{V, V2} <: SemImply
    identifier::V2
    n_par::V
end

############################################################################################
### Constructors
############################################################################################

function ImplyEmpty(;
        specification,
        kwargs...)

        ram_matrices = RAMMatrices(specification)
        identifier = StructuralEquationModels.identifier(ram_matrices)

        n_par = length(ram_matrices.parameters)

        return ImplyEmpty(identifier, n_par)
end

############################################################################################
### methods
############################################################################################

objective!(imply::ImplyEmpty, par, model) = nothing
gradient!(imply::ImplyEmpty, par, model) = nothing
hessian!(imply::ImplyEmpty, par, model) = nothing

############################################################################################
### Recommended methods
############################################################################################

identifier(imply::ImplyEmpty) = imply.identifier
n_par(imply::ImplyEmpty) = imply.n_par

update_observed(imply::ImplyEmpty, observed::SemObserved; kwargs...) = imply


text:
export Transition;

struct Transition
  move::Move
  label::String
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
using Jchemo
using Test


text:
# Core functionality for the `BlockDiagonal` type

"""
    BlockDiagonal{T, V<:AbstractMatrix{T}} <: AbstractMatrix{T}

A matrix with matrices on the diagonal, and zeros off the diagonal.
"""
struct BlockDiagonal{T, V<:AbstractMatrix{T}} <: AbstractMatrix{T}
    blocks::Vector{V}

    function BlockDiagonal{T, V}(blocks::Vector{V}) where {T, V<:AbstractMatrix{T}}
        return new{T, V}(blocks)
    end
end

function BlockDiagonal(blocks::Vector{V}) where {T, V<:AbstractMatrix{T}}
    return BlockDiagonal{T, V}(blocks)
end

function ChainRulesCore.rrule(::Type{<:BlockDiagonal}, blocks::Vector{V}) where {V}
    BlockDiagonal_pullback(Δ::Composite) = (NO_FIELDS, Δ.blocks)
    return BlockDiagonal(blocks), BlockDiagonal_pullback
end

BlockDiagonal(B::BlockDiagonal) = B

is_square(A::AbstractMatrix) = size(A, 1) == size(A, 2)

"""
    blocks(B::BlockDiagonal{T, V}) -> Vector{V}

Return the on-diagonal blocks of B.
"""
blocks(B::BlockDiagonal) = B.blocks

# BlockArrays-like functions
"""
    blocksizes(B::BlockDiagonal) -> Vector{Tuple}

Return the size of each on-diagonal block in order.

# Example
```jldoctest; setup = :(using BlockDiagonals)
julia> B = BlockDiagonal([rand(2, 2), rand(3, 3)]);

julia> blocksizes(B)
2-element Array{Tuple{Int64,Int64},1}:
 (2, 2)
 (3, 3)
```
See also [`blocksize`](@ref) for accessing the size of a single block efficiently.
"""
blocksizes(B::BlockDiagonal) = map(size, blocks(B))

"""
    blocksize(B::BlockDiagonal, p::Integer, q::Integer=p) -> Tuple

Return the size of the p^th on-diagonal block. Optionally specify `q` to return the
size of block `p, q`.

# Example
```jldoctest; setup = :(using BlockDiagonals)
julia> X = rand(2, 2); Y = rand(3, 3);

julia> B = BlockDiagonal([X, Y]);

julia> blocksize(B, 1)
(2, 2)

julia> blocksize(B, 1, 2)
(2, 3)
```
See also [`blocksizes`](@ref) for accessing the size of all on-diagonal blocks easily.
"""
blocksize(B::BlockDiagonal, p::Integer) = size(blocks(B)[p])
function blocksize(B::BlockDiagonal, p::Integer, q::Integer)
    return size(blocks(B)[p], 1), size(blocks(B)[q], 2)
end

"""
    nblocks(B::BlockDiagonal[, dim])

Return the number of on-diagonal blocks.

The total number of blocks in the matrix is `nblocks(B)^2`.
"""
nblocks(B::BlockDiagonal) = length(blocks(B))

getblock(B::BlockDiagonal, p::Integer) = blocks(B)[p]
function getblock(B::BlockDiagonal{T}, p::Integer, q::Integer) where T
    return p == q ? blocks(B)[p] : Zeros{T}(blocksize(B, p, q))
end

function setblock!(B::BlockDiagonal{T, V}, v::V, p::Integer) where {T, V}
    if blocksize(B, p) != size(v)
        throw(DimensionMismatch(
            "Cannot set block of size $(blocksize(B, p)) to block of size $(size(v))."
        ))
    end
    return blocks(B)[p] = v
end

function setblock!(B::BlockDiagonal{T, V}, v::V, p::Int, q::Int) where {T, V}
    p == q || throw(ArgumentError("Cannot set off-diagonal block ($p, $q) to non-zero value."))
    return setblock!(B, v, p)
end

## Base
function Base.Matrix(B::BlockDiagonal{T}) where {T}
    A = zeros(T, size(B))

    nrows = size.(B.blocks, 1)
    ncols = size.(B.blocks, 2)
    row_idxs = cumsum(nrows) .- nrows .+ 1
    col_idxs = cumsum(ncols) .- ncols .+ 1

    for n in eachindex(blocks(B))
        block_rows = row_idxs[n]:(row_idxs[n] + nrows[n] - 1)
        block_cols = col_idxs[n]:(col_idxs[n] + ncols[n] - 1)
        A[block_rows, block_cols] .= blocks(B)[n]
    end
    return A
end

function ChainRulesCore.rrule(::Type{<:Base.Matrix}, B::T) where {T<:BlockDiagonal}
    nrows = size.(B.blocks, 1)
    ncols = size.(B.blocks, 2)
    function Matrix_pullback(Δ::Matrix)
        row_idxs = cumsum(nrows) .- nrows .+ 1
        col_idxs = cumsum(ncols) .- ncols .+ 1

        Δblocks = map(eachindex(nrows)) do n
            block_rows = row_idxs[n]:(row_idxs[n] + nrows[n] - 1)
            block_cols = col_idxs[n]:(col_idxs[n] + ncols[n] - 1)
            return Δ[block_rows, block_cols]
        end
        return (NO_FIELDS, Composite{T}(blocks=Δblocks))
    end
    return Matrix(B), Matrix_pullback
end

Base.size(B::BlockDiagonal) = sum(first∘size, blocks(B)), sum(last∘size, blocks(B))
Base.similar(B::BlockDiagonal) = BlockDiagonal(map(similar, blocks(B)))
Base.parent(B::BlockDiagonal) = B.blocks

@propagate_inbounds function Base.setindex!(B::BlockDiagonal, v, i::Integer, j::Integer)
    p, i_, j_ = _block_indices(B, i, j)
    if p > 0
        @inbounds getblock(B, p)[i_, end+j_] = v
    elseif !iszero(v)
        throw(ArgumentError(
            "Cannot set entry ($i, $j) in off-diagonal-block to nonzero value $v."
        ))
    end
    return v
end

@propagate_inbounds function Base.getindex(B::BlockDiagonal{T}, i::Integer, j::Integer) where T
    p, i, j = _block_indices(B, i, j)
    # if not in on-diagonal block `p` then value at `i, j` must be zero
    @inbounds return p > 0 ? getblock(B, p)[i, end + j] : zero(T)
end

# Transform indices `i, j` (identifying entry `Matrix(B)[i, j]`) into indices `p, i, j` such
# that the same entry is available via `getblock(B, p)[i, end+j]`; `p = -1` if no such `p`.
function _block_indices(B::BlockDiagonal, i::Integer, j::Integer)
    all((0, 0) .< (i, j) .<= size(B)) || throw(BoundsError(B, (i, j)))
    nrows = size.(blocks(B), 1)
    ncols = size.(blocks(B), 2)
    # find the on-diagonal block `p` in column `j`
    p = 0
    while j > 0
        p += 1
        j -= ncols[p]
    end
    i -= sum(nrows[1:(p-1)])
    # if row `i` outside of block `p`, set `p` to place-holder value `-1`
    if i <= 0 || i > nrows[p]
        p = -1
    end
    return p, i, j
end

function Base.copy!(dest::BlockDiagonal, src::BlockDiagonal)
    isequal_blocksizes(dest, src) || throw(DimensionMismatch("dest and src have different block sizes"))
    for i in eachindex(blocks(dest))
        @inbounds copyto!(dest.blocks[i], src.blocks[i])
    end
    return dest
end


text:
module utils
using Random

#Import existing push!() and pop!() method definitions to qualify our push!() and pop()! methods for export.
import Base.push!,
        Base.pop!,
        Base.iterate,
        Base.argmin,
        Base.argmax,
        Base.length,
        Base.delete!;

export if_, Queue, FIFOQueue, Stack, PQueue, push!, pop!, extend!, delete!,
        iterate, length,
        MemoizedFunction, eval_memoized_function,
        AbstractProblem,
        argmin, argmax, argmin_random_tie, argmax_random_tie,
        weighted_sampler, weighted_sample_with_replacement,
        distance, distance2,
        RandomDeviceInstance,
        isfunction, removeall,
        normalize_probability_distribution,
        mode, sigmoid, sigmoid_derivative,
        combinations, iterable_cartesian_product,
        weighted_choice;

function if_(boolean_expression::Bool, ans1::Any, ans2::Any)
    if (boolean_expression)
        return ans1;
    else
        return ans2;
    end
end

function distance(p1::Tuple{Number, Number}, p2::Tuple{Number, Number})
    return sqrt(((Float64(p1[1]) - Float64(p2[1]))^2) + ((Float64(p1[2]) - Float64(p2[2]))^2));
end

function distance2(p1::Tuple{Number, Number}, p2::Tuple{Number, Number})
    return (Float64(p1[1]) - Float64(p2[1]))^2 + (Float64(p1[2]) - Float64(p2[2]))^2;
end

function null_index(v::AbstractVector)
    local i::Int64 = 0;
    for element in v
        i = i + 1;
        if (element === nothing)
            return i;
        end
    end
    return -1;          #couldn't find the item in the array
end

function index(v::Array{T, 1}, item::T) where {T <: Any}
    local i::Int64 = 0;
    for element in v
        i = i + 1;
        if (element == item)
            return i;
        end
    end
    return -1;          #couldn't find the item in the array
end

function turn_heading(heading::Tuple{Any, Any}, inc::Int64)
    local o = [(1, 0), (0, 1), (-1, 0), (0, -1)];
    # 4 (for negative increments) - 1 (adjust index) = 3 offset
    return o[((index(o, heading) + inc + 3) % length(o)) + 1];
end

function vector_add_tuples(a::Tuple, b::Tuple)
    return map(+, a, b);
end

abstract type AbstractProblem end;

RandomDeviceInstance = RandomDevice();

#=

    Define a Queue as an abstract DataType.

    FIFOQueue, PriorityQueue, Stack are implementations of the Queue DataType.

=#

abstract type Queue end;

#=

    Stack is a Last In First Out (LIFO) Queue implementation.

=#
struct Stack <: Queue
    array::Array{Any, 1}

    function Stack()
        return new(Array{Any, 1}());
    end
end

# Map length function calls to underlying array in Stack
length(s::Stack) = length(s.array);
isempty(s::Stack) = isempty(s.array);

# Map iterator function calls to underlying array in Stack
iterate(s::Stack) = iterate(s.array);
iterate(s::Stack, i) = iterate(s.array, i);

#=

    FIFOQueue is a First In First Out (FIFO) Queue implementation.

=#
struct FIFOQueue <: Queue
    array::Array{Any, 1}

    function FIFOQueue()
        return new(Array{Any, 1}());
    end
end

# Map length function calls to underlying array in FIFOQueue
length(fq::FIFOQueue) = length(fq.array);
isempty(fq::FIFOQueue) = isempty(fq.array);

# Map iterator function calls to underlying array in FIFOQueue
iterate(fq::FIFOQueue) = iterate(fq.array);
iterate(fq::FIFOQueue, i) = iterate(fq.array, i);

#=

    PQueue is a Priority Queue implementation.

    The array must consist of Tuple{Any, Any} such that,

        -the first element is the priority of the item.

        -the second element is the item.

=#
struct PQueue <: Queue
    array::Array{Tuple{Any, Any}, 1}
    order::Base.Order.Ordering

    function PQueue(;order::Base.Order.Ordering=Base.Order.Forward)
        return new(Array{Tuple{Any, Any}, 1}(), order);
    end
end

# Map length function calls to underlying array in PQueue
length(pq::PQueue) = length(pq.array);
isempty(pq::PQueue) = isempty(pq.array);

# Map iterator function calls to underlying array in PQueue
iterate(pq::PQueue) = iterate(pq.array);
iterate(pq::PQueue, i) = iterate(pq.array, i);

#=

    MemoizedFunction is a DataType that wraps the original function (.f) with a dictionary

    (.values) containing the previous given arguments as keys to their computed values.

=#

struct MemoizedFunction
    f::Function     #original function
    values::Dict{Tuple{Vararg}, Any}

    function MemoizedFunction(f::Function)
        return new(f, Dict{Tuple{Vararg}, Any}());
    end
end

function eval_memoized_function(mf::MemoizedFunction, args::Vararg{Any})
    if (haskey(mf.values, args))
        return mf.values[args];
    else
        mf.values[args] = mf.f(args...);
        return mf.values[args];
    end
end

#=

    Define method definitions of push!(), pop()!, and extend()! for Queue implementations.

=#

"""
    push!(s::Stack, i::Any)

Push the given item 'i' to the end of the collection.
"""
function push!(s::Stack, i::Any)
    push!(s.array, i);
    nothing;
end

"""
    push!(fq::FIFOQueue, i::Any)

Push the given item 'i' to the end of the collection.
"""
function push!(fq::FIFOQueue, i::Any)
    push!(fq.array, i);
    nothing;
end

"""
    pop!(s::Stack)

Delete the last item of the collection and return the deleted item.
"""
function pop!(s::Stack)
    return pop!(s.array);
end

"""
    pop!(fq::FIFOQueue)

Delete the first item of the collection and return the deleted item.
"""
function pop!(fq::FIFOQueue)
    return popfirst!(fq.array);
end

"""
    extend!(s1::Stack, s2::Queue)
    extend!(s1::Stack, s2::AbstractVector)

Push item(s) of s2 to the end of s1.
"""
function extend!(s1::Stack, s2::T) where {T <: Queue}
    if (!(typeof(s2) <: PQueue))
        for e in s2.array
            push!(s1, e);
        end
    else
        for e in s2.array
            push!(s1, getindex(e, 2));
        end
    end
    nothing;
end

function extend!(s1::Stack, s2::AbstractVector)
    for e in s2
        push!(s1, e);
    end
    nothing;
end

"""
    extend!(fq1::FIFOQueue, fq2::Queue)
    extend!(fq1::FIFOQueue, fq2::AbstractVector)

Push item(s) of fq2 to the end of fq1.
"""
function extend!(fq1::FIFOQueue, fq2::T) where {T <: Queue}
    if (!(typeof(fq2) <: PQueue))
        for e in fq2.array
            push!(fq1, e);
        end
    else
        for e in fq2.array
            push!(fq1, getindex(e, 2));
        end
    end
    nothing;
end

function extend!(fq1::FIFOQueue, fq2::AbstractVector)
    for e in fq2
        push!(fq1, e);
    end
    nothing;
end

# Modified sorted binary search for array of tuples
#   https://github.com/JuliaLang/julia/blob/master/base/sort.jl
#       searchsortedfirst(), searchsortedlast(), and searchsorted()
#
#       These 3 functions were renamed to avoid confusion.
#
#
# Base.Order.Forward will make the PQueue ordered by minimums.
# Base.Order.Reverse will make the PQueue ordered by maximums.
function bisearchfirst(v::Array{T, 1}, x::T, lo::Int, hi::Int, o::Base.Sort.Ordering) where {T <: Tuple{Any, Any}}
    lo = lo-1;
    hi = hi+1;
    @inbounds while (lo < hi-1)
        m = (lo+hi) >>> 1
        if (Base.Order.lt(o, getindex(v[m], 1), getindex(x, 1)))
            lo = m;
        else
            hi = m;
        end
    end
    return hi;
end

function bisearchlast(v::Array{T, 1}, x::T, lo::Int, hi::Int, o::Base.Sort.Ordering) where {T <: Tuple{Any, Any}}
    lo = lo-1;
    hi = hi+1;
    @inbounds while (lo < hi-1)
        m = (lo+hi) >>> 1;
        if (Base.Order.lt(o, getindex(x, 1), getindex(v[m], 1)))
            hi = m;
        else
            lo = m;
        end
    end
    return lo;
end

function bisearch(v::Array{T, 1}, x::T, ilo::Int, ihi::Int, o::Base.Sort.Ordering) where {T <: Tuple{Any, Any}}
    lo = ilo-1;
    hi = ihi+1;
    @inbounds while (lo < hi-1)
        m = (lo+hi) >>> 1;
        if (Base.Order.lt(o, getindex(v[m], 1), getindex(x, 1)))
            lo = m;
        elseif (Base.Order.lt(o, getindex(x, 1), getindex(v[m], 1)))
            hi = m;
        else
            a = bisearchfirst(v, x, max(lo, ilo), m, o)
            b = bisearchlast(v, x, m, min(hi, ihi), o)
            return a : b;
        end
    end
    return (lo + 1) : (hi - 1);
end

"""
    push!(pq::PQueue, i::Tuple{Any, Tuple})

Push the given item 'i' to the index after existing entries with the same priority as getitem(i, 1).
"""
function push!(pq::PQueue, item::Tuple{Any, Any})
    bsi = bisearch(pq.array, item, 1, length(pq), pq.order);

    if (pq.order == Base.Order.Forward)
        insert!(pq.array, bsi.stop + 1, item);
    else
        insert!(pq.array, bsi.start, item);
    end
    nothing;
end

function push!(pq::PQueue, item::Any, mf::MemoizedFunction)
    local item_tuple = (eval_memoized_function(mf, item), item);
    bsi = bisearch(pq.array, item_tuple, 1, length(pq), pq.order);

    if (pq.order == Base.Order.Forward)
        insert!(pq.array, bsi.stop + 1, item_tuple);
    else
        insert!(pq.array, bsi.start, item_tuple);
    end
    nothing;
end

function push!(pq::PQueue, item::Any, mf::Function)
    local item_tuple = (mf(item), item);
    bsi = bisearch(pq.array, item_tuple, 1, length(pq), pq.order);

    if (pq.order == Base.Order.Forward)
        insert!(pq.array, bsi.stop + 1, item_tuple);
    else
        insert!(pq.array, bsi.start, item_tuple);
    end
    nothing;
end

function push!(pq::PQueue, item::Any, mf::MemoizedFunction, problem::T) where {T <: AbstractProblem}
    local item_tuple = (eval_memoized_function(mf, problem, item), item);
    bsi = bisearch(pq.array, item_tuple, 1, length(pq), pq.order);

    if (pq.order == Base.Order.Forward)
        insert!(pq.array, bsi.stop + 1, item_tuple);
    else
        if (bsi.start != 0)
            insert!(pq.array, bsi.start, item_tuple);
        else
            insert!(pq.array, 1, item_tuple);
        end
    end
    nothing;
end

"""
    pop!(pq::PQueue)

Delete the lowest/highest item based on ordering of the collection and return the deleted item.
"""
function pop!(pq::PQueue)
    return getindex(popfirst!(pq.array), 2);   #return lowest/highest priority tuple by pq.order
end

"""
    extend!(pq1::PQueue, pq2::Queue, pv::Function)
    extend!(pq1::PQueue, pq2::AbstractVector, pv::Function)

Push item(s) of pq2 to pq1 by the priority of the item(s) returned by pv().
"""
function extend!(pq1::PQueue, pq2::T, pv::Function) where {T <: Queue}
    if (!(typeof(pq2) <: PQueue))
        for e in pq2.array
            push!(pq1, (pv(e), e));
        end
    else
        for e in pq2.array
            push!(pq1, (pv(getindex(e, 2)), getindex(e, 2)));
        end
    end
    nothing;
end

#pv - function that returns priority value
function extend!(pq1::PQueue, pq2::AbstractVector, pv::Function)
    for e in pq2
        push!(pq1, (pv(e), e));
    end
    nothing;
end

#mpv - function that returns memoized priority value
function extend!(pq1::PQueue, pq2::AbstractVector, mpv::MemoizedFunction)
    for e in pq2
        push!(pq1, (eval_memoized_function(mpv, e), e));
    end
    nothing;
end

function extend!(pq1::PQueue, pq2::AbstractVector, mpv::MemoizedFunction, problem::T) where {T <: AbstractProblem}
    for e in pq2
        push!(pq1, (eval_memoized_function(mpv, problem, e), e));
    end
    nothing;
end

"""
    delete!(pq::PQueue, item::Any)

Remove the item if it already exists in pq.array.
"""
function delete!(pq::PQueue, item::Any)
    for (i, entry) in enumerate(pq.array)
        if (item == getindex(entry, 2))
            deleteat!(pq.array, i);
            return nothing;
        end
    end
    return nothing;
end

function removeall(v::String, item)
    return replace(v, item, "");
end

function removeall(v::AbstractVector, item)
    return collect(x for x in v if (x != item));
end

"""
    weighted_sample_with_replacement(seq, weights, n)

Return an array of 'n' elements that are chosen from 'seq' at random with replacement, with
the probability of picking each element based on its corresponding weight in 'weights'.
"""
function weighted_sample_with_replacement(seq::T1, weights::T2, n::Int64) where {T1 <: Vector, T2 <: Vector}
    local sample = weighted_sampler(seq, weights);
    return collect(sample() for i in 1:n);
end

function weighted_sample_with_replacement(seq::String, weights::T, n::Int64) where {T <: Vector}
    local sample = weighted_sampler(seq, weights);
    return collect(sample() for i in 1:n);
end

"""
    weighted_sampler(seq, weights)

Return a random sample function that chooses an element from 'seq' based on its corresponding
weight in 'weight'.
"""
function weighted_sampler(seq::T1, weights::T2) where {T1 <: Vector, T2 <: Vector}
    local totals::Array{Float64, 1} = Array{Float64, 1}();
    for w in weights
        if (length(totals) != 0)
            push!(totals, (w + totals[length(totals)]));
        else
            push!(totals, w);
        end
    end
    return (function(;sequence=seq, totals_array=totals)
                element = rand(RandomDeviceInstance)*totals_array[end];
                bsi = searchsorted(totals_array, element);
                if (bsi.stop == length(seq))  # Prevent indices out of bounds.
                    return seq[bsi.stop];
                else
                    return seq[bsi.stop + 1];
                end
            end);
end

function weighted_sampler(seq::String, weights::T) where {T <: Vector}
    local totals = Array{Any, 1}();
    for w in weights
        if (length(totals) != 0)
            push!(totals, (w + totals[length(totals)]));
        else
            push!(totals, w);
        end
    end
    return (function(;sequence=seq, totals_array=totals)
                bsi = searchsorted(totals_array,
                                (rand(RandomDeviceInstance)*totals_array[length(totals_array)]),
                                1,
                                length(totals_array),
                                Base.Order.Forward);
                if (bsi.stop == length(seq))  # Prevent indices out of bounds.
                    return seq[bsi.stop];
                else
                    return seq[bsi.stop + 1];
                end
            end);
end

"""
    argmin(seq, fn)

Applies fn() to each element in seq and returns the element that has the lowest fn() value. argmin()
is similar to mapreduce(fn, min, seq) in computing the best score, but returns the corresponding element.
"""
function argmin(seq::T, fn::Function) where {T <: Vector}
    local best_element = seq[1];
    local best_score = fn(best_element);
    for element in seq
        element_score = fn(element);
        if (element_score < best_score)
            best_element = element;
            best_score = element_score;
        end
    end
    return best_element;
end

function argmin_random_tie(seq::T, fn::Function) where {T <: Vector}
    local best_score = fn(seq[1]);
    local n::Int64 = 0;
    local best_element = seq[1];
    for element in seq
        element_score = fn(element);
        if (element_score < best_score)
            best_element = element;
            best_score = element_score;
        elseif (element_score == best_score)
            n = n + 1;
            if (rand(RandomDeviceInstance, 1:n) == 1)
                best_element = element;
            end
        end
    end
    return best_element;
end

"""
    argmax(seq, fn)

Applies fn() to each element in seq and returns the element that has the highest fn() value. argmax()
is similar to mapreduce(fn, max, seq) in computing the best score, but returns the corresponding element.
"""
function argmax(seq::T, fn::Function) where {T <: Vector}
    local best_element = seq[1];
    local best_score = fn(best_element);
    for element in seq
        element_score = fn(element);
        if (element_score > best_score)
            best_element = element;
            best_score = element_score;
        end
    end
    return best_element;
end

function argmax_random_tie(seq::T, fn::Function) where {T <: Vector}
    local best_score = fn(seq[1]);
    local n::Int64 = 1;
    local best_element = seq[1];
    for element in seq
        element_score = fn(element);
        if (element_score > best_score)
            best_element = element;
            best_score = element_score;
        elseif (element_score == best_score)
            n = n + 1;
            if (rand(RandomDeviceInstance, 1:n) == 1)
                best_element = element;
            end
        end
    end
    return best_element;
end

"""
    isfunction(var)

Check if 'var' is callable as a function.
"""
function isfunction(var)
    return (typeof(var) <: Function);
end

"""
    normalize_probability_distribution(d)

Return a collection such that each value is the corresponding value in 'd' divided
by the sum of all values in 'd'.
"""
function normalize_probability_distribution(d::Dict)
    local total::Float64 = sum(values(d));
    for key in keys(d)
        d[key] = d[key] / total;
        if (!(0.0 <= d[key] <= 1.0))
            error("normalize_probability_distribution(): ", d[key], " is not a valid probability.");
        end
    end
    return dist;
end

function normalize_probability_distribution(d::AbstractVector)
    local total::Float64 = sum(d);
    return collect((i / total) for i in d);
end

function mode_reverse_isless(p1::Tuple, p2::Tuple)
    return (p1[2] > p2[2]);
end

function mode(v::AbstractVector)
    local sorted::AbstractVector = sort!(collect((i, count(j->(j == i), v)) for i in Set(v)),
                                        lt=mode_reverse_isless);
    if (length(sorted) == 0)
        error("mode(): There is no mode for an empty array!");
    else
        return getindex(getindex(sorted, 1), 1);
    end
end

function mode(iter::Base.Generator)
    local sorted::AbstractVector = sort!(collect((i, count(j->(j == i), iter)) for i in Set(iter)),
                                        lt=mode_reverse_isless);
    if (length(sorted) == 0)
        error("mode(): There is no mode for an empty array!");
    else
        return getindex(getindex(sorted, 1), 1);
    end
end

"""
    sigmoid(x::Number)

Return the activation value of 'x' by using a sigmoid function 'S(x)' as the activation function.
"""
function sigmoid(x::Number)
    return (Float64(1)/(Float64(1) + exp(-x)));
end

"""
    sigmoid_derivative(val::Number)

Return the derivative of the sigmoid function 'S(x)', where x = 'val'.
"""
function sigmoid_derivative(val::Number)
    return (Float64(val) * (Float64(1) - Float64(val)));
end

# The combinations() function below was adapted from
# https://github.com/JuliaMath/Combinatorics.jl/blob/master/src/combinations.jl

"""
    combinations(array::AbstractVector, l::Integer)
    combinations(tuple::Tuple, l::Integer)
    combinations(set::Set, l::Integer)

Return the 'l' length subsequences of the elements in the given collection of items.
"""
function combinations(array::AbstractVector, l::Integer)
    local indices::AbstractVector = collect(1:l);
    local visited::Tuple = ();
    local current_combination::AbstractVector;
    if (l == 0)
        return ([],);
    end

    if (binomial(length(array), l) > 0)
        while (indices[1] <= length(array) - l + 1)
            current_combination = collect(array[subseq_i] for subseq_i in indices);
            visited = (visited..., current_combination);
            indices = copy(indices);
            for i in reverse(1:length(indices))
                indices[i] = indices[i] + 1;
                if (indices[i] > (length(array) - (length(indices) - i)))
                    continue;
                end
                for j in (i + 1):lastindex(indices)
                    indices[j] = indices[j - 1] + 1;
                end
                break;
            end
        end
        return visited;
    else
        return visited;
    end
end

function combinations(tuple::Tuple, l::Integer)
    local indices::AbstractVector = collect(1:l);
    local visited::Tuple = ();
    local current_combination::AbstractVector;
    if (l == 0)
        return ([],);
    end

    if (binomial(length(tuple), l) > 0)
        while (indices[1] <= length(tuple) - l + 1)
            current_combination = collect(tuple[subseq_i] for subseq_i in indices);
            visited = (visited..., current_combination);
            indices = copy(indices);
            for i in reverse(1:length(indices))
                indices[i] = indices[i] + 1;
                if (indices[i] > (length(tuple) - (length(indices) - i)))
                    continue;
                end
                for j in (i + 1):lastindex(indices)
                    indices[j] = indices[j - 1] + 1;
                end
                break;
            end
        end
        return visited;
    else
        return visited;
    end
end

function combinations(set::Set, l::Integer)
    local array::AbstractVector = collect(set);
    local indices::AbstractVector = collect(1:l);
    local visited::Tuple = ();
    local current_combination::AbstractVector;
    if (l == 0)
        return ([],);
    end

    if (binomial(length(array), l) > 0)
        while (indices[1] <= length(array) - l + 1)
            current_combination = collect(array[subseq_i] for subseq_i in indices);
            visited = (visited..., current_combination);
            indices = copy(indices);
            for i in reverse(1:length(indices))
                indices[i] = indices[i] + 1;
                if (indices[i] > (length(array) - (length(indices) - i)))
                    continue;
                end
                for j in (i + 1):lastindex(indices)
                    indices[j] = indices[j - 1] + 1;
                end
                break;
            end
        end
        return visited;
    else
        return visited;
    end
end

function iterable_cartesian_product(iterable_items::AbstractVector, current_index::Int64, current_permutation::AbstractVector, product_array::AbstractVector)
    if (current_index == length(iterable_items))
        push!(product_array, current_permutation);
    elseif (current_index > length(iterable_items))
        error("iterable_cartesian_product(): The current index ", current_index, " exceeds the length of the given array!");
    else
        if ((typeof(iterable_items[current_index + 1]) <: Vector)
            || (typeof(iterable_items[current_index + 1]) <: Tuple)
            || (typeof(iterable_items[current_index + 1]) <: Set))
            for item in iterable_items[current_index + 1]
                iterable_cartesian_product(iterable_items, (current_index + 1), vcat(current_permutation, item), product_array);
            end
        else
            error("iterable_cartesian_product(): iterable_items[", current_index + 1, "] is not iterable!");
        end
    end
end

function iterable_cartesian_product(iterable_items::Tuple, current_index::Int64, current_permutation::AbstractVector, product_array::AbstractVector)
    if (current_index == length(iterable_items))
        push!(product_array, current_permutation);
    elseif (current_index > length(iterable_items))
        error("iterable_cartesian_product(): The current index ", current_index, " exceeds the length of the given array!");
    else
        if ((typeof(iterable_items[current_index + 1]) <: Vector)
            || (typeof(iterable_items[current_index + 1]) <: Tuple)
            || (typeof(iterable_items[current_index + 1]) <: Set))
            for item in iterable_items[current_index + 1]
                iterable_cartesian_product(iterable_items, (current_index + 1), vcat(current_permutation, item), product_array);
            end
        else
            error("iterable_cartesian_product(): iterable_items[", current_index + 1, "] is not iterable!");
        end
    end
end

"""
    iterable_cartesian_product(iterable_items)

Return the cartesian product of given items in 'iterable_items' as an array.
"""
function iterable_cartesian_product(iterable_items::AbstractVector)
    local product_array::AbstractVector = [];
    iterable_cartesian_product(iterable_items, 0, [], product_array);
    return product_array;
end

function iterable_cartesian_product(iterable_items::Tuple)
    local product_array::AbstractVector = [];
    iterable_cartesian_product(iterable_items, 0, [], product_array);
    return product_array;
end

function iterable_cartesian_product(iterable_items::Set)
    local product_array::AbstractVector = [];
    iterable_cartesian_product((iterable_items...), 0, [], product_array);
    return product_array;
end

"""
    weighted_choice(choices::AbstractVector)

Return an element from the given array 'choices' based on the choice and its corresponding weight.
"""
function weighted_choice(choices::AbstractVector)
    local total::Float64 = sum(collect(choice[2] for choice in choices));
    local r::Float64 = rand(RandomDeviceInstance) * total;
    local upto::Float64 = 0.0;
    for (choice, weight) in choices
        if ((upto + weight) >= r)
            return (choice, weight);
        end
        upto = upto + weight;
    end
end

end

max_stars_repo_name: mikhail-j/aima-julia-aimacode-fork

max_issues_repo_name: mikhail-j/aima-julia-aimacode-fork

max_forks_repo_name: mikhail-j/aima-julia-aimacode-fork