text:
module ExpressionExplorer
export compute_symbolreferences, try_compute_symbolreferences, compute_usings, SymbolsState, FuncName, join_funcname_parts

import ..PlutoRunner
import Markdown
import Base: union, union!, ==, push!

###
# TWO STATE OBJECTS
###

# TODO: use GlobalRef instead
FuncName = Array{Symbol,1}

"SymbolsState trickles _down_ the ASTree: it carries referenced and defined variables from endpoints down to the root."
mutable struct SymbolsState
    references::Set{Symbol}
    assignments::Set{Symbol}
    funccalls::Set{FuncName}
    funcdefs::Dict{FuncName,SymbolsState}
end

SymbolsState(references, assignments, funccalls) = SymbolsState(references, assignments, funccalls, Dict{FuncName,SymbolsState}())
SymbolsState(references, assignments) = SymbolsState(references, assignments, Set{Symbol}())
SymbolsState() = SymbolsState(Set{Symbol}(), Set{Symbol}())

function Base.show(io::IO, s::SymbolsState)
    print(io, "SymbolsState([")
    join(io, s.references, ", ")
    print(io, "], [")
    join(io, s.assignments, ", ")
    print(io, "], [")
    join(io, s.funccalls, ", ")
    print(io, "], [")
    if isempty(s.funcdefs)
        print(io, "]")
    else
        println(io)
        for (k, v) in s.funcdefs
            print(io, "    ", k, ": ", v)
            println(io)
        end
        print(io, "]")
    end
    print(io, ")")
end

"ScopeState moves _up_ the ASTree: it carries scope information up towards the endpoints."
mutable struct ScopeState
    inglobalscope::Bool
    exposedglobals::Set{Symbol}
    hiddenglobals::Set{Symbol}
    definedfuncs::Set{Symbol}
end

# The `union` and `union!` overloads define how two `SymbolsState`s or two `ScopeState`s are combined.

function union(a::Dict{FuncName,SymbolsState}, bs::Dict{FuncName,SymbolsState}...)
    union!(Dict{FuncName,SymbolsState}(), a, bs...)
end

function union!(a::Dict{FuncName,SymbolsState}, bs::Dict{FuncName,SymbolsState}...)
    for b in bs
        for (k, v) in b
            if haskey(a, k)
                a[k] = union!(a[k], v)
            else
                a[k] = v
            end
        end
        a
    end
    return a
end

function union(a::SymbolsState, b::SymbolsState)
    SymbolsState(a.references ∪ b.references, a.assignments ∪ b.assignments, a.funccalls ∪ b.funccalls, a.funcdefs ∪ b.funcdefs)
end

function union!(a::SymbolsState, bs::SymbolsState...)
    union!(a.references, (b.references for b in bs)...)
    union!(a.assignments, (b.assignments for b in bs)...)
    union!(a.funccalls, (b.funccalls for b in bs)...)
    union!(a.funcdefs, (b.funcdefs for b in bs)...)
    return a
end

function union!(a::Tuple{FuncName,SymbolsState}, bs::Tuple{FuncName,SymbolsState}...)
    a[1], union!(a[2], (b[2] for b in bs)...)
end

function union(a::ScopeState, b::ScopeState)
    SymbolsState(a.inglobalscope && b.inglobalscope, a.exposedglobals ∪ b.exposedglobals, a.hiddenglobals ∪ b.hiddenglobals)
end

function union!(a::ScopeState, bs::ScopeState...)
    a.inglobalscope &= all((b.inglobalscope for b in bs)...)
    union!(a.exposedglobals, (b.exposedglobals for b in bs)...)
    union!(a.hiddenglobals, (b.hiddenglobals for b in bs)...)
    union!(a.definedfuncs, (b.definedfuncs for b in bs)...)
    return a
end

function ==(a::SymbolsState, b::SymbolsState)
    a.references == b.references && a.assignments == b.assignments && a.funccalls == b.funccalls && a.funcdefs == b.funcdefs
end

Base.push!(x::Set) = x

###
# HELPER FUNCTIONS
###

# from the source code: https://github.com/JuliaLang/julia/blob/master/src/julia-parser.scm#L9
const modifiers = [:(+=), :(-=), :(*=), :(/=), :(//=), :(^=), :(÷=), :(%=), :(<<=), :(>>=), :(>>>=), :(&=), :(⊻=), :(≔), :(⩴), :(≕)]
const modifiers_dotprefixed = [Symbol('.' * String(m)) for m in modifiers]

function will_assign_global(assignee::Symbol, scopestate::ScopeState)::Bool
    (scopestate.inglobalscope || assignee ∈ scopestate.exposedglobals) && (assignee ∉ scopestate.hiddenglobals || assignee ∈ scopestate.definedfuncs)
end

function will_assign_global(assignee::Array{Symbol,1}, scopestate::ScopeState)::Bool
    if length(assignee) == 0
        false
    elseif length(assignee) > 1
        scopestate.inglobalscope
    else
        will_assign_global(assignee[1], scopestate)
    end
end

function get_global_assignees(assignee_exprs, scopestate::ScopeState)::Set{Symbol}
    global_assignees = Set{Symbol}()
    for ae in assignee_exprs
        if isa(ae, Symbol)
            will_assign_global(ae, scopestate) && push!(global_assignees, ae)
        else
            if ae.head == :(::)
                will_assign_global(ae.args[1], scopestate) && push!(global_assignees, ae.args[1])
            else
                @warn "Unknown assignee expression" ae
            end
        end
    end
    return global_assignees
end

function get_assignees(ex::Expr)::FuncName
    if ex.head == :tuple
        # e.g. (x, y) in the ex (x, y) = (1, 23)
        union!(Symbol[], get_assignees.(ex.args)...)
        # filter(s->s isa Symbol, ex.args)
    elseif ex.head == :(::)
        # TODO: type is referenced
        Symbol[ex.args[1]]
    elseif ex.head == :ref || ex.head == :(.)
        Symbol[]
    else
        @warn "unknow use of `=`. Assignee is unrecognised." ex
        Symbol[]
    end
end

# e.g. x = 123
get_assignees(ex::Symbol) = Symbol[ex]

# When you assign to a datatype like Int, String, or anything bad like that
# e.g. 1 = 2
# This is parsable code, so we have to treat it
get_assignees(::Any) = Symbol[]

# TODO: this should return a FuncName, and use `split_funcname`.
"Turn :(A{T}) into :A."
function uncurly!(ex::Expr, scopestate::ScopeState)::Symbol
    @assert ex.head == :curly
    push!(scopestate.hiddenglobals, (a for a in ex.args[2:end] if a isa Symbol)...)
    Symbol(ex.args[1])
end

uncurly!(ex::Expr)::Symbol = ex.args[1]

uncurly!(s::Symbol, scopestate=nothing)::Symbol = s

"Turn `:(Base.Submodule.f)` into `[:Base, :Submodule, :f]` and `:f` into `[:f]`."
function split_funcname(funcname_ex::Expr)::FuncName
    if funcname_ex.head == :(.)
        vcat(split_funcname.(funcname_ex.args)...)
    else
        # a call to a function that's not a global, like calling an array element: `funcs[12]()`
        # TODO: explore symstate!
        Symbol[]
    end
end

function split_funcname(funcname_ex::QuoteNode)::FuncName
    split_funcname(funcname_ex.value)
end

function split_funcname(funcname_ex::Symbol)::FuncName
    Symbol[funcname_ex |> without_dotprefix |> without_dotsuffix]
end

# this includes GlobalRef - it's fine that we don't recognise it, because you can't assign to a globalref?
function split_funcname(::Any)::FuncName
    Symbol[]
end

"Turn :(.+) into :(+)"
function without_dotprefix(funcname::Symbol)::Symbol
    fn_str = String(funcname)
    if length(fn_str) > 0 && fn_str[1] == '.'
        Symbol(fn_str[2:end])
    else
        funcname
    end
end

"Turn :(sqrt.) into :(sqrt)"
function without_dotsuffix(funcname::Symbol)::Symbol
    fn_str = String(funcname)
    if length(fn_str) > 0 && fn_str[end] == '.'
        Symbol(fn_str[1:end - 1])
    else
        funcname
    end
end

"""Turn `Symbol[:Module, :func]` into Symbol("Module.func").

This is **not** the same as the expression `:(Module.func)`, but is used to identify the function name using a single `Symbol` (like normal variables).
This means that it is only the inverse of `ExpressionExplorer.split_funcname` iff `length(parts) ≤ 1`."""
function join_funcname_parts(parts::FuncName)::Symbol
	join(parts .|> String, ".") |> Symbol
end

###
# MAIN RECURSIVE FUNCTION
###

# Spaghetti code for a spaghetti problem 🍝

# Possible leaf: value
# Like: a = 1
# 1 is a value (Int64)
function explore!(value, scopestate::ScopeState)::SymbolsState
    # includes: LineNumberNode, Int64, String,
    return SymbolsState()
end

# Possible leaf: symbol
# Like a = x
# x is a symbol
# We handle the assignment separately, and explore!(:a, ...) will not be called.
# Therefore, this method only handles _references_, which are added to the symbolstate, depending on the scopestate.
function explore!(sym::Symbol, scopestate::ScopeState)::SymbolsState
    if sym ∈ scopestate.hiddenglobals
        SymbolsState(Set{Symbol}(), Set{Symbol}(), Set{Symbol}(), Dict{FuncName,SymbolsState}())
    else
        SymbolsState(Set([sym]), Set{Symbol}(), Set{Symbol}(), Dict{FuncName,SymbolsState}())
    end
end

# General recursive method. Is never a leaf.
# Modifies the `scopestate`.
function explore!(ex::Expr, scopestate::ScopeState)::SymbolsState
    if ex.head == :(=)
        # Does not create scope

        if ex.args[1] isa Expr && (ex.args[1].head == :call || ex.args[1].head == :where || (ex.args[1].head == :(::) && ex.args[1].args[1] isa Expr && ex.args[1].args[1].head == :call))
            # f(x, y) = x + y
            # Rewrite to:
            # function f(x, y) x + y end
            return explore!(Expr(:function, ex.args...), scopestate)
        end
        assignees = get_assignees(ex.args[1])
        val = ex.args[2]

        global_assignees = get_global_assignees(assignees, scopestate)

        symstate = innersymstate = explore!(val, scopestate)
        # If we are _not_ assigning a global variable, then this symbol hides any global definition with that name
        push!(scopestate.hiddenglobals, setdiff(assignees, global_assignees)...)
        assigneesymstate = explore!(ex.args[1], scopestate)

        push!(scopestate.hiddenglobals, global_assignees...)
        push!(symstate.assignments, global_assignees...)
        push!(symstate.references, setdiff(assigneesymstate.references, global_assignees)...)

        return symstate
    elseif ex.head in modifiers
        # We change: a[1] += 123
        # to:        a[1] = a[1] + 123
        # We transform the modifier back to its operator
        # for when users redefine the + function

        operator = let
            s = string(ex.head)
            Symbol(s[1:prevind(s, lastindex(s))])
        end
        expanded_expr = Expr(:(=), ex.args[1], Expr(:call, operator, ex.args[1], ex.args[2]))
        return explore!(expanded_expr, scopestate)
    elseif ex.head in modifiers_dotprefixed
        # We change: a[1] .+= 123
        # to:        a[1] .= a[1] + 123

        operator = Symbol(string(ex.head)[2:end - 1])
        expanded_expr = Expr(:(.=), ex.args[1], Expr(:call, operator, ex.args[1], ex.args[2]))
        return explore!(expanded_expr, scopestate)
    elseif ex.head == :let || ex.head == :for || ex.head == :while
        # Creates local scope

        # Because we are entering a new scope, we create a copy of the current scope state, and run it through the expressions.
        innerscopestate = deepcopy(scopestate)
        innerscopestate.inglobalscope = false

        return mapfoldl(a -> explore!(a, innerscopestate), union!, ex.args, init=SymbolsState())
    elseif ex.head == :macrocall
        # Does not create sccope

        funcname = ex.args[1] |> split_funcname
        # Macros can transform the expression into anything - the best way to treat them is to macroexpand
        # the problem is that the macro is only available on the worker process, see https://github.com/fonsp/Pluto.jl/issues/196
        if (length(funcname) == 1 || (length(funcname) >= 2 && funcname[1] == :Base))
            if funcname[end] == Symbol("@md_str") || funcname[end] == Symbol("@bind") || funcname[end] == Symbol("@gensym") || funcname[end] == Symbol("@kwdef")
                # we macroexpand these, and recurse
                expanded = macroexpand(PlutoRunner, ex; recursive=false)
                return explore!(Expr(:call, ex.args[1], expanded), scopestate)
            elseif funcname[end] == Symbol("@enum")
                # we could do macroexpand, but the expanded macro defines typemin and typemax methods for the new enum type, and because of
                # https://github.com/fonsp/Pluto.jl/issues/177
                # this would mean that you can only define one enum per notebook :(
                syms = filter(x -> x isa Symbol, ex.args[2:end])
                rest = setdiff(ex.args[2:end], syms)

                return mapfoldl(a -> explore!(a, scopestate), union!, rest, init=SymbolsState(Set{Symbol}(), Set{Symbol}(syms), Set{FuncName}([[Symbol("@enum")]])))
            end
        end

        return explore!(Expr(:call, ex.args...), scopestate)
    elseif ex.head == :call
        # Does not create scope

        funcname = ex.args[1] |> split_funcname
        symstate = if length(funcname) == 0
            explore!(ex.args[1], scopestate)
        elseif length(funcname) == 1
            if funcname[1] ∈ scopestate.hiddenglobals
                SymbolsState()
            else
            SymbolsState(Set{Symbol}(), Set{Symbol}(), Set{FuncName}([funcname]))
            end
        else
            SymbolsState(Set{Symbol}([funcname[end - 1]]), Set{Symbol}(), Set{FuncName}([funcname]))
        end
        # Explore code inside function arguments:
        union!(symstate, explore!(Expr(:block, ex.args[2:end]...), scopestate))
        return symstate
    elseif ex.head == :kw
        return explore!(ex.args[2], scopestate)
    elseif ex.head == :struct
        # Creates local scope

        structname = ex.args[2]
        structfields = ex.args[3].args

        equiv_func = Expr(:function, Expr(:call, structname, structfields...), Expr(:block, nothing))

        return explore!(equiv_func, scopestate)
    elseif ex.head == :generator
        # Creates local scope

        # In a `generator`, a single expression is followed by the iterator assignments.
        # In a `for`, this expression comes at the end.

        # This is not strictly the normal form of a `for` but that's okay
        return explore!(Expr(:for, ex.args[2:end]..., ex.args[1]), scopestate)
    elseif ex.head == :function || ex.head == :abstract
        symstate = SymbolsState()
        # Creates local scope

        funcroot = ex.args[1]

        # Because we are entering a new scope, we create a copy of the current scope state, and run it through the expressions.
        innerscopestate = deepcopy(scopestate)
        innerscopestate.inglobalscope = false

        funcname, innersymstate = explore_funcdef!(funcroot, innerscopestate)

        union!(innersymstate, explore!(Expr(:block, ex.args[2:end]...), innerscopestate))

        if will_assign_global(funcname, scopestate)
            symstate.funcdefs[funcname] = innersymstate
            if length(funcname) == 1
                push!(scopestate.definedfuncs, funcname[end])
                push!(scopestate.hiddenglobals, funcname[end])
            elseif length(funcname) > 1
                push!(symstate.references, funcname[end - 1]) # reference the module of the extended function
            end
        else
            # The function is not defined globally. However, the function can still modify the global scope or reference globals, e.g.

            # let
            #     function f(x)
            #         global z = x + a
            #     end
            #     f(2)
            # end

            # so we insert the function's inner symbol state here, as if it was a `let` block.
            symstate = innersymstate
        end

        return symstate
    elseif ex.head == :(->)
        # Creates local scope

        tempname = Symbol("anon", rand(UInt64))

        # We will rewrite this to a normal function definition, with a temporary name
        funcroot = ex.args[1]
        args_ex = if funcroot isa Symbol || (funcroot isa Expr && funcroot.head == :(::))
            [funcroot]
        elseif funcroot.head == :tuple
            funcroot.args
        else
            @error "Unknown lambda type"
        end

        equiv_func = Expr(:function, Expr(:call, tempname, args_ex...), ex.args[2])

        return explore!(equiv_func, scopestate)
    elseif ex.head == :global
        # Does not create scope

        # We have one of:
        # global x;
        # global x = 1;
        # global x += 1;

        # where x can also be a tuple:
        # global a,b = 1,2

        globalisee = ex.args[1]

        if isa(globalisee, Symbol)
            push!(scopestate.exposedglobals, globalisee)
            return SymbolsState()
        elseif isa(globalisee, Expr)
            # temporarily set inglobalscope to true
            old = scopestate.inglobalscope
            scopestate.inglobalscope = true
            result = explore!(globalisee, scopestate)
            scopestate.inglobalscope = old
            return result
        else
            @error "unknow global use" ex
            return explore!(globalisee, scopestate)
        end

        return symstate
    elseif ex.head == :local
        # Does not create scope

        localisee = ex.args[1]

        if isa(localisee, Symbol)
            push!(scopestate.hiddenglobals, localisee)
            return SymbolsState()
        elseif isa(localisee, Expr) && (localisee.head == :(=) || localisee.head in modifiers)
            push!(scopestate.hiddenglobals, get_assignees(localisee.args[1])...)
            return explore!(localisee, scopestate)
        else
            @warn "unknow local use" ex
            return explore!(localisee, scopestate)
        end
    elseif ex.head == :tuple
        # Does not create scope

        # There are three (legal) cases:
        # 1. Creating a tupe:
        #   (a, b, c)

        # 2. Creating a named tuple:
        #   (a=1, b=2, c=3)

        # 3. Multiple assignments
        # a,b,c = 1,2,3
        # This parses to:
        # head = :tuple
        # args = [:a, :b, :(c=1), :2, :3]
        #
        # 🤔
        # we turn it into two expressions:
        #
        # (a, b) = (2, 3)
        # (c = 1)
        #
        # and explore those :)

        indexoffirstassignment = findfirst(a -> isa(a, Expr) && a.head == :(=), ex.args)
        if indexoffirstassignment !== nothing
            # we have one of two cases, see next `if`
            indexofsecondassignment = findnext(a -> isa(a, Expr) && a.head == :(=), ex.args, indexoffirstassignment + 1)

            if length(ex.args) == 1 || indexofsecondassignment !== nothing
                # 2.
                # we have a named tuple, e.g. (a=1, b=2)
                new_args = map(ex.args) do a
                    (a isa Expr && a.head == :(=)) ? a.args[2] : a
                end
                return explore!(Expr(:block, new_args...), scopestate)
            else
                # 3.
                # we have a tuple assignment, e.g. `a, (b, c) = [1, [2, 3]]`
                before = ex.args[1:indexoffirstassignment - 1]
                after = ex.args[indexoffirstassignment + 1:end]

                symstate_middle = explore!(ex.args[indexoffirstassignment], scopestate)
                symstate_outer = explore!(Expr(:(=), Expr(:tuple, before...), Expr(:block, after...)), scopestate)

                return union!(symstate_middle, symstate_outer)
            end
        else
            # 1.
            # good 'ol tuple
            return explore!(Expr(:block, ex.args...), scopestate)
        end
    elseif ex.head == :(.) && ex.args[2] isa Expr && ex.args[2].head == :tuple
        # pointwise function call, e.g. sqrt.(nums)
        # we rewrite to a regular call

        return explore!(Expr(:call, ex.args[1], ex.args[2].args...), scopestate)
    elseif ex.head == :using || ex.head == :import
        if scopestate.inglobalscope
            imports = if ex.args[1].head == :(:)
                ex.args[1].args[2:end]
            else
            ex.args
            end

            packagenames = map(e -> e.args[end], imports)

            return SymbolsState(Set{Symbol}(), Set{Symbol}(packagenames))
        else
            return SymbolsState(Set{Symbol}(), Set{Symbol}())
        end
    elseif ex.head == :quote
        # We ignore contents

        return SymbolsState()
    elseif ex.head == :module
        # We ignore contents; the module name is a definition

        return SymbolsState(Set{Symbol}(), Set{Symbol}([ex.args[2]]))
    else
        # fallback, includes:
        # begin, block, do, toplevel, const
        # (and hopefully much more!)

        # Does not create scope (probably)

        return mapfoldl(a -> explore!(a, scopestate), union!, ex.args, init=SymbolsState())
    end
end

"Return the function name and the SymbolsState from argument defaults. Add arguments as hidden globals to the `scopestate`.

Is also used for `struct` and `abstract`."
function explore_funcdef!(ex::Expr, scopestate::ScopeState)::Tuple{FuncName,SymbolsState}
    if ex.head == :call
        # get the function name
        name, symstate = explore_funcdef!(ex.args[1], scopestate)
        # and explore the function arguments
        return mapfoldl(a -> explore_funcdef!(a, scopestate), union!, ex.args[2:end], init=(name, symstate))

    elseif ex.head == :(::) || ex.head == :kw || ex.head == :(=)
        # recurse
        name, symstate = explore_funcdef!(ex.args[1], scopestate)
        if length(ex.args) > 1
            # use `explore!` (not `explore_funcdef!`) to explore the argument's default value - these can contain arbitrary expressions
            union!(symstate, explore!(ex.args[2], scopestate))
        end
        return name, symstate

    elseif ex.head == :where
        # function(...) where {T, S <: R, U <: A.B}
        # supertypes `R` and `A.B` are referenced
        supertypes_symstate = SymbolsState()
        for a in ex.args[2:end]
            name, inner_symstate = explore_funcdef!(a, scopestate)
            if length(name) == 1
                push!(scopestate.hiddenglobals, name[1])
            end
            union!(supertypes_symstate, inner_symstate)
        end
        # recurse
        name, symstate = explore_funcdef!(ex.args[1], scopestate)
        union!(symstate, supertypes_symstate)
        return name, symstate

    elseif ex.head == :(<:)
        # for use in `struct` and `abstract`
        name = uncurly!(ex.args[1], scopestate)
        symstate = if length(ex.args) == 1
            SymbolsState()
        else
            explore!(ex.args[2], scopestate)
        end
        return Symbol[name], symstate

    elseif ex.head == :curly
        name = uncurly!(ex, scopestate)
        return Symbol[name], SymbolsState()

    elseif ex.head == :parameters || ex.head == :tuple
        return mapfoldl(a -> explore_funcdef!(a, scopestate), union!, ex.args, init=(Symbol[], SymbolsState()))

    elseif ex.head == :(.)
        return split_funcname(ex), SymbolsState()

    else
        return Symbol[], explore!(ex, scopestate)
    end
end

function explore_funcdef!(ex::QuoteNode, scopestate::ScopeState)::Tuple{FuncName,SymbolsState}
    explore_funcdef!(ex.value, scopestate)
end

function explore_funcdef!(ex::Symbol, scopestate::ScopeState)::Tuple{FuncName,SymbolsState}
    push!(scopestate.hiddenglobals, ex)
    Symbol[ex |> without_dotprefix |> without_dotsuffix], SymbolsState()
end

function explore_funcdef!(::Any, ::ScopeState)::Tuple{FuncName,SymbolsState}
    Symbol[], SymbolsState()
end

###
# UTILITY FUNCTIONS
###

"Get the global references, assignment, function calls and function defintions inside an arbitrary expression."
function compute_symbolreferences(ex::Any)::SymbolsState
    symstate = explore!(ex, ScopeState(true, Set{Symbol}(), Set{Symbol}(), Set{Symbol}()))

    # We do something special to account for recursive functions:
    # If a function `f` calls a function `g`, and both are defined inside this cell, the reference to `g` inside the symstate of `f` will be deleted.
    # The motivitation is that normally, an assignment (or function definition) will add that symbol to a list of 'hidden globals' - any future references to that symbol will be ignored. i.e. the _local definition hides a global_.
    # In the case of functions, you can reference functions and variables that do not yet exist, and so they won't be in the list of hidden symbols when the function definition is analysed.
    # Of course, our method will fail if a referenced function is defined both inside the cell **and** in another cell. However, this will lead to a MultipleDefinitionError before anything bad happens.
    for (func, inner_symstate) in symstate.funcdefs
        inner_symstate.references = setdiff(inner_symstate.references, keys(symstate.funcdefs))
        inner_symstate.funccalls = setdiff(inner_symstate.funccalls, keys(symstate.funcdefs))
    end
    symstate
end

function try_compute_symbolreferences(ex::Any)
	try
		compute_symbolreferences(ex)
	catch e
		@error "Expression explorer failed on: " ex
		showerror(stderr, e, stacktrace(catch_backtrace()))
		SymbolsState(Set{Symbol}([:fake_reference_to_prevent_it_from_looking_like_a_text_only_cell]), Set{Symbol}())
	end
end

# TODO: this can be done during the `explore` recursion
"Get the set of `using Module` expressions that are contained in this expression."
function compute_usings(ex::Any)::Set{Expr}
    if isa(ex, Expr)
        if ex.head == :using
        Set{Expr}([ex])
        else
            union!(Set{Expr}(), compute_usings.(ex.args)...)
        end
    else
        Set{Expr}()
    end
end

"Return whether the expression is of the form `Expr(:toplevel, LineNumberNode(..), any)`."
function is_toplevel_expr(ex::Expr)::Bool
    (ex.head == :toplevel) && (length(ex.args) == 2) && (ex.args[1] isa LineNumberNode)
end

is_toplevel_expr(::Any)::Bool = false

"If the expression is a (simple) assignemnt at its root, return the assignee as `Symbol`, return `nothing` otherwise."
function get_rootassignee(ex::Expr, recurse::Bool=true)::Union{Symbol,Nothing}
    if is_toplevel_expr(ex) && recurse
        get_rootassignee(ex.args[2], false)
    elseif ex.head == :(=) && ex.args[1] isa Symbol
        ex.args[1]
    else
        nothing
    end
end

get_rootassignee(ex::Any, recuse::Bool=true)::Union{Symbol,Nothing} = nothing

end


text:
@doc raw"""
    Sphere{N} <: AbstractEmbeddedManifold{DefaultEmbeddingType}

The unit sphere manifold $𝕊^n$ represented by $n+1$-Tuples, i.e. in by
vectors in $ℝ^{n+1}$ of unit length. A sphere is an embedded manifold
and inherits the inner product of its embedding $\mathbb R^{n+1}$.
The tangent space at point p is given by

````math
T_p𝕊^n := \bigl\{ X ∈ ℝ^{n+1} : ⟨p,X⟩ = 0 \bigr \},
````

where $⟨\cdot,\cdot⟩$ denotes the inner product in the embedding $\mathbb R^{n+1}$.

This manifold is modeled as an embedded manifold to the [`Euclidean`](@ref), i.e.
several functions like the [`inner`](@ref inner(::Euclidean, ::Any...)) product and the
[`zero_tangent_vector`](@ref zero_tangent_vector(::Euclidean, ::Any...)) are inherited from the embedding.

# Constructor

    Sphere(n)

Generate $𝕊^{n} ⊂ ℝ^{n+1}$.
"""
struct Sphere{N} <: AbstractEmbeddedManifold{DefaultIsometricEmbeddingType} end

Sphere(n::Int) = Sphere{n}()

"""
    check_manifold_point(M, p; kwargs...)

Check whether `p` is a valid point on the [`Sphere`](@ref) `M`, i.e. is a vector
of length [`manifold_dimension`](@ref)`(M)+1` (approximately) of unit length.
The tolerance for the last test can be set using the `kwargs...`.
"""
function check_manifold_point(M::Sphere{N}, p; kwargs...) where {N}
    mpv = invoke(check_manifold_point, Tuple{supertype(typeof(M)), typeof(p)}, M, p; kwargs...)
    mpv === nothing || return mpv
    if !isapprox(norm(p), 1.0; kwargs...)
        return DomainError(
            norm(p),
            "The point $(p) does not lie on the sphere $(M) since its norm is not 1.",
        )
    end
    return nothing
end

"""
    check_tangent_vector(M, p, X; check_base_point = true, kwargs... )

Check whether `X` is a tangent vector to `p` on the [`Sphere`](@ref) `M`, i.e.
after [`check_manifold_point`](@ref)`(M,p)`, `X` has to be of same dimension as `p`
and orthogonal to `p`.
The optional parameter `check_base_point` indicates, whether to call [`check_manifold_point`](@ref)  for `p`.
The tolerance for the last test can be set using the `kwargs...`.
"""
function check_tangent_vector(
    M::Sphere{N},
    p,
    X;
    check_base_point = true,
    kwargs...,
) where {N}
    if check_base_point
        mpe = check_manifold_point(M, p; kwargs...)
        mpe === nothing || return mpe
    end
    mpv = invoke(
        check_tangent_vector,
        Tuple{supertype(typeof(M)), typeof(p), typeof(X)},
        M,
        p,
        X;
        check_base_point = false, # already checked above
        kwargs...
    )
    mpv === nothing || return mpv
    if !isapprox(abs(dot(p, X)), 0.0; kwargs...)
        return DomainError(
            abs(dot(p, X)),
            "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not orthogonal in the embedding.",
        )
    end
    return nothing
end

decorated_manifold(M::Sphere) = Euclidean(representation_size(M)...; field=ℝ)

@doc raw"""
    distance(M::Sphere, p, q)

Compute the geodesic distance betweeen `p` and `q` on the [`Sphere`](@ref) `M`.
The formula is given by the (shorter) great arc length on the (or a) great circle
both `p` and `q` lie on.

````math
d_{𝕊^n}(p,q) = \arccos(⟨p,q⟩).
````
"""
distance(::Sphere, x, y) = acos(clamp(dot(x, y), -1, 1))

embed!(::Sphere, q, p) = (q .= p)

embed!(::Sphere, Y, p, X) = (Y .= X)

@doc raw"""
    exp(M::Sphere, p, X)

Compute the exponential map from `p` in the tangent direction `X` on the [`Sphere`](@ref)
`M` by following the great arc eminating from `p` in direction `X`.

````math
\exp_p X = \cos(\lVert X \rVert_p)p + \sin(\lVert X \rVert_p)\frac{X}{\lVert X \rVert_p}X,
````
where $\lVert X \rVert_p$ is the [`norm`](@ref norm(::Sphere,p,X)) on the
[`Sphere`](@ref) `M`.
"""
exp(::Sphere, ::Any...)

function exp!(M::Sphere, q, p, X)
    θ = norm(M, p, X)
    q .= cos(θ) .* p .+ usinc(θ) .* X
    return q
end

flat!(M::Sphere, ξ::CoTFVector, p, X::TFVector) = copyto!(ξ, X)

function get_basis(M::Sphere{N}, p, B::DiagonalizingOrthonormalBasis) where {N}
    A = zeros(N + 1, N + 1)
    A[1, :] = transpose(p)
    A[2, :] = transpose(B.frame_direction)
    V = nullspace(A)
    κ = ones(N)
    if !iszero(B.frame_direction)
        # if we have a nonzero direction for the geodesic, add it and it gets curvature zero from the tensor
        V = cat(B.frame_direction / norm(M, p, B.frame_direction), V; dims = 2)
        κ[1] = 0 # no curvature along the geodesic direction, if x!=y
    end
    Ξ = [V[:, i] for i = 1:N]
    return CachedBasis(B, κ, Ξ)
end

@doc raw"""
    get_coordinates(M::Sphere, p, X, B::DefaultOrthonormalBasis)

Represent the tangent vector `X` at point `p` from the [`Sphere`](@ref) `M` in
an orthonormal basis by rotating the vector `X` using the rotation matrix
$2\frac{q q^\mathrm{T}}{q^\mathrm{T} q} - I$ where $q = p + (1, 0, …, 0)$.
"""
function get_coordinates(M::Sphere{N}, p, X, B::DefaultOrthonormalBasis) where {N}
    if isapprox(p[1], 1)
        return X[2:end]
    else
        xp1 = p .+ ntuple(i -> ifelse(i == 1, 1, 0), N + 1)
        return (2*xp1*dot(xp1, X)/dot(xp1, xp1)-X)[2:end]
    end
end

function get_coordinates!(M::Sphere, Y, p, X, B::DefaultOrthonormalBasis)
    return copyto!(Y, get_coordinates(M, p, X, B))
end

function get_vector(M::Sphere{N}, p, X, B::DefaultOrthonormalBasis) where {N}
    p[1] ≈ 1 && return vcat(0, X)
    xp1 = p .+ ntuple(i -> ifelse(i == 1, 1, 0), N + 1)
    X0 = vcat(0, X)
    return 2 * xp1 * dot(xp1, X0) / dot(xp1, xp1) - X0
end

function get_vector!(M::Sphere, Y::AbstractVector, p, X, B::DefaultOrthonormalBasis)
    return copyto!(Y, get_vector(M, p, X, B))
end

@doc raw"""
    injectivity_radius(M::Sphere[, p])

Return the injectivity radius for the [`Sphere`](@ref) `M`, which is globally $π$.

    injectivity_radius(M::Sphere, x, ::ProjectionRetraction)

Return the injectivity radius for the [`ProjectionRetraction`](@ref) on the
[`Sphere`](@ref), which is globally $\frac{π}{2}$.
"""
injectivity_radius(::Sphere) = π
injectivity_radius(::Sphere, ::ExponentialRetraction) = π
injectivity_radius(::Sphere, ::ProjectionRetraction) = π / 2
injectivity_radius(::Sphere, ::Any) = π
injectivity_radius(::Sphere, ::Any, ::ExponentialRetraction) = π
injectivity_radius(::Sphere, ::Any, ::ProjectionRetraction) = π / 2
eval(quote
    @invoke_maker 1 Manifold injectivity_radius(M::Sphere, rm::AbstractRetractionMethod)
end)

@doc raw"""
    inverse_retract(M::Sphere, p, q, ::ProjectionInverseRetraction)

Compute the inverse of the projection based retraction on the [`Sphere`](@ref),
i.e. rearranging $p+X = q\lVert p+X\rVert_2$ yields
since $⟨p,X⟩ = 0$ and when $d_{𝕊^2}(p,q) ≤ \frac{π}{2}$ that

````math
\operatorname{retr}_p^{-1}(q) = \frac{q}{⟨p, q⟩} - p.
````
"""
inverse_retract(::Sphere, ::Any, ::Any, ::ProjectionInverseRetraction)

function inverse_retract!(::Sphere, X, p, q, ::ProjectionInverseRetraction)
    return (X .= q ./ dot(p, q) .- p)
end

@doc raw"""
    log(M::Sphere, p, q)

Compute the logarithmic map on the [`Sphere`](@ref) `M`, i.e. the tangent vector,
whose geodesic starting from `p` reaches `q` after time 1.
The formula reads for $x ≠ -y$

````math
\log_p q = d_{𝕊^n}(p,q) \frac{q-⟨p,q⟩ p}{\lVert q-⟨p,q⟩ p \rVert_2},
````

and a deterministic choice from the set of tangent vectors is returned if $x=-y$, i.e. for
opposite points.
"""
log(::Sphere, ::Any...)

function log!(S::Sphere, X, p, q)
    cosθ = dot(p, q)
    if cosθ ≈ -1 # appr. opposing points, return deterministic choice from set-valued log
        fill!(X, 0)
        if p[1] ≈ 1
            X[2] = 1
        else
            X[1] = 1
        end
        copyto!(X, X .- dot(p, X) .* p)
        X .*= π / norm(X)
    else
        cosθ = cosθ > 1 ? one(cosθ) : cosθ
        θ = acos(cosθ)
        X .= (q .- cosθ .* p) ./ usinc(θ)
    end
    return project!(S, X, p, X)
end

@doc raw"""
    manifold_dimension(M::Sphere)

Return the dimension of the [`Sphere`](@ref)`(n) `M`, i.e. $𝕊^n$, which is $\dim(𝕊^n) = n$.
"""
manifold_dimension(S::Sphere{N}) where {N} = N

"""
    mean(
        S::Sphere,
        x::AbstractVector,
        [w::AbstractWeights,]
        method = GeodesicInterpolationWithinRadius(π/2);
        kwargs...,
    )

Compute the Riemannian [`mean`](@ref mean(M::Manifold, args...)) of `x` using
[`GeodesicInterpolationWithinRadius`](@ref).
"""
mean(::Sphere, ::Any...)

function mean!(S::Sphere, p, x::AbstractVector, w::AbstractVector; kwargs...)
    return mean!(S, p, x, w, GeodesicInterpolationWithinRadius(π / 2); kwargs...)
end

"""
    normal_tvector_distribution(S::Sphere, p, σ)

Normal distribution in ambient space with standard deviation `σ`
projected to tangent space at `p`.
"""
function normal_tvector_distribution(S::Sphere, p, σ)
    d = Distributions.MvNormal(zero(p), σ)
    return ProjectedFVectorDistribution(TangentBundleFibers(S), p, d, project_vector!, p)
end

@doc raw"""
    project(M::Sphere, p)

Project the point `p` from the embedding onto the [`Sphere`](@ref) `M`.

````math
\operatorname{proj}_{𝕊^n}(p) = \frac{p}{\lVert p \rVert_2}.
````
"""
project(::Sphere, ::Any)

project!(S::Sphere, q, p) = copyto!(q, p./ norm(p))

@doc raw"""
    project(M::Sphere, p, X)

Project the point `X` onto the tangent space at `p` on the [`Sphere`](@ref) `M`.

````math
\operatorname{proj}_{p}(X) = X - ⟨p, X⟩p
````
"""
project(::Sphere, ::Any, ::Any)

project!(S::Sphere, Y, p, X) = (Y .= X .- dot(p, X) .* p)

@doc raw"""
    representation_size(M::Sphere)

Return the size points on the [`Sphere`](@ref) `M` are represented as, i.e.
for the `n`-dimensional [`Sphere`](@ref) it is vectors of size `(n+1,)`.
"""
@generated representation_size(::Sphere{N}) where {N} = (N + 1,)

@doc raw"""
    retract(M::Sphere, p, X, ::ProjectionRetraction)

Compute the retraction that is based on projection, i.e.

````math
\operatorname{retr}_p(X) = \frac{p+X}{\lVert p+X \rVert_2}
````
"""
retract(::Sphere, ::Any, ::Any, ::ProjectionRetraction)

function retract!(M::Sphere, q, p, X, ::ProjectionRetraction)
    q .= p .+ X
    return project!(M, q, q)
end

show(io::IO, ::Sphere{N}) where {N} = print(io, "Sphere($(N))")

"""
    uniform_distribution(M::Sphere, p)

Uniform distribution on given [`Sphere`](@ref) `M`. Generated points will be of
similar type as `p`.
"""
function uniform_distribution(M::Sphere, p)
    d = Distributions.MvNormal(zero(p), 1.0)
    return ProjectedPointDistribution(M, d, project!, p)
end

@doc doc"""
    vector_transport_to(M, p, X, q, ::ParallelTransport)

Compute the paralllel transport of the tangent vector `X` at `p` to `q`,
provided, the [`geodesic`](@ref) between `p` and `q` is unique

````math
P_{p←q}(X) = X - \frac{\langle \log_p q,X\rangle_p}{d^2_{𝕊^n}(p,q)}
\bigl(\log_xy + \log_yx \bigr).
````
"""
vector_transport_to(::Sphere, ::Any, ::Any, ::Any, ::Any, ::ParallelTransport)

function vector_transport_to!(M::Sphere, Y, p, X, q, ::ParallelTransport)
    X_pq = log(M, p, q)
    Xl = norm(M, p, X_pq)
    copyto!(Y, X)
    if Xl > 0
        factor = 2 * dot(X, q) / (norm(p + q)^2)
        Y .-= factor .* (p .+ q)
    end
    return Y
end

ext: jl


text:
# module for representing observation edge

using LinearAlgebra

mutable struct ObsrvEdge
  t1
  t2
  z1
  z2
  x1
  x2
  omega_upper_right
  omega_upper_left
  omega_lower_right
  omega_lower_left
  xi_upper
  xi_lower

  # init
  function ObsrvEdge(t1, t2, z1, z2, pose_list;
                     sensor_noise_rate=[0.14,0.05])
    @assert z1[1] == z2[1] # prevent id is different

    self = new()

    self.t1 = t1 # time
    self.t2 = t2
    self.z1 = z1[2] # observation [distance, direction orientation]
    self.z2 = z2[2]
    self.x1 = pose_list[string(t1)] # pose [x, y, theta]
    self.x2 = pose_list[string(t2)]

    sin1 = sin(self.x1[3] + self.z1[2])
    cos1 = cos(self.x1[3] + self.z1[2])
    sin2 = sin(self.x2[3] + self.z2[2])
    cos2 = cos(self.x2[3] + self.z2[2])

    # calculate error
    error_pose = self.x2[1:2] - self.x1[1:2]
    error_obsrv = [
      self.z2[1]*cos2-self.z1[1]*cos1,
      self.z2[1]*sin2-self.z1[1]*sin1
    ]
    error = error_pose + error_obsrv

    # precision matrix for edge
    Q1 = diagm(0 => [(self.z1[1]*sensor_noise_rate[1])^2,
                     sensor_noise_rate[2]^2])
    R1 = -[cos1 -self.z1[1]*sin1;
           sin1  self.z1[1]*cos1]
    Q2 = diagm(0 => [(self.z2[1]*sensor_noise_rate[1])^2,
                     sensor_noise_rate[2]^2])
    R2 = [cos2 -self.z2[1]*sin2;
          sin2  self.z2[1]*cos2]
    sigma_edge = R1*Q1*R1' + R2*Q2*R2'
    omega_edge = inv(sigma_edge)

    # precision matrix for graph
    B1 = -[1 0 -self.z1[1]*sin1;
           0 1  self.z1[1]*cos1]
    B2 = [1 0 -self.z2[1]*sin2;
          0 1  self.z2[1]*cos2]
    self.omega_upper_left = B1' * omega_edge * B1
    self.omega_upper_right = B1' * omega_edge * B2
    self.omega_lower_left = B2' * omega_edge * B1
    self.omega_lower_right = B2' * omega_edge * B2
    # coefficient vector for graph
    self.xi_upper = -B1' * omega_edge * error
    self.xi_lower = -B2' * omega_edge * error

    return self
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
"""
`ParameterTable`s contain the specification of a structural equation model.

# Constructor

    (1) ParameterTable(;graph, observed_vars, latent_vars, ...)

    (2) ParameterTable(ram_matrices)

Return a `ParameterTable` constructed from (1) a graph or (2) RAM matrices.

# Arguments
- `graph`: graph defined via `@StenoGraph`
- `observed_vars::Vector{Symbol}`: observed variable names
- `latent_vars::Vector{Symbol}`: latent variable names
- `ram_matrices::RAMMatrices`: a `RAMMatrices` object

# Examples
See the online documentation on [Model specification](@ref) and the [ParameterTable interface](@ref).

# Extended help
## Additional keyword arguments
- `parname::Symbol = :θ`: prefix for automatically generated parameter labels
"""
function ParameterTable end

"""
`EnsembleParameterTable`s contain the specification of an ensemble structural equation model.

# Constructor

    (1) EnsembleParameterTable(;graph, observed_vars, latent_vars, groups)

    (2) EnsembleParameterTable(args...; groups)

Return an `EnsembleParameterTable` constructed from (1) a graph or (2) multiple RAM matrices.

# Arguments
- `graph`: graph defined via `@StenoGraph`
- `observed_vars::Vector{Symbol}`: observed variable names
- `latent_vars::Vector{Symbol}`: latent variable names
- `groups::Vector{Symbol}`: group names
- `args...`: `RAMMatrices` for each model

# Examples
See the online documentation on [Multigroup models](@ref).
"""
function EnsembleParameterTable end

"""
`RAMMatrices` contain the specification of a structural equation model.

# Constructor

    (1) RAMMatrices(partable::ParameterTable)

    (2) RAMMatrices(;A, S, F, M = nothing, parameters, colnames)

    (3) RAMMatrices(partable::EnsembleParameterTable)

Return `RAMMatrices` constructed from (1) a parameter table or (2) individual matrices.

(3) Return a dictionary of `RAMMatrices` from an `EnsembleParameterTable` (keys are the group names).

# Arguments
- `partable`
- `A`: matrix of directed effects
- `S`: matrix of undirected effects
- `F`: filter matrix
- `M`: vector of mean effects
- `parameters::Vector{Symbol}`: parameter labels
- `colnames::Vector{Symbol}`: variable names corresponding to the A, S and F matrix columns

# Examples
See the online documentation on [Model specification](@ref) and the [RAMMatrices interface](@ref).
"""
function RAMMatrices end

"""
    fixed(args...)

Fix parameters to a certain value.
For ensemble models, multiple values (one for each submodel/group) are needed.

# Examples
```julia
graph = @StenoGraph begin
    x → fixed(1)*y
end
```
"""
function fixed end

"""
    start(args...)

Define starting values for parameters.
For ensemble models, multiple values (one for each submodel/group) are needed.

# Examples
```julia
graph = @StenoGraph begin
    x → start(1)*y
end
```
"""
function start end

"""
    label(args...)

Label parameters.
For ensemble models, multiple values (one for each submodel/group) are needed.

# Examples
```julia
graph = @StenoGraph begin
    x → label(:a)*y
end
```
"""
function label end


text:
using WinKnownPaths
using UUIDs
using Test

import WinKnownPaths: FOLDERID
@testset "WinKnownPaths.jl" begin
    if Sys.iswindows()
        @test !isempty(WinKnownPaths.path(FOLDERID.System))
        @test !isempty(WinKnownPaths.path(FOLDERID.LocalAppData))
        @test_throws SystemError WinKnownPaths.path(UUID(UInt128(0)))
    end
end


text:
# Contains functions for generating synthetic training data

"""
    generate_training_data(unknown_fn)

Generates synthetic data from an unknown function by uniformly sampling datapoints over the domain. Returns the input and output points.

    Note that you cannot use the `measurement_dist` and `known_fn` together.

# Arguments
- `unknown_fn`: the handle of the unknown function.
- `domain::Dict`: the dictionary indicating the state space domain with keys `"x1","x2"...` and corresponding min-max bounds.
- `data_num::Int`: the number of datapoints to generate.
- `known_fn=nothing`: the known part of the function.
- `random_seed::Int=11`: random seed for the uniform sample (and noise processes, if any).
- `n_dims_out::Int=-1`: dimension of the output of the function, defaults to same dimension as input.
- `process_dist=nothing`: process noise distristribution.
- `measurement_dist=nothing`: measurement noise distristribution.
"""
function generate_training_data(unknown_fn, domain::Dict, data_num::Int;
                                known_fn=nothing, random_seed::Int=11, n_dims_out::Int=-1, process_dist=nothing, measurement_dist=nothing)

    if !isnothing(known_fn)
        f_sub = (x) -> unknown_fn(x) + known_fn(x)
    else
        f_sub = (x) -> unknown_fn(x)
    end

    n_dims_in = length(domain)
    n_dims_out = n_dims_out > 0 ? n_dims_out : n_dims_in
    mt = MersenneTwister(random_seed)

    x_train = vcat([rand(mt, Uniform(domain["x$i"][1], domain["x$i"][2]), 1, data_num) for i=1:n_dims_in]...)
    y_train = mapslices(f_sub, x_train, dims=1)

    # Account for process noise...
    if !isnothing(process_dist)
        y_train += rand(mt, process_dist, (n_dims_out, data_num))
    end

    # Account for measurement noise...
    if !isnothing(measurement_dist)
        y_train += rand(mt, measurement_dist, (n_dims_out, data_num))
    end

    # Account for the known part, if any...
    if !isnothing(known_fn)
        y_train = y_train - mapslices(known_fn, x_train, dims=1)
    end

    @assert size(x_train, 2) == size(y_train, 2)

    return x_train, y_train
end

function save_training_data(x_train, y_train, exp_dir::String)
    training_data = Dict()
    training_data[:input] = x_train
    training_data[:output] = y_train
    bson("$exp_dir/training_data.bson", training_data)
end

function load_training_data(filename::String)
    training_data = BSON.load(filename)
    return training_data[:input], training_data[:output]
end


text:
using Convex, MultiConvex, Compat, SCS
set_default_solver(SCSSolver(verbose=false))
using Base.Test
import Convex.hash

x = Variable()
y = Variable()
z = Variable()
x.value = 1
y.value = 2
z.value = 3

# variables in
@assert length(variablesin(x*y*z)) == 3
@assert length(variablesin(x*y+y*z+x+y+z)) == 3
@assert length(variablesin(log(x) + 7*x + 5)) == 1

# partitioning sets of integers
n = Dict{Int, Set{Int}}({1=>Set{Int}({2,3,4}),
	                             2=>Set{Int}({1,3}),
	                             3=>Set{Int}({1,2}),
	                             4=>Set{Int}({1}),
	                             5=>Set{Int}({6,7}),
	                             6=>Set{Int}({5,7}),
	                             7=>Set{Int}({6,5})})
p = stablepartition!(n)
@assert length(p) == 3
@assert sum(map(length, p)) == 7
@assert length(union(p...)) == 7

# partitioning given initial sets
p = stablepartition!(n, Set{Int}[Set{Int}({1,5})])
@assert length(p) == 3
@assert sum(map(length, p)) == 7
@assert length(union(p...)) == 7
@assert p[1] == Set({1,5})

# multivexity: (convex) addition
conflicts = Set{@compat Tuple{Uint64, Uint64}}()
@assert multivexity(x+y, conflicts) == AffineVexity()
@assert length(conflicts) == 0

# multivexity: (biconvex) multiplication
conflicts = Set{@compat Tuple{Uint64, Uint64}}()
v = multivexity(x*y, conflicts)
@assert v == AffineVexity()
@assert length(conflicts) == 1
@assert pop!(conflicts) == (x.id_hash, y.id_hash)
# note hash(x::Variable) is not defined to grab id_hash except *internally* in the multiconvex module

# multivexity: addition and multiplication
conflicts = Set{@compat Tuple{Uint64, Uint64}}()
v = multivexity(x*y + y*z + z*x, conflicts)
@assert v == AffineVexity()
@assert length(conflicts) == 3

# partitioning variables given conflict graphs
conflicts = Set{@compat Tuple{Uint64, Uint64}}()
m = x*y + y*z + z*x
v = multivexity(m, conflicts)
vars = variablesin(m)
parts = stablepartition!(conflicts)
@assert length(parts) == 3
@assert all(map(x -> length(x) == 1, parts))

# multivexity: problems with no constraints
x = Variable()
y = Variable()
z = Variable()
problem = minimize(exp(x*y*z) + x^2 + y^2 + z^2)
vex, stablesetvars, vars = multivexity(problem)
@assert length(vars) == 3
@assert length(stablesetvars) == 3
@assert sum(map(length, stablesetvars)) == 3
@assert vex == ConvexVexity()

# multivexity: problems with constraints
x = Variable()
y = Variable()
z = Variable()
w = Variable()
problem = minimize(exp(x*z) + exp(w*y) + x^2 + y^2 + z^2, x+y >= 0)
vex, stablesetvars, vars = multivexity(problem)
@assert length(vars) == 4
@assert length(stablesetvars) == 2
@assert map(length, stablesetvars) == [2, 2]
@assert vex == ConvexVexity()

# solving a simple problem
x = Variable()
y = Variable()
z = Variable()
# solution to this problem should be (0,0,0)
problem = minimize(exp(x*y*z) + x^2 + y^2 + z^2, x>=0, y>=0, z>=0)
altmin!(problem)
@show t1 = (x.value, y.value, z.value)
# solution to this problem should be closer to (some permutation of) (0,1,1)
problem = minimize(exp(x*y*z) + 5*(1-x)^2 + 5*(1-y)^2 + 5*(1-z)^2, x>=0, y>=0, z>=0)
altmin!(problem, warmstart = false, maxAMiters=100)
@show t2 = (x.value[1], y.value[1], z.value[1])
@assert sort([t2...])[2] > maximum(t1)


text:
println("Test OrderedNormalMeanPrior")

@testset "OrderedNormalMeanPrior" begin
  @testset "GMM (ordered prior)" begin
    Random.seed!(0)
    K = 3
    mm = MixtureModel([Normal(-2, .3), Normal(0, .2), Normal(2, .5)], [.3, .4, .3])
    y = rand(mm, 500)

    onm_prior = OrderedNormalMeanPrior(K, Normal(0, 6), truncated(Normal(1, .3), 0, Inf))
    model = MixSkewT(y, K, mu=onm_prior, tdist=false, skew=false)
    print_model_info(model)
    spl = make_sampler(model)

    nsamps = 1000
    chain, metrics = mcmc(spl, nsamps, nburn=2000, thin=1)
    @test length(unique(chain)) == nsamps
    @test isapprox(mean.(mm.components), mean(getindex.(chain, :mu)), atol=0.2)
    @test isapprox(std.(mm.components), mean(getindex.(chain, :sigma)), atol=0.1)
  end

  @testset "MixSkewT (ordered prior)" begin
    Random.seed!(0)
    K = 3
    mm = MixtureModel([SkewT(-5, .2, 10, -7),
                       SkewT(0, .3, 10, -7),
                       SkewT(5, .2, 10, -7)], [.3, .4, .3])
    y = rand(mm, 1000)

    # NOTE:
    # - Prior for psi needs to be loose (e.g. Normal(0, 3). Normal(-2, 1) is
    #   too strict and skew may be extremely negative in posterior.)
    # - OrderedNormalMeanPrior needs to somewhat loose: if too strict, then
    #   difficult to escape bad local modes.
    onm_prior = OrderedNormalMeanPrior(K, Normal(-5, 1), truncated(Normal(1, 3), 0, Inf))
    model = MixSkewT(y, K, mu=onm_prior, psi=Normal(0, 3), eta=Dirichlet(K, 1))
    print_model_info(model)

    # Find best seed.
    seed, min_loglike = find_good_seed(model, nsamps=100, nburn=200, thin=2, seeds=41:60)
    println("Best seed: $seed ($min_loglike)")

    # Define callback.
    nsamps, nburn, thin = (2000, 3000, 2)
    callback = make_callback(model, nsamps=nsamps, nburn=nburn, thin=thin)

    # Run official model.
    Random.seed!(seed)
    spl = make_sampler(model)
    chain, metrics = mcmc(spl, nsamps, nburn=nburn, thin=thin, callback=callback)

    # Tests.
    @test length(unique(chain)) == nsamps
    @test isapprox(getfield.(mm.components, :loc), mean(getindex.(chain, :mu)), atol=0.1)
    @test isapprox(getfield.(mm.components, :scale), mean(getindex.(chain, :sigma)), atol=0.1)
    @test let
      # Test that the 99% CI contains truth.
      nu_post = hcat(getindex.(chain, :nu)...)
      nu_true = getfield.(mm.components, :df)
      # nu_lower = vec(MCMC.quantiles(nu_post, .005, dims=2))
      # nu_upper = vec(MCMC.quantiles(nu_post, .995, dims=2))
      nu_sd = vec(std(nu_post, dims=2))
      nu_mean = vec(mean(nu_post, dims=2))
      nu_lower = nu_mean - 3 * nu_sd
      nu_upper = nu_mean + 3 * nu_sd
      all(nu_lower .< nu_true .< nu_upper)
    end
    @test isapprox(getfield.(mm.components, :skew), mean(getindex.(chain, :phi)), atol=4)
    @test isapprox(mm.prior.p, mean(getindex.(chain, :eta)), atol=.1)

    # TODO: check that the 95% CI contains truth.
  end
end


text:
export rule

@rule GCV(:y, Marginalisation) (m_x::UniNormalOrExpLinQuad, q_z::Any, q_κ::Any, q_ω::Any) = begin

    x_mean, x_var = mean_var(m_x)
    z_mean, z_var = mean_var(q_z)
    κ_mean, κ_var = mean_var(q_κ)
    ω_mean, ω_var = mean_var(q_ω)

    ksi = κ_mean ^ 2 * z_var + z_mean ^ 2 * κ_var + z_var * κ_var
    A = exp(-ω_mean + ω_var / 2)
    B = exp(-κ_mean * z_mean + ksi / 2)

    return NormalMeanVariance(x_mean, x_var + inv(A * B))
end

@rule GCV(:y, Marginalisation) (q_x::Any, q_z::Any, q_κ::Any, q_ω::Any) = begin

    x_mean        = mean(q_x)
    z_mean, z_var = mean_var(q_z)
    κ_mean, κ_var = mean_var(q_κ)
    ω_mean, ω_var = mean_var(q_ω)

    ksi = κ_mean ^ 2 * z_var + z_mean ^ 2 * κ_var + z_var * κ_var
    A = exp(-ω_mean + ω_var / 2)
    B = exp(-κ_mean * z_mean + ksi / 2)

    return NormalMeanVariance(x_mean, inv(A * B))
end


text:
# This fitting function assumes
function fit_main(Y, X, α_init, model;
                  exposure = nothing,
                  penalty = true, pen_α = 5.0, pen_params = nothing,
                  ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
                  grad_jump = true, grad_seq = nothing,
                  print_steps = true)

    # Make variables accessible within the scope of `let`
    let α_em, gate_em, model_em, model_em_expo, ll_em_list, ll_em, ll_em_np, ll_em_old, ll_em_np_old, iter, z_e_obs, z_e_lat, k_e, params_old
        # Initial loglik
        gate_init = LogitGating(α_init, X)
        model_expo = exposurize_model(model, exposure = exposure) # exposurize
        ll_np_list = loglik_np(Y, gate_init, model_expo)
        ll_init_np = ll_np_list.ll
        ll_penalty = penalty ? (penalty_α(α_init, pen_α) + penalty_params(model, pen_params)) : 0.0
        ll_init = ll_init_np + ll_penalty

        if print_steps
            @info("Initial loglik: $(ll_init_np) (no penalty), $(ll_init) (with penalty)")
        end

        # start em
        α_em = copy(α_init)
        model_em = copy(model)
        model_em_expo = exposurize_model(model_em, exposure = exposure)
        gate_em = LogitGating(α_em, X)
        ll_em_list = loglik_np(Y, gate_em, model_em_expo)
        ll_em_np = ll_em_list.ll
        ll_em = ll_init
        ll_em_old = -Inf
        iter = 0

        while (ll_em - ll_em_old > ϵ) & (iter < ecm_iter_max)
            # Update counter and loglikelihood
            iter = iter + 1
            ll_em_np_old = ll_em_np
            ll_em_old = ll_em

            # E-Step
            z_e_obs = EM_E_z_obs(ll_em_list.gate_expert_ll_comp, ll_em_list.gate_expert_ll)
            z_e_lat = EM_E_z_lat(ll_em_list.gate_expert_tn_bar_comp_z_lat, ll_em_list.gate_expert_tn_bar_z_lat)
            # k_e = EM_E_k(ll_em_list.gate_expert_tn)
            k_e = EM_E_k(ll_em_list.gate_expert_tn_bar_k)

            # M-Step: α
            ll_em_temp = ll_em
            α_em = EM_M_α(X, α_em, z_e_obs, z_e_lat, k_e, α_iter_max = α_iter_max, penalty = penalty, pen_α = pen_α)
            gate_em = LogitGating(α_em, X)
            ll_em_list = loglik_np(Y, gate_em, model_em_expo)
            ll_em_np = ll_em_list.ll
            ll_em_penalty = penalty ? (penalty_α(α_em, pen_α) + penalty_params(model_em, pen_params)) : 0.0
            ll_em = ll_em_np + ll_em_penalty

            s = ll_em - ll_em_temp > 0 ? "+" : "-"
            pct = abs((ll_em - ll_em_temp) / ll_em_temp) * 100
            if print_steps
                @info("Iteration $(iter), updating α: $(ll_em_old) ->  $(ll_em), ( $(s) $(pct) % )")
            end
            ll_em_temp = ll_em

            # M-Step: component distributions
            for j in 1:size(model)[1] # by dimension
            # for j in 1:1
                for k in 1:size(model)[2] # by component

                    params_old = params(model_em[j,k])

                    model_em[j,k] = EM_M_expert(model_em[j,k],
                                                Y[:, 4*(j-1)+1], Y[:, 4*(j-1)+2], Y[:, 4*(j-1)+3],Y[:, 4*(j-1)+4],
                                                exposure,
                                                # ll_em_list.expert_ll_pos_dim_comp[j][:,k],
                                                # ll_em_list.expert_tn_pos_dim_comp[j][:,k],
                                                # ll_em_list.expert_tn_bar_pos_dim_comp[j][:,k],
                                                vec(z_e_obs[:,k]), vec(z_e_lat[:,k]), vec(k_e),
                                                penalty = penalty, pen_pararms_jk = pen_params[j][k])

                    model_em_expo = exposurize_model(model_em, exposure = exposure)
                    ll_em_list = loglik_np(Y, gate_em, model_em_expo)
                    ll_em_np = ll_em_list.ll
                    ll_em_penalty = penalty ? (penalty_α(α_em, pen_α) + penalty_params(model_em, pen_params)) : 0.0
                    ll_em = ll_em_np + ll_em_penalty

                    if print_steps
                        s = ll_em - ll_em_temp > 0 ? "+" : "-"
                        pct = abs((ll_em - ll_em_temp) / ll_em_temp) * 100
                        @info("Iteration $(iter), updating model[$j, $k]: $(ll_em_temp) ->  $(ll_em), ( $(s) $(pct) % )")
                        if s=="-" @info("Intended update of params: $(params_old) ->  $(params(model_em[j,k]))") end
                    end
                    ll_em_temp = ll_em
                end
            end

            α_em = α_em
            model_em = model_em
            model_em_expo = exposurize_model(model_em, exposure = exposure)
            gate_em = LogitGating(α_em, X)
            ll_em_list = loglik_np(Y, gate_em, model_em_expo)
            ll_em_np = ll_em_list.ll
            ll_em_penalty = penalty ? (penalty_α(α_em, pen_α) + penalty_params(model_em, pen_params)) : 0.0
            ll_em = ll_em_np + ll_em_penalty

        end

        converge = (ll_em - ll_em_old > ϵ) ? false : true
        AIC = -2.0*ll_em_np + 2*(_count_α(α_em) + _count_params(model_em))
        BIC = -2.0*ll_em_np + log(size(Y)[1])*(_count_α(α_em) + _count_params(model_em))

        return (α_fit = α_em, model_fit = model_em,
                converge = converge, iter = iter,
                ll_np = ll_em_np, ll = ll_em,
                AIC = AIC, BIC = BIC)
    end

end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
export ValueTrans

"""
Values obtained from the `inner` family are transformed by
`trans_function(x)`, which should be invertible when used with `fit`.

Basis is not defined, since the transformation is generally not
linear.
"""
immutable ValueTrans{T, S} <: ParametricFamily
    trans_function::T
    inner::S
end

@forward ValueTrans.inner degf, domain, points

function basis!(::ValueTrans, ::Any, ::Any)
    throw(ArgumentError("Value transformations are generally not linear, basis not defined."))
end

evaluate(p::ValueTrans, θ, x) = p.trans_function(x)(evaluate(p.inner, θ, x))

function fit!(p::ValueTrans, ys::AbstractVector, θ)
    zs = [inv(p.trans_function(x))(y) for (x,y) in zip(points(p.inner), ys)]
    fit!(p.inner, zs, θ)
end


text:
## Overlay for animation
@with_kw struct BeliefOverlay <: SceneOverlay
    # belief::SingleOCFBelief = SingleOCFBelief()
     b_dict::Dict{Int64, SingleOCFBelief} = Dict{Int64, SingleOCFBelief}()
     ego_vehicle::Vehicle = Vehicle(VehicleState(VecSE2(0., 0., 0.), 0.), VehicleDef(), 1)
     color::Colorant = RGBA(0.0, 0.0, 1.0, 0.2)
 end

 function AutoViz.render!(rendermodel::RenderModel, overlay::BeliefOverlay, scene::Scene, roadway::R) where R

     b_dict = overlay.b_dict
     ego = overlay.ego_vehicle

     for b_id in keys(b_dict)
         for (s, prob) in weighted_iterator( b_dict[b_id])

             ped = Vehicle(VehicleState(VecSE2(s.ped_s+ego.def.length/2+AutomotivePOMDPs.PEDESTRIAN_DEF.width/2, s.ped_T, s.ped_theta), 0.), AutomotivePOMDPs.PEDESTRIAN_DEF, 1)
             if prob > 1e-6
                 prob_color = RGBA(0.0, 0.0, 1.0, prob*10)
                 add_instruction!(rendermodel, render_vehicle, (ego.state.posG.x+ped.state.posG.x, ego.state.posG.y+ped.state.posG.y, ego.state.posG.θ + ped.state.posG.θ, ped.def.length, ped.def.width, prob_color))
             end
         end
     end

     return rendermodel
 end

function AutomotivePOMDPs.animate_record(scenes::Vector{Scene},dt::Float64, env::CrosswalkEnv, sensor::GaussianSensor, sensor_o::Vector{Vector{Vehicle}}, risk::Vector{Float64}, belief_dict::Vector{Dict{Int64, SingleOCFBelief}}, ego_vehicle::Vector{Vehicle}, action_pomdp::Vector{SingleOCFAction}, prediction::Vector{Array{Float64}}, cam=FitToContentCamera(0.0))

    duration =length(scenes)*dt::Float64
    fps = Int(1/dt)
    function render_rec(t, dt)

        frame_index = Int(floor(t/dt)) + 1
        text_to_visualize = [string("v-ego: ", round(ego_vehicle[frame_index].state.v*3.6,digits=1), " km/h" ,
                                " y: ", ego_vehicle[frame_index].state.posG.y, " m",
                                " ax: ", action_pomdp[frame_index].acc, " m/s²")]
        sensor_overlay = GaussianSensorOverlay(sensor=sensor, o=sensor_o[frame_index])
        occlusion_overlay = OcclusionOverlay(obstacles=env.obstacles)

        text_overlay = TextOverlay(text=text_to_visualize,pos=VecE2(80.,10.),incameraframe=true,color=colorant"white",font_size=20)
        belief_overlay = BeliefOverlay(b_dict=belief_dict[frame_index], ego_vehicle=ego_vehicle[frame_index])
      #  max_speed = 14.0
      #  histogram_overlay = HistogramOverlay(pos = VecE2(65.0, 10.0), val=ego_vehicle[frame_index].state.v/max_speed, label="v speed")

        if length(prediction) > 0
            prediction_overlay = EmergencyBrakingSystem.PredictionOverlay(prediction=prediction[frame_index])
            return AutoViz.render(scenes[frame_index], env, [prediction_overlay, belief_overlay, sensor_overlay, occlusion_overlay, text_overlay], cam=cam)
        else
            return AutoViz.render(scenes[frame_index], env, [belief_overlay, sensor_overlay, occlusion_overlay, text_overlay], cam=cam)
        end
    end
    return duration, fps, render_rec
end


text:
isblock(x) = length(x) == 1 && x[1] isa Expr && x[1].head == :block
function flatten_expr!(x)
    isb = isblock(x)
    if isb
        x = MacroTools.striplines(x[1])
        filter!(z->z isa Symbol || z.head != :line, x.args)
        x = (x.args...,)
    end
    x
end
function build_expr(head::Symbol, args)
    ex = Expr(head)
    append!(ex.args, args)
    ex
end

# Sym / Term --> Symbol
Base.Symbol(x::Union{Num,Symbolic}) = tosymbol(x)
tosymbol(x; kwargs...) = x
tosymbol(x::Sym; kwargs...) = nameof(x)
tosymbol(t::Num; kwargs...) = tosymbol(value(t); kwargs...)

"""
    diff2term(x::Term) -> Term
    diff2term(x) -> x

Convert a differential variable to a `Term`. Note that it only takes a `Term`
not a `Num`.

```julia
julia> ModelingToolkit.diff2term(ModelingToolkit.value(D(D(x))))
xˍtt(t)
```
"""
function diff2term(O)
    istree(O) || return O
    if is_derivative(O)
        (x, t, order) = flatten_differential(O)
        return lower_varname(x, t, order)
    end
    return Term{Real}(operation(O), diff2term.(arguments(O)))
end

function flatten_differential(O::Term)
    @assert is_derivative(O) "invalid differential: $O"
    is_derivative(arguments(O)[1]) || return (arguments(O)[1], operation(O).x, 1)
    (x, t, order) = flatten_differential(arguments(O)[1])
    isequal(t, operation(O).x) || throw(ArgumentError("non-matching differentials on lhs: $t, $(operation(O).x)"))
    return (x, t, order + 1)
end

"""
    tosymbol(x::Union{Num,Symbolic}; states=nothing, escape=true) -> Symbol

Convert `x` to a symbol. `states` are the states of a system, and `escape`
means if the target has escapes like `val"y⦗t⦘"`. If `escape` then it will only
output `y` instead of `y⦗t⦘`.

# Examples

```julia
julia> @parameters t; @variables z(t)
(z(t),)

julia> ModelingToolkit.tosymbol(z)
Symbol("z⦗t⦘")
```
"""
function tosymbol(t::Term; states=nothing, escape=true)
    if operation(t) isa Sym
        if states !== nothing && !(any(isequal(t), states))
            return nameof(operation(t))
        end
        op = nameof(operation(t))
        args = arguments(t)
    elseif operation(t) isa Differential
        term = diff2term(t)
        op = Symbol(operation(term))
        args = arguments(term)
    else
        @goto err
    end

    return escape ? Symbol(op, "⦗", join(args, ", "), "⦘") : op
    @label err
    error("Cannot convert $t to a symbol")
end

function lower_varname(var::Symbolic, idv, order)
    order == 0 && return var
    name = string(nameof(operation(var)))
    underscore = 'ˍ'
    idx = findlast(underscore, name)
    append = string(idv)^order
    if idx === nothing
        newname = Symbol(name, underscore, append)
    else
        nidx = nextind(name, idx)
        newname = Symbol(name[1:idx], name[nidx:end], append)
    end
    return Sym{symtype(operation(var))}(newname)(arguments(var)[1])
end

function lower_varname(t::Symbolic, iv)
    var, order = var_from_nested_derivative(t)
    lower_varname(var, iv, order)
end
lower_varname(t::Sym, iv) = t


text:
# Autogenerated wrapper script for Glib_jll for armv6l-linux-musleabihf
export libgio, libglib, libgmodule, libgobject, libgthread

using Libiconv_jll
using Libffi_jll
using Gettext_jll
using PCRE_jll
using Zlib_jll
using Libmount_jll
JLLWrappers.@generate_wrapper_header("Glib")
JLLWrappers.@declare_library_product(libgio, "libgio-2.0.so.0")
JLLWrappers.@declare_library_product(libglib, "libglib-2.0.so.0")
JLLWrappers.@declare_library_product(libgmodule, "libgmodule-2.0.so.0")
JLLWrappers.@declare_library_product(libgobject, "libgobject-2.0.so.0")
JLLWrappers.@declare_library_product(libgthread, "libgthread-2.0.so.0")
function __init__()
    JLLWrappers.@generate_init_header(Libiconv_jll, Libffi_jll, Gettext_jll, PCRE_jll, Zlib_jll, Libmount_jll)
    JLLWrappers.@init_library_product(
        libgio,
        "lib/libgio-2.0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libglib,
        "lib/libglib-2.0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgmodule,
        "lib/libgmodule-2.0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgobject,
        "lib/libgobject-2.0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgthread,
        "lib/libgthread-2.0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "nash_conv" begin
    env = OpenSpielEnv("kuhn_poker")
    p = TabularRandomPolicy()

    @test RLZoo.nash_conv(p, env) == 11 / 12

    p = get_optimal_kuhn_policy()
    @test RLZoo.nash_conv(p, env) == 0.0

    env = OpenSpielEnv("leduc_poker")
    p = TabularRandomPolicy()
    @test isapprox(RLZoo.nash_conv(p, env), 4.747222222222222; atol = 0.0001)

    env = OpenSpielEnv("kuhn_poker(players=3)")
    p = TabularRandomPolicy()
    @test RLZoo.nash_conv(p, env) ≈ 2.0624999999999996

    env = OpenSpielEnv("kuhn_poker(players=4)")
    p = TabularRandomPolicy()
    @test RLZoo.nash_conv(p, env) ≈ 3.4760416666666663

    env = KuhnPokerEnv()
    p = TabularRandomPolicy()
    @test RLZoo.nash_conv(p, env) == 11 / 12

    p = get_optimal_kuhn_policy(env)
    @test RLZoo.nash_conv(p, env) == 0.0
end


text:
# ============================================================
# Balance Functions
# ============================================================

function bc_balance!(
    ctr::T,
    ctr0::T,
    ctr1::T,
) where {T<:Union{ControlVolume1D,ControlVolume2D,ControlVolumeUS}}
    ctr.w = 2.0 * ctr0.w - ctr1.w
    @. ctr.prim = 2.0 * ctr0.prim - ctr1.prim

    return nothing
end

function bc_balance!(
    ctr::T,
    ctr0::T,
    ctr1::T,
) where {T<:Union{ControlVolume1D1F,ControlVolume2D1F,ControlVolumeUS1F}}
    @. ctr.w = 2.0 * ctr0.w - ctr1.w
    @. ctr.prim = 2.0 * ctr0.prim - ctr1.prim
    @. ctr.f = 2.0 * ctr0.f - ctr1.f

    return nothing
end

function bc_balance!(
    ctr::T,
    ctr0::T,
    ctr1::T,
) where {T<:Union{ControlVolume1D2F,ControlVolume2D2F,ControlVolumeUS2F}}
    @. ctr.w = 2.0 * ctr0.w - ctr1.w
    @. ctr.prim = 2.0 * ctr0.prim - ctr1.prim
    @. ctr.h = 2.0 * ctr0.h - ctr1.h
    @. ctr.b = 2.0 * ctr0.b - ctr1.b

    return nothing
end

function bc_balance!(
    ctr::T,
    ctr0::T,
    ctr1::T,
) where {T<:Union{ControlVolume1D3F,ControlVolume2D3F}}
    @. ctr.w = 2.0 * ctr0.w - ctr1.w
    @. ctr.prim = 2.0 * ctr0.prim - ctr1.prim
    @. ctr.h0 = 2.0 * ctr0.h0 - ctr1.h0
    @. ctr.h1 = 2.0 * ctr0.h1 - ctr1.h1
    @. ctr.h2 = 2.0 * ctr0.h2 - ctr1.h2
    @. ctr.E = 2.0 * ctr0.E - ctr1.E
    @. ctr.B = 2.0 * ctr0.B - ctr1.B
    ctr.ϕ = 2.0 * ctr0.ϕ - ctr1.ϕ
    ctr.ψ = 2.0 * ctr0.ψ - ctr1.ψ
    @. ctr.lorenz = 2.0 * ctr0.lorenz - ctr1.lorenz

    return nothing
end

function bc_balance!(ctr::T, ctr0::T, ctr1::T) where {T<:ControlVolume1D4F}
    @. ctr.w = 2.0 * ctr0.w - ctr1.w
    @. ctr.prim = 2.0 * ctr0.prim - ctr1.prim
    @. ctr.h0 = 2.0 * ctr0.h0 - ctr1.h0
    @. ctr.h1 = 2.0 * ctr0.h1 - ctr1.h1
    @. ctr.h2 = 2.0 * ctr0.h2 - ctr1.h2
    @. ctr.E = 2.0 * ctr0.E - ctr1.E
    @. ctr.B = 2.0 * ctr0.B - ctr1.B
    ctr.ϕ = 2.0 * ctr0.ϕ - ctr1.ϕ
    ctr.ψ = 2.0 * ctr0.ψ - ctr1.ψ
    @. ctr.lorenz = 2.0 * ctr0.lorenz - ctr1.lorenz

    return nothing
end


text:
module Grammar
export tuplejoin, isTerminalIndex, getTerminalValue, getProductionRuleIndexAndProb, getNodeIndexFromSymbol, getSymbolsOfProductionRule, add, subtract, mul, divide

"""
Grammar:
NODE                                := DEFINITION
------------------------------------------------------------------------------------------------
(node_index 1)  expr                := line \n expr | line
(node_index 2)  line                := exo_line | endo_line
(node_index 3)  exo_line            := bool_var '~' dist | float_var '~' float_dist
(node_index 4)  endo_line           := unary_endo_expr | binary_endo_expr | ternary_endo_expr
*(node_index 5)  unary_endo_expr    := bool_var = '!' bool_var | bool_var = bool_var | int_var = int_var | float_var = float_var
*(node_index 6)  binary_endo_expr   := bool_var = bool_binary_op bool_var bool_var_list | int_var = num_binary_op int_var int_var_list | float_var = num_binary_op float_var float_var_list
*(node_index 7)  ternary_endo_expr  := bool_var = bool_var '?' bool_var ':' bool_var | int_var = bool_var '?' int_var ':' int_var | float_var = bool_var '?' float_var ':' float_var

(node_index 8) bool_var_list        := bool_var ',' bool_var_list | bool_var
(node_index 9) int_var_list         := int_var ',' int_var_list '|' int_var
(node_index 10) float_var_list      := float_var ',' float_var_list '|' float_var

(node_index 11) var                 := bool_var | int_var | float_var
(node_index 12) bool_dist           := Bernoulli( bernoulli_params )

/** BELOW ARE TERMINAL NODES */
(node_index 13)  bool_binary_op     := '&' | '|'
(node_index 14)  num_binary_op      := '+' | '-' | '*' | '\'

(node_index 15) bool_var            := 'bool_var_1' | 'bool_var_2' | ... | 'bool_var_5'
(node_index 16) int_var             := 'int_var_1' | 'int_var_2' | ... | 'int_var_5'
(node_index 17) float_var           := 'float_var_1' | 'float_var_2' | ... | 'float_var_5'

(node_index 18) bernoulli_params    := '0.0' | '0.1' | '0.2' | '0.3' | ... | '1.0'
(node_index 19) normal_params       := /* finite number of (mean, variance) pairs */
(node_index 20) uniform_params      := /* finite number of [start, end] intervals */

(node_index 21) float_dist          := Normal ( normal_params ) | Uniform( uniform_params )

"""

""" Binary Operators """
function add(args...)
    total = 0
    for arg in args
        total += arg
    end
    total
end

function subtract(args...)
    total = args[1]
    for (i,) in enumerate(args)
        if i != 1
            total -= args[i]
        end
    end
    total
end

function mul(args...)
    total = 1
    for arg in args
        total *= arg
    end
    total
end

function divide(args...)
    total = args[1]
    for (i,) in enumerate(args)
        if i != 1
            if args[i] != 0
                total /= args[i]
            else
                if total >= 0
                    Inf
                else
                    -Inf
                end
            end
        end
    end
    total
end

""" ----- GRAMMAR-RELATED UTILITY FUNCTIONS ----- """

""" Concatenates tuples """
tuplejoin(t1::Tuple, t2::Tuple, t3...) = tuplejoin((t1..., t2...), t3...)
tuplejoin(t::Tuple) = t

""" Returns whether node_index refers to a terminal node """
function isTerminalIndex(node_index)
    if (node_index >= 13 && node_index <= 20)
        true
    elseif ((node_index >= 1 && node_index < 13) || node_index == 21)
        false
    else
        throw(ArgumentError(""))
    end
end

""" Returns value associated with terminal_node_index and production_rule_index """
function getTerminalValue(terminal_node_index, production_rule_index)
    if isTerminalIndex(terminal_node_index)
        if (terminal_node_index == 13)
            values = [:and, :or]
        elseif (terminal_node_index == 14)
            values = [:add, :subtract, :mul, :divide]
        elseif (terminal_node_index == 15)
            values = [:bool_var_1, :bool_var_2, :bool_var_3, :bool_var_4, :bool_var_5]
        elseif (terminal_node_index == 16)
            values = [:int_var_1, :int_var_2, :int_var_3, :int_var_4, :int_var_5]
        elseif (terminal_node_index == 17)
            values = [:float_var_1, :float_var_2, :float_var_3, :float_var_4, :float_var_5]
        elseif (terminal_node_index == 18)
            values = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
        elseif (terminal_node_index == 19)
            values = [(0, 1), (0, 0.5)]
        elseif (terminal_node_index == 20)
            values = [(0, 1), (0, 0.5)]
        end
        values[production_rule_index]
    else
        throw(ArgumentError(""))
    end
end

""" Returns ordered node indices associated with node_index and production_rule_index """
function getNodeIndexFromSymbol(symbol)
    symbols_to_node_indices[symbol]
end

function getSymbolsOfProductionRule(node_index, production_rule_index)
    nonterminal_indices_to_child_symbols[node_index][production_rule_index]
end

""" Returns randomly sampled production_rule_index for node_index and probability """
function getProductionRuleIndexAndProb(node_index)
    if node_index in [12]
        production_rule_indices = [1]
    elseif node_index in [1, 2, 3, 8, 9, 10, 13, 19, 20, 21]
        production_rule_indices = [1, 2]
    elseif node_index in [4, 6, 7, 11]
        production_rule_indices = [1,2,3]
    elseif node_index in [14, 5]
        production_rule_indices = [1,2,3,4]
    elseif node_index in [15, 16, 17]
        production_rule_indices = [1,2,3,4,5]
    elseif node_index in [18]
        production_rule_indices = [1,2,3,4,5,6,7,8,9,10]
    else
        throw(ArgumentError(""))
    end
    return rand(production_rule_indices), 1/length(production_rule_indices)
end

""" ----- GRAMMAR DATA STRUCTURES ----- """

nonterminal_indices_to_child_symbols = Dict(
    1 => [["line", "expr"],
          ["line"]],
    2 => [["exo_line"],
          ["endo_line"]],
    3 => [["bool_var", "bool_dist"],
          ["float_var","float_dist"]],
    4 => [["unary_endo_expr"],
          ["binary_endo_expr"],
          ["ternary_endo_expr"]],
    5 => [["bool_var", "bool_var"],
          ["bool_var", "bool_var"],
          ["int_var", "int_var"],
          ["float_var", "float_var"]],
    6 => [["bool_var", "bool_binary_op", "bool_var", "bool_var_list"],
          ["int_var", "num_binary_op", "int_var", "int_var_list"],
          ["float_var", "num_binary_op", "float_var", "float_var_list"]],
    7 => [["bool_var", "bool_var", "bool_var", "bool_var"],
          ["int_var", "bool_var", "int_var", "int_var"],
          ["float_var", "bool_var", "float_var", "float_var"]],
    8 => [["bool_var", "bool_var_list"],
          ["bool_var"]],
    9 => [["int_var", "int_var_list"],
          ["int_var"]],
    10 => [["float_var", "float_var_list"],
          ["float_var"]],
    11 => [["bool_var"],
           ["int_var"],
           ["float_var"]],
    12 => [["bernoulli_params"]],
    21 => [["normal_params"], ["uniform_params"]]
)

symbols_to_node_indices = Dict(
    "expr" => 1,
    "line" => 2,
    "exo_line" => 3,
    "endo_line" => 4,
    "unary_endo_expr" => 5,
    "binary_endo_expr" => 6,
    "ternary_endo_expr" => 7,
    "bool_var_list" => 8,
    "int_var_list" => 9,
    "float_var_list" => 10,
    "var" => 11,
    "bool_dist" => 12,
    "bool_binary_op" => 13,
    "num_binary_op" => 14,
    "bool_var" => 15,
    "int_var" => 16,
    "float_var" => 17,
    "bernoulli_params" => 18,
    "normal_params" => 19,
    "uniform_params" => 20,
    "float_dist" => 21,
)

end


text:
factorial2 = n -> n == 1 ? 1 : n * factorial2(n-1)
binomial2 = (n, p) -> sum(p .> rand(n))
pi2 = sum([x^2 + y^2 <= 1 for (x, y) in zip(rand(10000), rand(10000))])/10000*4

function pay()
  s = 0
  for i in 1:10
    1/2 > rand() ? s += 1 : s = 0
    s >= 3 && return true
  end
  return false
end

#println(pi2)
#println(pay())


text:
lines(s) = split(s, "\n")

function memoize(f)
  mem = d()
  (args...) -> Base.@get!(mem, args, f(args...))
end

function memoize_debounce(f, expiry = 1)
  mem = d()
  function (args...)
    if haskey(mem, args)
      result, t = mem[args]
      if t + expiry > time()
        mem[args] = result, time()
        return result
      end
    end
    result = f(args...)
    mem[args] = result, time()
    return result
  end
end


text:
include("../paths.jl")
using Policy

function testgradll(p, n, m, N)
   dt     = eltype(p.theta)
   feat   = rand(dt, n, N)
   act    = rand(dt, m, N)
   gradll = zeros(dt, p.nparam, N)

   @time Policy.gradLL!(gradll, p, feat, act)
   return gradll
end

#function testll(p, n, m, N)
#   dt   = eltype(p.theta)
#   feat = rand(dt, n, N)
#   act  = rand(dt, m, N)
#
#   ll   = zeros(dt, 1, N)
#
#   @time Policy.loglikelihood!(ll, p.theta, p, feat, act)
#end

function testclone(p, n, m, N)
   dt   = eltype(p.theta)
   feat = randn(dt, n, N)
   A = randn(m, n)
   b = rand(m)
   act  = A*feat .+ b

   @time Policy.clone(p, feat, act; epochs=30)
   act2 = p.model(feat).data

   println("BC:")
   println(mean(mean(act2-act, 2)))

   println(A)
   println(b)
end

n = 30
m = 20
N = 10000
hidden = 64

dt = Float64
#for dt in [Float32, Float64]
#   info(dt)
#   #glp = Policy.FLP{dt}(n, m)
#   nn  = Policy.FNN{dt}(n, m, 64, 2)
#
#   info("Gradll")
#   #info("GLP")
#   #testgradll(glp, n, m, 12)
#   #testgradll(glp, n, m, N)
#   info("NN")
#   testgradll(nn, n, m, 12)
#   testgradll(nn, n, m, N)
#
#   info("LL")
#   #info("GLP")
#   #testll(glp, n, m, 12)
#   #testll(glp, n, m, N)
#   info("NN")
#   testll(nn, n, m, 12)
#   testll(nn, n, m, N)
#end

feat = rand(dt, n, N)
act  = rand(dt, m, N)
ll   = zeros(dt, 1, N)
#glp = Policy.GLP{dt}(n, m)
fnn  = Policy.NN{dt}(n, m, hidden)
nn   = Policy.GLP{dt}(n, m)
#nn  = Policy.NN{dt}(n, m, hidden, 2)

#Policy.updatetheta(nn, fnn.theta)

#g1 = testgradll(fnn, n, m, N);
#
#printlin("GradLL diff: $(mean(g2-g1))")

testclone(fnn, n, m, N)


text:
using Manifolds, Manopt, Test, ManifoldsBase

using Random
Random.seed!(42)
# Test the additional manifold functions
#
@testset "Additional Manifold functions" begin
    @testset "mid point & reflect" begin
        M = Sphere(2)
        p = [1.0, 0.0, 0.0]
        q = [0.0, 1.0, 0.0]

        r = mid_point(M,p,q)
        r2 = similar(r)
        mid_point!(M,r2,p,q)
        r3 = shortest_geodesic(M,p,q,0.5)
        r4 = mid_point(M,p,q,q)
        @test isapprox(M,r,r2)
        @test isapprox(M,r2,r3)
        @test isapprox(M,r3,r4)
        r5 = similar(r4)
        mid_point!(M,r5,p,q,q)
        @test isapprox(M,r4,r5)

        r4 = mid_point(M,p,-p,q)
        r5 = similar(r4)
        mid_point!(M,r5,p,-p,q)
        @test isapprox(M,r4,q)
        @test isapprox(M,r4,r5)

        @test isapprox(M,reflect(M,p,q),-q)
        f = x->x
        @test reflect(M,f,q) == q

        M2 = Euclidean(2)
        p2 = [1.0, 0.0]
        q2 = [0.0, 1.0]
        s = mid_point(M2,p2,q2)
        s2 = similar(s)
        mid_point!(M2,s2,p2,q2)
        @test s==s2
        @test s==(p2+q2)/2
        s = mid_point(M2,p2,q2,s)
        @test s==s2
        s2 = similar(s)
        mid_point!(M2,s2,p2,q2,s)
        @test s==s2

    end
    @testset "random" begin
        Mc = Circle()
        pc = random_point(Mc)
        @test is_manifold_point(Mc, pc, true)
        Xc = random_tangent(Mc,pc)
        @test is_tangent_vector(Mc, pc, Xc, true)

        Me = Euclidean(3)
        pe = random_point(Me)
        @test is_manifold_point(Me, pe, true)
        Xe = random_tangent(Me,pe)
        @test is_tangent_vector(Me, pe, Xe, true)

        Mg = Grassmann(3,2)
        pg = random_point(Mg)
        @test is_manifold_point(Mg, pg, true)
        Xg = random_tangent(Mg,pg)
        @test is_tangent_vector(Mg, pg, Xg, true; atol=10^(-14))

        Mp = ProductManifold(Mg,Me)
        pp = random_point(Mp)
        @test is_manifold_point(Mp, pp, true)
        Xp = random_tangent(Mp,pp)
        @test is_tangent_vector(Mp, pp, Xp,true; atol=10^(-15))

        Mp2 = PowerManifold(Me, NestedPowerRepresentation(),4)
        pp2 = random_point(Mp2)
        @test is_manifold_point(Mp2,pp2)
        Xp2 = random_tangent(Mp2,pp2)
        @test is_tangent_vector(Mp2,pp2,Xp2)

        Mp3 = PowerManifold(Me, ArrayPowerRepresentation(), 4)
        pp3 = random_point(Mp3)
        @test is_manifold_point(Mp3,pp3)
        Xp3 = random_tangent(Mp3,pp3)
        @test is_tangent_vector(Mp3,pp3,Xp3)

        Mr = Rotations(3)
        pr = random_point(Mr)
        @test is_manifold_point(Mr,pr)
        Xr = random_tangent(Mr,pr)
        @test is_tangent_vector(Mr,pr,Xr)

        Mspd = SymmetricPositiveDefinite(3)
        pspd = random_point(Mspd)
        @test is_manifold_point(Mspd,pspd;atol=10^(-14))
        Xspd = random_tangent(Mspd,pspd)
        @test is_tangent_vector(Mspd,pspd,Xspd, true; atol=10^(-15))
        Xspd2 = random_tangent(Mspd,pspd, Val(:Rician))
        @test is_tangent_vector(Mspd,pspd,Xspd2,true;atol=10^(-15))

        Mst = Stiefel(3,2)
        pst = random_point(Mst)
        @test is_manifold_point(Mst, pst, true)

        Msp = Sphere(2)
        psp = random_point(Msp)
        @test is_manifold_point(Msp, psp, true)
        Xsp = random_tangent(Msp,psp)
        @test is_tangent_vector(Msp,psp,Xsp,true;atol=10^(-15))

        Mh = Hyperbolic(2)
        ph = [0.0,0.0,1.0]
        Xh = random_tangent(Mh, ph)
        @test is_tangent_vector(Mh, ph, Xh, true)
    end
end


text:
# module Bukdu

export routes

export get, post, delete, patch, put
import Base: get
import .System: HaltedError, NotApplicableError, InternalError, SystemController, internal_error, halted_error, not_applicable, info_response

struct DirectRequest
    _req
end

const routing_verbs = [:get, :post, :delete, :patch, :put]

# HTTP.URIs: queryparams
# HTTP.jl - URIs.jl
function queryparams(q::AbstractString)
    Dict(Deps.HTTP.URIs.unescapeuri(k) => Deps.HTTP.URIs.unescapeuri(v)
        for (k,v) in ([split(e, "=")..., ""][1:2]
            for e in split(q, "&", keepempty=false)))
end

function fetch_query_params(req::Deps.Request)::Vector{Pair{String,String}}
    params = queryparams(Deps.HTTP.URIs.URI(req.target).query)
    collect(params)
end

function _build_conn_and_pipelines(route::Route, req::Deps.Request)
    body_params = Plug.Parsers.fetch_body_params(req)
    query_params = fetch_query_params(req)
    path_params = route.path_params
    params = merge(body_params, query_params, path_params)
    halted = false
    conn = Conn(req, Assoc.((body_params, query_params, path_params, params))..., halted)
    for pipefunc in route.pipelines
        pipefunc(conn)
        halted = conn.halted
        halted && break
    end
    if halted
        err = HaltedError("halted on pipelines")
        rou = Route(SystemController, halted_error, route.path_params, route.pipelines)
        obj = halted_error(SystemController(conn, err))
        (rou, obj)
    else
        controller = route.C(conn)
        if applicable(route.action, controller)
            obj = route.action(controller)
            (route, obj)
        else
            err = NotApplicableError(string(route.action, "(::", route.C, ")"))
            rou = Route(SystemController, not_applicable, route.path_params, route.pipelines)
            obj = not_applicable(SystemController(conn, err))
            (rou, obj)
        end
    end
end

function _catch_internal_error(block, route, req)
    try
        block(route, req)
    catch ex
        stackframes = stacktrace(catch_backtrace())
        err = InternalError(ex, stackframes)
        conn = Conn(req)
        rou = Route(SystemController, internal_error, route.path_params, route.pipelines)
        obj = internal_error(SystemController(conn, err))
        (rou, obj)
    end
end

function _proc_request(route::Route, req::Deps.Request)
    System.catch_request(route, req)           # System
    req.response.status = 200
    _catch_internal_error(route, req) do route, req
    # begin #
        _build_conn_and_pipelines(route, req)
    end
end

function _proc_response(route::Route, req::Deps.Request)
    info_response(route, req, req.response)
    System.catch_response(route, req.response) # System
end

function request_handler(route::Route, dreq::DirectRequest)
    (rou, obj) = _proc_request(route, dreq._req)
    _proc_response(rou, dreq._req)
    (got=obj, resp=dreq._req.response, route=rou)
end

function request_handler(route::Route, req::Deps.Request)
    (rou, obj) = _proc_request(route, req)
    if obj isa AbstractRender
        data = obj.body
        push!(req.response.headers, Pair("Content-Type", obj.content_type))
    else
        data = unsafe_wrap(Vector{UInt8}, string(obj))
    end
    req.response.body = data
    _proc_response(rou, req)
    (got=obj, resp=req.response, route=rou)
end

"""
    routes(block::Function)
"""
function routes(block::Function)
    block()
end

"""
    routes(block::Function, pipe::Symbol)
"""
function routes(block::Function, pipe::Symbol)
    Routing.store[:pipe] = pipe
    block()
    Routing.store[:pipe] = nothing
end

"""
    get(url::String, C::Type{<:ApplicationController}, action)
"""
function get
end

"""
    post(url::String, C::Type{<:ApplicationController}, action)
"""
function post
end

"""
    delete(url::String, C::Type{<:ApplicationController}, action)
"""
function delete
end

"""
    patch(url::String, C::Type{<:ApplicationController}, action)
"""
function patch
end

"""
    put(url::String, C::Type{<:ApplicationController}, action)
"""
function put
end

for verb in routing_verbs
    @eval ($verb)(url::String, C::Type{<:ApplicationController}, action) = Routing.add_route($verb, url, C, action)
end

# module Bukdu

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
__precompile__(true)
module OctreeBuilder

using Images
using Printf

export build_octree

# load image from filename
load_data(filename) = rawview(channelview(load(filename)))

# save image to filename
function save_data(filename, data)
    @printf "Writing %s\n" filename
    save(filename, data)
end

# 2nd brightest of the 8 pixels
# equivalent to sort(vec(arg))[7] but half the time and a third the memory usage
function downsampling_function(arg::Array{UInt16,3})
    m0::UInt16 = 0x0000
    m1::UInt16 = 0x0000
    for i = 1:8
        @inbounds tmp::UInt16 = arg[i]
        if tmp > m0
            m1 = m0
            m0 = tmp
        elseif tmp > m1
            m1 = tmp
        end
    end
    m1
end

function downsample!(out_tile_jl, coord, shape_leaf_px, scratch)
    iy = ((coord - 1) >> 1) & 1 * shape_leaf_px[1] >> 1
    ix = ((coord - 1) >> 0) & 1 * shape_leaf_px[2] >> 1
    iz = ((coord - 1) >> 2) & 1 * shape_leaf_px[3] >> 1
    for z = 1:2:shape_leaf_px[3] - 1
        tmpz = iz + (z + 1) >> 1
        for x = 1:2:shape_leaf_px[2] - 1
            tmpx = ix + (x + 1) >> 1
            for y = 1:2:shape_leaf_px[1] - 1
                tmpy = iy + (y + 1) >> 1
                out_tile_jl[tmpy, tmpx, tmpz] = downsampling_function(scratch[y:y + 1, x:x + 1, z:z + 1])
            end
        end
    end
end

function build_octree(sourcefilename, targetpath, nlevels, channel)
    @printf "Loading image from %s\n" sourcefilename
    img = load_data(sourcefilename)

    @printf "Will generate octree with %d levels to %s\n" nlevels targetpath
    @printf "Image size: %s\n" size(img)

    # crop image to a power-of-two multiple of the octree leaf size
    dim = [size(img)[1:3]...]
    while(rem(dim[1], 2^(nlevels - 1)) > 0)  dim[1] -= 1;  end
    while(rem(dim[2], 2^(nlevels - 1)) > 0)  dim[2] -= 1;  end
    while(rem(dim[3], 2^(nlevels - 1)) > 0)  dim[3] -= 1;  end
    dim_leaf = [x >> (nlevels - 1) for x in dim]

    @printf "Adjusted image size: %s, Dim leaf: %s\n" dim dim_leaf

    function octree_division(relpath, img_down::Array{T,3}, img_view::SubArray{T,3}) where T <: UInt16
        morton = split(relpath, Base.Filesystem.path_separator)
        level = relpath == "" ? 1 : length(morton) + 1
        @printf "Level: %d, Current path: %s/%s\n" level targetpath relpath
        mkpath(joinpath(targetpath, relpath))
        if level < nlevels
            img_down_next = Array{T}(undef, dim_leaf...)
            dim_view = size(img_view)[1:3]
            for z = 1:2, y = 1:2, x = 1:2
                octant_path = string(x + 2 * (y - 1) + 4 * (z - 1))
                start_x = x == 1 ? 1 : dim_view[2] >> 1 + 1
                end_x = x == 1 ? dim_view[2] >> 1 : dim_view[2]
                start_y = y == 1 ? 1 : dim_view[1] >> 1 + 1
                end_y = y == 1 ? dim_view[1] >> 1 : dim_view[1]
                start_z = z == 1 ? 1 : dim_view[3] >> 1 + 1
                end_z = z == 1 ? dim_view[3] >> 1 : dim_view[3]

                @printf "octant:(%d,%d,%d) -> %s/%s (%d:%d, %d:%d, %d:%d)\n" x y z relpath octant_path start_x end_x start_y end_y start_z end_z
                octree_division(joinpath(relpath, octant_path), img_down_next, view(img_view, start_y:end_y, start_x:end_x, start_z:end_z))
            end
        end
        saveme = level == nlevels ? img_view : img_down_next
        filename = @sprintf "default.%d.tif" channel
        save_data(joinpath(targetpath, relpath, filename), saveme)
        if (level > 1)
            downsample!(img_down, parse(Int, morton[end]), dim_leaf, saveme)
            nothing
        end
    end

    octree_division("", Array{eltype(img),3}(undef, dim_leaf...), view(img, 1:dim[1], 1:dim[2], 1:dim[3]))
end

function write_result(targetpath, voxelsize_um, nlevels)
    @printf "Write transform.txt to %s\n" targetpath
    fid = open(joinpath(targetpath, "transform.txt"), "w")
    println(fid, "ox: 10000")
    println(fid, "oy: 10000")
    println(fid, "oz: 10000")
    println(fid, "sx: ", voxelsize_um[1] * 2^(nlevels - 1))
    println(fid, "sy: ", voxelsize_um[2] * 2^(nlevels - 1))
    println(fid, "sz: ", voxelsize_um[3] * 2^(nlevels - 1))
    println(fid, "nl: ", nlevels)
    close(fid)
end

greet() = "Hello"

end # module

max_stars_repo_name: JaneliaSciComp/jacs-tools-docker

max_issues_repo_name: JaneliaSciComp/jacs-tools-docker

max_forks_repo_name: JaneliaSciComp/jacs-tools-docker


text:
module QGA
export minimize

using Printf
#using StatsBase
using Random
using Statistics
using LinearAlgebra

function ent(t, f, fmin, w)
    w .= exp2.(t.*(fmin.-f))
    z = sum(w)
    w .*= t.*(fmin.-f)
    log2(z)-sum(w)/z
end

function selection_entropy!(w, k, t, f, fmin)
    if ent(t, f, fmin, w) > k
        while true
            t = t*1.001
            if ent(t, f, fmin, w) < k
                return t
            end
#             if t > 1e15
#                 @warn "probably infinite loop" f
#             end
        end
    else
        smax = log2(length(f))
        while true
            t = t/1.001
            s = ent(t, f, fmin, w)
            if s > k
                return t
            elseif s >= smax - 1e-2
                @warn "reached maximum entropy $smax"
                return t
            end
        end
    end
end

function minimize(fitness, x, f, w, r, xm, xnew, s, n,
        FminTarget, sdmin, trace, MaxEvals, mumin; evals = 0)
    fmin, fmini = findmin(f)
    beta = 0.1/std(f)
    sd = 0.0
    dim = size(x, 1)
    local halt
    if trace < Inf
        println("Iteration     Fmin     Beta")
    end
    while true
        beta = selection_entropy!(w, s, beta, f, fmin)
        w .= exp2.(beta.*(fmin.-f))
        z = sum(w)
        z2 = sum(abs2, w)
        fm = mapreduce(*, +, f, w)/z
        if sdmin > 0.0
            sd = mapreduce((f,w)-> (f-fm)*w, +, f, w)/(z - z2/z)
            if sd < sdmin
                halt = :sdmin
                break
            end
        end
        if evals % trace == 0
            @printf("%9i %1.2e %1.2e\n",
                    evals, fmin, beta)
        end
        if mumin
            xm .= x[:, fmini]
        else
            for i = 1:dim
                xm[i] = 0.0
                for j = 1:n
                    xm[i] += x[i, j]*w[j]
                end
                xm[i] /= z
            end
        end

        randn!(r)
        r .= r.*sqrt.(w)/sqrt(z-z2/z)
        mul!(xnew, x, r)
        xnew .+= xm.*(1-sum(r))

        j = argmin(w)
        f[j] = fitness(xnew)
        x[:, j] .= xnew
        evals += 1
        if f[j] < fmin
            fmin = f[j]
            fmini = j
        end
        if fmin < FminTarget
            halt = :FminTarget
            break
        elseif evals >= MaxEvals
            halt = :MaxEvals
            break
        elseif sum(isequal(f[j]), f) > 1
            halt = :DupFit
            break
        end
    end
    Dict(:fmin => fmin, :halt => halt, :xmin => x[:,fmini],
         :n => n, :mumin => mumin, :sd => sd,
         :evals => evals, :s=> s)
end

function minimize(fitness, xm0, xsd0, s;
        FminTarget = -Inf, sdmin = 1e-16,
        trace = Inf, MaxEvals = 10000*length(xm0),
        mumin = true,
        n = floor(Int64, exp2(s+1)))
    dim = length(xm0)
    x = randn(dim, n).*xsd0 .+ xm0
    f = vec(mapslices(fitness, x; dims=1))
    w = zeros(n)
    r = zeros(n)
    xm = zeros(dim)
    xnew = zeros(dim)
    minimize(fitness, x, f, w, r, xm, xnew, s, n,
        FminTarget, sdmin, trace, MaxEvals, mumin, evals = n)
end

end

max_stars_repo_name: jotwin/qga

max_issues_repo_name: jotwin/qga

max_forks_repo_name: jotwin/qga


text:
using BinaryProvider
# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = Product[
    # Instantiate products here, e.g.:
    LibraryProduct(prefix, "libgumbo", :libgumbo),
 ]
# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaWeb/GumboBuilder/releases/download/v0.1.0"
# Listing of files generated by BinaryBuilder:
download_info = Dict(
    BinaryProvider.Linux(:aarch64, :glibc) => ("$bin_prefix/Gumbo.aarch64-linux-gnu.tar.gz", "14e50e3705352ead8e189375c93886028b3a947e61e31b0e8f8068d8f6a7e84b"),
    BinaryProvider.Linux(:armv7l, :glibc) => ("$bin_prefix/Gumbo.arm-linux-gnueabihf.tar.gz", "d76c42f3f3a8f75ab8fa976f8e309bc7bd30e2f2188844df785c4b36462545fb"),
    BinaryProvider.Linux(:i686, :glibc) => ("$bin_prefix/Gumbo.i686-linux-gnu.tar.gz", "0d1de0621fc6f23ae3119f2fc8102f0af83bc1f82287e943d23159f60ab0c503"),
    BinaryProvider.Windows(:i686) => ("$bin_prefix/Gumbo.i686-w64-mingw32.tar.gz", "48f80a41d504775a7d79acffd597d5921bf3fd14f2d046f229f643c43656c9ce"),
    BinaryProvider.Linux(:powerpc64le, :glibc) => ("$bin_prefix/Gumbo.powerpc64le-linux-gnu.tar.gz", "7b7862ee7d0192dd139e75342c00993cb78c9b5b1237bf0d20a42fde855fcaf0"),
    BinaryProvider.MacOS() => ("$bin_prefix/Gumbo.x86_64-apple-darwin14.tar.gz", "717b52a0361a70e7f7637ad5cbc970f8c3648777cf382e609631319adef6abc0"),
    BinaryProvider.Linux(:x86_64, :glibc) => ("$bin_prefix/Gumbo.x86_64-linux-gnu.tar.gz", "6cf72df5dbbea1de1e0b777c70f88b9ac41a0ef728a22940aed3b910713d3ccc"),
    BinaryProvider.Windows(:x86_64) => ("$bin_prefix/Gumbo.x86_64-w64-mingw32.tar.gz", "7a2994a966ae4a408c34df408a9721760836e0115839dc42e6c3b140f9b509a8"),
)
# First, check to see if we're all satisfied
if any(!satisfied(p; verbose=verbose) for p in products)
    if haskey(download_info, platform_key())
        # Download and install binaries
        url, tarball_hash = download_info[platform_key()]
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    else
        # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
        # Alternatively, you could attempt to install from a separate provider,
        # build from source or something more even more ambitious here.
        error("Your platform $(Sys.MACHINE) is not supported by this package!")
    end
end
# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products)