text:
"""
    pc = finish!(recurse, frame, istoplevel=false)
    pc = finish!(frame, istoplevel=false)

Run `frame` until execution terminates. `pc` is either `nothing` (if execution terminates
when it hits a `return` statement) or a reference to a breakpoint.
In the latter case, `leaf(frame)` returns the frame in which it hit the breakpoint.

`recurse` controls call evaluation; `recurse = Compiled()` evaluates :call expressions
by normal dispatch, whereas the default `recurse = finish_and_return!` uses recursive interpretation.
"""
function finish!(@nospecialize(recurse), frame::Frame, istoplevel::Bool=false)
    while true
        pc = step_expr!(recurse, frame, istoplevel)
        (pc === nothing || isa(pc, BreakpointRef)) && return pc
        shouldbreak(frame, pc) && return BreakpointRef(frame.framecode, pc)
    end
end
finish!(frame::Frame, istoplevel::Bool=false) = finish!(finish_and_return!, frame, istoplevel)

"""
    ret = finish_and_return!(recurse, frame, istoplevel::Bool=false)
    ret = finish_and_return!(frame, istoplevel::Bool=false)

Call [`JuliaInterpreter.finish!`](@ref) and pass back the return value `ret`. If execution
pauses at a breakpoint, `ret` is the reference to the breakpoint.
"""
function finish_and_return!(@nospecialize(recurse), frame::Frame, istoplevel::Bool=false)
    pc = finish!(recurse, frame, istoplevel)
    isa(pc, BreakpointRef) && return pc
    return get_return(frame)
end
finish_and_return!(frame::Frame, istoplevel::Bool=false) = finish_and_return!(finish_and_return!, frame, istoplevel)

"""
    bpref = dummy_breakpoint(recurse, frame::Frame, istoplevel)

Return a fake breakpoint. `dummy_breakpoint` can be useful as the `recurse` argument to
`evaluate_call!` (or any of the higher-order commands) to ensure that you return immediately
after stepping into a call.
"""
dummy_breakpoint(@nospecialize(recurse), frame::Frame, istoplevel) = BreakpointRef(frame.framecode, 0)

"""
    ret = finish_stack!(recurse, frame, rootistoplevel=false)
    ret = finish_stack!(frame, rootistoplevel=false)

Unwind the callees of `frame`, finishing each before returning to the caller.
`frame` itself is also finished. `rootistoplevel` should be true if the root frame is top-level.

`ret` is typically the returned value. If execution hits a breakpoint, `ret` will be a
reference to the breakpoint.
"""
function finish_stack!(@nospecialize(recurse), frame::Frame, rootistoplevel::Bool=false)
    frame0 = frame
    frame = leaf(frame)
    while true
        istoplevel = rootistoplevel && frame.caller === nothing
        ret = finish_and_return!(recurse, frame, istoplevel)
        isa(ret, BreakpointRef) && return ret
        frame === frame0 && return ret
        frame = return_from(frame)
        frame === nothing && return ret
        pc = frame.pc
        if isassign(frame, pc)
            lhs = SSAValue(pc)
            do_assignment!(frame, lhs, ret)
        else
            stmt = pc_expr(frame, pc)
            if isexpr(stmt, :(=))
                lhs = stmt.args[1]
                do_assignment!(frame, lhs, ret)
            end
        end
        pc += 1
        frame.pc = pc
        shouldbreak(frame, pc) && return BreakpointRef(frame.framecode, pc)
    end
end
finish_stack!(frame::Frame, istoplevel::Bool=false) = finish_stack!(finish_and_return!, frame, istoplevel)

"""
    pc = next_until!(predicate, recurse, frame, istoplevel=false)
    pc = next_until!(predicate, frame, istoplevel=false)

Execute the current statement. Then step through statements of `frame` until the next
statement satisfies `predicate(frame)`. `pc` will be the index of the statement at which
evaluation terminates, `nothing` (if the frame reached a `return`), or a `BreakpointRef`.
"""
function next_until!(@nospecialize(predicate), @nospecialize(recurse), frame::Frame, istoplevel::Bool=false)
    pc = step_expr!(recurse, frame, istoplevel)
    while pc !== nothing && !isa(pc, BreakpointRef)
        shouldbreak(frame, pc) && return BreakpointRef(frame.framecode, pc)
        predicate(frame) && return pc
        pc = step_expr!(recurse, frame, istoplevel)
    end
    return pc
end
next_until!(predicate, frame::Frame, istoplevel::Bool=false) =
    next_until!(predicate, finish_and_return!, frame, istoplevel)

"""
    pc = maybe_next_until!(predicate, recurse, frame, istoplevel=false)
    pc = maybe_next_until!(predicate, frame, istoplevel=false)

Like [`next_until!`](@ref) except checks `predicate` before executing the current statment.

"""
function maybe_next_until!(@nospecialize(predicate), @nospecialize(recurse), frame::Frame, istoplevel::Bool=false)
    predicate(frame) && return frame.pc
    return next_until!(predicate, recurse, frame, istoplevel)
end
maybe_next_until!(@nospecialize(predicate), frame::Frame, istoplevel::Bool=false) =
    maybe_next_until!(predicate, finish_and_return!, frame, istoplevel)

"""
    pc = next_call!(recurse, frame, istoplevel=false)
    pc = next_call!(frame, istoplevel=false)

Execute the current statement. Continue stepping through `frame` until the next
`:return` or `:call` expression.
"""
next_call!(@nospecialize(recurse), frame::Frame, istoplevel::Bool=false) =
    next_until!(frame -> is_call_or_return(pc_expr(frame)), recurse, frame, istoplevel)
next_call!(frame::Frame, istoplevel::Bool=false) = next_call!(finish_and_return!, frame, istoplevel)

"""
    pc = maybe_next_call!(recurse, frame, istoplevel=false)
    pc = maybe_next_call!(frame, istoplevel=false)

Return the current program counter of `frame` if it is a `:return` or `:call` expression.
Otherwise, step through the statements of `frame` until the next `:return` or `:call` expression.
"""
maybe_next_call!(@nospecialize(recurse), frame::Frame, istoplevel::Bool=false) =
    maybe_next_until!(frame -> is_call_or_return(pc_expr(frame)), recurse, frame, istoplevel)
maybe_next_call!(frame::Frame, istoplevel::Bool=false) = maybe_next_call!(finish_and_return!, frame, istoplevel)

"""
    pc = through_methoddef_or_done!(recurse, frame)
    pc = through_methoddef_or_done!(frame)

Runs `frame` at top level until it either finishes (e.g., hits a `return` statement)
or defines a new method.
"""
function through_methoddef_or_done!(@nospecialize(recurse), frame::Frame)
    predicate(frame) = (stmt = pc_expr(frame); isexpr(stmt, :method, 3) || isexpr(stmt, :thunk))
    pc = next_until!(predicate, recurse, frame, true)
    (pc === nothing || isa(pc, BreakpointRef)) && return pc
    return step_expr!(recurse, frame, true)  # define the method and return
end
through_methoddef_or_done!(@nospecialize(recurse), t::Tuple{Module,Expr,Frame}) =
    through_methoddef_or_done!(recurse, t[end])
through_methoddef_or_done!(@nospecialize(recurse), modex::Tuple{Module,Expr,Expr}) = Core.eval(modex[1], modex[3])
through_methoddef_or_done!(@nospecialize(recurse), ::Nothing) = nothing
through_methoddef_or_done!(arg) = through_methoddef_or_done!(finish_and_return!, arg)

# Sentinel to see if the call was a wrapper call
struct Wrapper end

"""
    pc = next_line!(recurse, frame, istoplevel=false)
    pc = next_line!(frame, istoplevel=false)

Execute until reaching the first call of the next line of the source code.
Upon return, `pc` is either the new program counter, `nothing` if a `return` is reached,
or a `BreakpointRef` if it encountered a wrapper call. In the latter case, call `leaf(frame)`
to obtain the new execution frame.
"""
function next_line!(@nospecialize(recurse), frame::Frame, istoplevel::Bool=false)
    pc = frame.pc
    initialline, initialfile = linenumber(frame, pc), getfile(frame, pc)
    predicate(frame) = isexpr(pc_expr(frame), :return) || (linenumber(frame) != initialline || getfile(frame) != initialfile)
    pc = next_until!(predicate, recurse, frame, istoplevel)
    (pc === nothing || isa(pc, BreakpointRef)) && return pc
    maybe_step_through_kwprep!(recurse, frame, istoplevel)
    maybe_next_call!(recurse, frame, istoplevel)
end
next_line!(frame::Frame, istoplevel::Bool=false) = next_line!(finish_and_return!, frame, istoplevel)

"""
    pc = until_line!(recurse, frame, line=nothing istoplevel=false)
    pc = until_line!(frame, line=nothing, istoplevel=false)

Execute until the current frame reaches a line greater than `line`. If `line == nothing`
execute until the current frame reaches any line greater than the current line.
"""
function until_line!(@nospecialize(recurse), frame::Frame, line::Union{Nothing, Integer}=nothing, istoplevel::Bool=false)
    pc = frame.pc
    initialline, initialfile = linenumber(frame, pc), getfile(frame, pc)
    line === nothing && (line = initialline + 1)
    predicate(frame) = isexpr(pc_expr(frame), :return) || (linenumber(frame) >= line && getfile(frame) == initialfile)
    pc = next_until!(predicate, frame, istoplevel)
    (pc === nothing || isa(pc, BreakpointRef)) && return pc
    maybe_step_through_kwprep!(recurse, frame, istoplevel)
    maybe_next_call!(recurse, frame, istoplevel)
end
until_line!(frame::Frame, line::Union{Nothing, Integer}=nothing, istoplevel::Bool=false) = until_line!(finish_and_return!, frame, line, istoplevel)

"""
    cframe = maybe_step_through_wrapper!(recurse, frame)
    cframe = maybe_step_through_wrapper!(frame)

Return the new frame of execution, potentially stepping through "wrapper" methods like those
that supply default positional arguments or handle keywords. `cframe` is the leaf frame from
which execution should start.
"""
function maybe_step_through_wrapper!(@nospecialize(recurse), frame::Frame)
    code = frame.framecode
    stmts, scope = code.src.code, code.scope::Method
    length(stmts) < 2 && return frame
    last = stmts[end-1]
    isexpr(last, :(=)) && (last = last.args[2])
    comp = VERSION < v"1.4.0-DEV.215" ? startswith : endswith
    is_kw = isa(scope, Method) && comp(String(Base.unwrap_unionall(Base.unwrap_unionall(scope.sig).parameters[1]).name.name), "#kw")
    if is_kw || isexpr(last, :call) && any(isequal(SlotNumber(1)), last.args)
        # If the last expr calls #self# or passes it to an implementation method,
        # this is a wrapper function that we might want to step through
        while frame.pc != length(stmts)-1
            pc = next_call!(recurse, frame, false)  # since we're in a Method we're not at toplevel
            if pc === nothing || isa(pc, BreakpointRef)
                return frame
            end
        end
        ret = evaluate_call!(dummy_breakpoint, frame, last)
        if !isa(ret, BreakpointRef) # Happens if next call is Compiled
            return frame
        end
        frame.framedata.ssavalues[frame.pc] = Wrapper()
        return maybe_step_through_wrapper!(recurse, callee(frame))
    end
    maybe_step_through_nkw_meta!(frame)
    return frame
end
maybe_step_through_wrapper!(frame::Frame) = maybe_step_through_wrapper!(finish_and_return!, frame)

"""
    frame = maybe_step_through_kwprep!(recurse, frame)
    frame = maybe_step_through_kwprep!(frame)

If `frame.pc` points to the beginning of preparatory work for calling a keyword-argument
function, advance forward until the actual call.
"""
function maybe_step_through_kwprep!(@nospecialize(recurse), frame::Frame, istoplevel::Bool=false)
    pc, src = frame.pc, frame.framecode.src
    n = length(src.code)
    stmt = pc_expr(frame, pc)
    if isa(stmt, Tuple{Symbol,Vararg{Symbol}})
        # Check to see if we're creating a NamedTuple followed by kwfunc call
        pccall = pc + 5
        if pccall <= n
            stmt1 = src.code[pc+1]
            # We deliberately check isexpr(stmt, :call) rather than is_call(stmt): if it's
            # assigned to a local, it's *not* kwarg preparation.
            if isexpr(stmt1, :call) && is_quotenode(stmt1.args[1], Core.apply_type) && is_quoted_type(stmt1.args[2], :NamedTuple)
                stmt4, stmt5 = src.code[pc+4], src.code[pc+5]
                if isexpr(stmt4, :call) && is_quotenode(stmt4.args[1], Core.kwfunc)
                    while pc < pccall
                        pc = step_expr!(recurse, frame, istoplevel)
                    end
                    return frame
                elseif isexpr(stmt5, :call) && is_quotenode(stmt5.args[1], Core.kwfunc) && pccall+1 <= n
                    # This happens when the call is scoped by a module
                    pccall += 1
                    while pc < pccall
                        pc = step_expr!(recurse, frame, istoplevel)
                    end
                    maybe_next_call!(recurse, frame, istoplevel)
                    return frame
                end
            end
        end
    elseif isexpr(stmt, :call) && is_quoted_type(stmt.args[1], :NamedTuple) && length(stmt.args) == 1
        # Creating an empty NamedTuple, now split by type (no supplied kwargs vs kwargs...)
        if pc + 1 <= n
            stmt1 = src.code[pc+1]
            if isexpr(stmt1, :call)
                f = stmt1.args[1]
                if is_quotenode(f, Base.pairs)
                    # No supplied kwargs
                    pccall = pc + 2
                    if pccall <= n
                        stmt2 = src.code[pccall]
                        if isexpr(stmt2, :call) && length(stmt2.args) >= 3 && stmt2.args[2] == SSAValue(pc+1) && stmt2.args[3] == SlotNumber(1)
                            while pc < pccall
                                pc = step_expr!(recurse, frame, istoplevel)
                            end
                        end
                    end
                elseif is_quotenode(f, Base.merge) && ((pccall = pc + 7) <= n)
                    stmtk = src.code[pccall-1]
                    if isexpr(stmtk, :call) && is_quotenode(stmtk.args[1], Core.kwfunc)
                        for i = 1:4
                            pc = step_expr!(recurse, frame, istoplevel)
                        end
                        stmti = src.code[pc]
                        if isexpr(stmti, :call) && is_quotenode(stmti.args[1], Core.kwfunc)
                            pc = step_expr!(recurse, frame, istoplevel)
                        end
                    end
                end
            end
        end
    end
    return frame
end
maybe_step_through_kwprep!(frame::Frame, istoplevel::Bool=false) =
    maybe_step_through_kwprep!(finish_and_return!, frame, istoplevel)

"""
    ret = maybe_reset_frame!(recurse, frame, pc, rootistoplevel)

Perform a return to the caller, or descend to the level of a breakpoint.
`pc` is the return state from the previous command (e.g., `next_call!` or similar).
`rootistoplevel` should be true if the root frame is top-level.

`ret` will be `nothing` if we have just completed a top-level frame. Otherwise,

    cframe, cpc = ret

where `cframe` is the frame from which execution should continue and `cpc` is the state
of `cframe` (the program counter, a `BreakpointRef`, or `nothing`).
"""
function maybe_reset_frame!(@nospecialize(recurse), frame::Frame, @nospecialize(pc), rootistoplevel::Bool)
    isa(pc, BreakpointRef) && return leaf(frame), pc
    if pc === nothing
        val = get_return(frame)
        frame = return_from(frame)
        frame === nothing && return nothing
        ssavals = frame.framedata.ssavalues
        is_wrapper = isassigned(ssavals, frame.pc) && ssavals[frame.pc] === Wrapper()
        maybe_assign!(frame, val)
        frame.pc += 1
        if is_wrapper
            return maybe_reset_frame!(recurse, frame, finish!(recurse, frame), rootistoplevel)
        end
        pc = maybe_next_call!(recurse, frame, rootistoplevel && frame.caller===nothing)
        return maybe_reset_frame!(recurse, frame, pc, rootistoplevel)
    end
    return frame, pc
end

# Unwind the stack until an exc is eventually caught, thereby
# returning the frame that caught the exception at the pc of the catch
# or rethrow the error
function unwind_exception(frame::Frame, exc)
    while frame !== nothing
        if !isempty(frame.framedata.exception_frames)
            # Exception caught
            frame.pc = frame.framedata.exception_frames[end]
            frame.framedata.last_exception[] = exc
            return frame
        end
        frame = return_from(frame)
    end
    rethrow(exc)
end

function maybe_step_through_nkw_meta!(frame)
    stmt = pc_expr(frame)
    if stmt === nothing || (isexpr(stmt, :meta) && stmt.args[1] == :nkw)
        @assert frame.pc == 1
        frame.pc += 1
    end
end

"""
    ret = debug_command(recurse, frame, cmd, rootistoplevel=false; line=nothing)
    ret = debug_command(frame, cmd, rootistoplevel=false; line=nothing)

Perform one "debugger" command. The keyword arguments are not used for all debug commands.
`cmd` should be one of:

- `:n`: advance to the next line
- `:s`: step into the next call
- `:until`: advance the frame to line `line` if given, otherwise advance to the line after the current line
- `:c`: continue execution until termination or reaching a breakpoint
- `:finish`: finish the current frame and return to the parent

or one of the 'advanced' commands

- `:nc`: step forward to the next call
- `:se`: execute a single statement
- `:si`: execute a single statement, stepping in if it's a call
- `:sg`: step into the generator of a generated function

`rootistoplevel` and `ret` are as described for [`JuliaInterpreter.maybe_reset_frame!`](@ref).
"""
function debug_command(@nospecialize(recurse), frame::Frame, cmd::Symbol, rootistoplevel::Bool=false; line=nothing)
    function nicereturn!(@nospecialize(recurse), frame, pc, rootistoplevel)
        if pc === nothing || isa(pc, BreakpointRef)
            return maybe_reset_frame!(recurse, frame, pc, rootistoplevel)
        end
        maybe_step_through_kwprep!(recurse, frame, rootistoplevel && frame.caller === nothing)
        return frame, frame.pc
    end

    istoplevel = rootistoplevel && frame.caller === nothing
    cmd0 = cmd
    is_si = false
    if cmd == :si
        stmt = pc_expr(frame)
        cmd = is_call(stmt) ? :s : :se
        is_si = true
    end
    try
        cmd == :nc && return nicereturn!(recurse, frame, next_call!(recurse, frame, istoplevel), rootistoplevel)
        cmd == :n && return maybe_reset_frame!(recurse, frame, next_line!(recurse, frame, istoplevel), rootistoplevel)
        cmd == :se && return maybe_reset_frame!(recurse, frame, step_expr!(recurse, frame, istoplevel), rootistoplevel)
        cmd == :until && return maybe_reset_frame!(recurse, frame, until_line!(recurse, frame, line, istoplevel), rootistoplevel)

        enter_generated = false
        if cmd == :sg
            enter_generated = true
            cmd = :s
        end
        if cmd == :s
            pc = maybe_next_call!(recurse, frame, istoplevel)
            (isa(pc, BreakpointRef) || pc === nothing) && return maybe_reset_frame!(recurse, frame, pc, rootistoplevel)
            is_si || maybe_step_through_kwprep!(recurse, frame, istoplevel)
            pc = frame.pc
            stmt0 = stmt = pc_expr(frame, pc)
            isexpr(stmt0, :return) && return maybe_reset_frame!(recurse, frame, nothing, rootistoplevel)
            if isexpr(stmt, :(=))
                stmt = stmt.args[2]
            end
            local ret
            try
                ret = evaluate_call!(dummy_breakpoint, frame, stmt; enter_generated=enter_generated)
            catch err
                ret = handle_err(recurse, frame, err)
                return isa(ret, BreakpointRef) ? (leaf(frame), ret) : ret
            end
            if isa(ret, BreakpointRef)
                newframe = leaf(frame)
                cmd0 == :si && return newframe, ret
                is_si || (newframe = maybe_step_through_wrapper!(recurse, newframe))
                is_si || maybe_step_through_kwprep!(recurse, newframe, istoplevel)
                return newframe, BreakpointRef(newframe.framecode, 0)
            end
            # if we got here, the call returned a value
            maybe_assign!(frame, stmt0, ret)
            frame.pc += 1
            return frame, frame.pc
        end
        if cmd == :c
            r = root(frame)
            ret = finish_stack!(recurse, r, rootistoplevel)
            return isa(ret, BreakpointRef) ? (leaf(r), ret) : nothing
        end
        cmd == :finish && return maybe_reset_frame!(recurse, frame, finish!(recurse, frame, istoplevel), rootistoplevel)
    catch err
        frame = unwind_exception(frame, err)
        if cmd == :c
            return debug_command(recurse, frame, :c, istoplevel)
        else
            return debug_command(recurse, frame, :nc, istoplevel)
        end
    end
    throw(ArgumentError("command $cmd not recognized"))
end
debug_command(frame::Frame, cmd::Symbol, rootistoplevel::Bool=false; kwargs...) =
    debug_command(finish_and_return!, frame, cmd, rootistoplevel; kwargs...)


text:
"""
    abstract type LatticeGeometry

A `LatticeGeometry` controls which sites in an [`AbstractFockAddress`](@ref) are considered
to be neighbours.

Currently only supported by [`HubbardRealSpace`](@ref).

# Available implementations

* [`PeriodicBoundaries`](@ref)
* [`HardwallBoundaries`](@ref)
* [`LadderBoundaries`](@ref)

# Interface to implement

* `Base.size`: return the lattice size.
* [`neighbour_site(::LatticeGeometry, ::Int, ::Int)`](@ref)
* [`num_neighbours(::LatticeGeometry)`](@ref)
"""
abstract type LatticeGeometry end

function Base.show(io::IO, geom::LatticeGeometry)
    print(io, nameof(typeof(geom)), size(geom))
end

"""
    onr(add::AbstractFockAddress, geom::LatticeGeometry)
Returns the occupation number representation of a Fock state address as an `SArray`
with the shape of the lattice geometry `geom`. For composite addresses, a tuple
of `onr`s is returned.
"""
function BitStringAddresses.onr(add, geom::LatticeGeometry)
    return reshape(onr(add), size(geom))
end
function BitStringAddresses.onr(add::CompositeFS, geom::LatticeGeometry)
    return map(fs -> onr(fs, geom), add.components)
end

"""
    neighbour_site(geom::LatticeGeometry, site, i)

Find the `i`-th neighbour of `site` in the geometry. If the move is illegal, return 0.
"""
neighbour_site

"""
    num_neighbours(geom::LatticeGeometry)

Return the number of neighbours each lattice site has in this geometry.

Note that for efficiency reasons, all sites are expected to have the same number of
neighbours. If some of the neighbours are invalid, this is handled by having
[`neighbour_site`](@ref) return 0.
"""
num_neighbours

"""
    PeriodicBoundaries(size...) <: LatticeGeometry

Rectangular lattice with periodic boundary conditions of size `size`.

The dimension of the lattice is controlled by the number of arguments given to its
constructor.

This is the default geometry used by [`HubbardRealSpace`](@ref).

## Example

```
julia> lattice = PeriodicBoundaries(5, 4) # 2D lattice of size 5 × 4
PeriodicBoundaries(5, 4)

julia> num_neighbours(lattice)
4

julia> neighbour_site(lattice, 1, 1)
2

julia> neighbour_site(lattice, 1, 2)
5

julia> neighbour_site(lattice, 1, 3)
6

julia> neighbour_site(lattice, 1, 4)
16
```

## See also

* [`LatticeGeometry`](@ref)
* [`HardwallBoundaries`](@ref)
* [`LadderBoundaries`](@ref)
* [`num_neighbours`](@ref)
* [`neighbour_site`](@ref)
"""
struct PeriodicBoundaries{D} <: LatticeGeometry
    size::NTuple{D,Int}
end

PeriodicBoundaries(args::Vararg{Int}) = PeriodicBoundaries(args)

Base.size(geom::PeriodicBoundaries) = geom.size

num_neighbours(::PeriodicBoundaries{D}) where {D} = 2D

function neighbour_site(geom::PeriodicBoundaries{D}, site, i) where {D}
    # Neighbour indexing
    # i |  x  y  z …
    # 0 | +1  0  0
    # 1 | -1  0  0
    # 2 |  0 +1  0
    # 3 |  0 -1  0
    # 4 |  0  0 +1
    #         ⋮    ⋱
    i -= 1
    cart_indices = CartesianIndices(size(geom))
    cart_index = Tuple(cart_indices[site])
    offset = ntuple(Val(D)) do k
        ifelse(2(k - 1) ≤ i < 2k, ifelse(iseven(i), 1, -1), 0)
    end
    new_index = CartesianIndex(mod1.(cart_index .+ offset, size(geom) .% UInt))
    return LinearIndices(cart_indices)[new_index]
end

"""
    HardwallBoundaries

Rectangular lattice with hard wall boundary conditions of size `size`.
[`neighbour_site()`](@ref) will return 0 for some neighbours of boundary sites.

The dimension of the lattice is controlled by the number of arguments given to its
constructor.

## Example

```
julia> lattice = HardwallBoundaries(5) # 1D lattice of size 5
HardwallBoundaries(5)

julia> neighbour_site(lattice, 1, 1)
2

julia> neighbour_site(lattice, 1, 2)
0

```

## See also

* [`LatticeGeometry`](@ref)
* [`PeriodicBoundaries`](@ref)
* [`LadderBoundaries`](@ref)
"""
struct HardwallBoundaries{D} <: LatticeGeometry
    size::NTuple{D,Int}
end

HardwallBoundaries(args::Vararg{Int}) = HardwallBoundaries(args)

Base.size(geom::HardwallBoundaries) = geom.size

num_neighbours(::HardwallBoundaries{D}) where {D} = 2D

function neighbour_site(geom::HardwallBoundaries{D}, site, i) where {D}
    i -= 1
    cart_indices = CartesianIndices(size(geom))
    cart_index = Tuple(cart_indices[site])
    offset = ntuple(Val(D)) do k
        ifelse(2(k - 1) ≤ i < 2k, ifelse(iseven(i), 1, -1), 0)
    end
    new_index = CartesianIndex(cart_index .+ offset)
    if new_index in cart_indices
        return LinearIndices(cart_indices)[new_index]
    else
        return 0
    end
end

"""
    LadderBoundaries(size...; subgeometry=PeriodicBoundaries) <: LatticeGeometry

Lattice geometry where the first dimension is of size 2 and has hardwall boundary conditions.
Using this geometry is more efficient than using [`HardwallBoundaries`](@ref) with a size of
2, as it does not generate rejected neighbours.

In other dimensions, it behaves like its subgeometry, which can be any
[`LatticeGeometry`](@ref).

## Example

```
julia> lattice = LadderBoundaries(2, 3, 4) # 3D lattice of size 2 × 3 × 4
LadderBoundaries(2, 3, 4)

julia> num_neighbours(lattice)
5

julia> neighbour_site(lattice, 1, 1)
2

julia> neighbour_site(lattice, 1, 2)
3

julia> neighbour_site(lattice, 1, 3)
5

julia> neighbour_site(lattice, 1, 4)
7

julia> neighbour_site(lattice, 1, 5)
19
```

## See also

* [`LatticeGeometry`](@ref)
* [`PeriodicBoundaries`](@ref)
* [`HardwallBoundaries`](@ref)
"""
struct LadderBoundaries{G<:LatticeGeometry} <: LatticeGeometry
    subgeometry::G
end

function LadderBoundaries(arg::Int, args::Vararg{Int}; subgeometry=PeriodicBoundaries)
    if arg ≠ 2
        error("First dimension must be of size 2")
    end
    return LadderBoundaries(subgeometry(args...))
end

Base.size(geom::LadderBoundaries) = (2, size(geom.subgeometry)...)

num_neighbours(geom::LadderBoundaries) = num_neighbours(geom.subgeometry) + 1

function neighbour_site(geom::LadderBoundaries, site, i)
    if i == 1
        # Make odd site even, or make even site odd.
        return site + site % 2 - (site + 1) % 2
    else
        i -= 1
        subneighbour = neighbour_site(geom.subgeometry, cld(site, 2), i)
        return ifelse(iszero(subneighbour), 0, 2 * subneighbour - site % 2)
    end
end

ext: jl


text:
## Filter the MCMC results for DBSCAN on the ϵ+Xdisp metric (see NBs)
##

using DataFrames , CSV , Statistics
using Query
using PyPlot
using TSne
using Printf

import CSV

rootdir = ENV["GAIA_ROOT"]
wdir    = "$rootdir/e2e_products"
plotdir = "$rootdir/e2e_products"

cd(wdir)

## Input dataframe
mcfull= CSV.read("votlist-mcmc_full-TEST.csv", delim = ";")
scfull= CSV.read("votlist-SCproperties_full-TEST.csv", delim = ";")
dfcomb= join(mcfull, scfull, on = :votname, makeunique=true)

## Building the metric
x = @from i in dfcomb begin
    @let votname=i.votname
    @let q1 = (i.epsm / i.mclm) * (i.xdisp / i.distance)

    @select {votname=votname,Q1=q1}
    @collect DataFrame
end

println(x)

## Plot the metric histogram
fig = figure(figsize=(7,5))
nbins=50
h = plt.hist(x.Q1,nbins,range = [0,0.15], color = "g",  label = "dbscan+weight density" , alpha=0.6 , density =true)
plt.xlabel("ϵ+Xdisp Metric")
plt.savefig("dbscanMetric.png")

## filtering
# Threshold maximum on Q1 to filter the OC. Check it in the plot dbscanMetric
#
THRESHOLD_Q1= 0.02

final = @from i in dfcomb begin
    @where (i.xdisp / i.distance)*(i.epsm / i.mclm) < THRESHOLD_Q1
    @select i
    @collect DataFrame
end

println("# N selected")
println(size(final))
fileout= @sprintf("votlist.finalSample_metric%3.3f.csv", THRESHOLD_Q1)
CSV.write(fileout, final, delim = ';')

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text: Dict{String,Any}("milliseconds" => Dict{String,Any}("value" => "1977-12-21T03:32:00.555+00:00","type" => "datetime"),"bestdayever" => Dict{String,Any}("value" => "1987-07-05T17:45:00Z","type" => "datetime"),"numoffset" => Dict{String,Any}("value" => "1977-06-28T12:32:00Z","type" => "datetime"))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using Gadfly, Distributions

set_default_plot_size(6inch, 3inch)

plot(x=rand(Rayleigh(2),1000), y=rand(Rayleigh(2),1000),
    Geom.density2d(levels = x->maximum(x)*0.5.^collect(1:2:8)), Geom.point,
    Theme(key_position=:none),
    Scale.color_continuous(colormap=x->colorant"red"))


text:
export run_model

function run_model(
    data::AbstractHLAData;
    keep_all_parameters::Bool = false,
    save_file::String = "",
    stan_kwargs...
) where T

    if save_file != ""
        save(save_file, Dict("result" => "test"))
        rm(save_file)
    end

    result = Escape.run(
        Escape.HLAModel{4}(), data; keep_all_parameters = keep_all_parameters, stan_kwargs...
    )

    if save_file != ""
        save(save_file, Dict("result" => result))
    end

    return result
end

function run_model(
    model::HLAModel{T}, data::AbstractHLAData;
    keep_all_parameters::Bool = false,
    save_file::String = "",
    stan_kwargs...
) where T

    if save_file != ""
        save(save_file, Dict("result" => "test"))
        rm(save_file)
    end

    result = Escape.run(
        model, data;
        keep_all_parameters = keep_all_parameters, stan_kwargs...
    )

    if save_file != ""
        save(save_file, Dict("result" => result))
    end

    return result
end

function Escape.run(
    model::HLAModel{T}, data::AbstractHLAData;
    keep_all_parameters::Bool = false,
    stan_kwargs...
) where T

    !ismissing(stan_input(data)) || error("HLAData does not contain a stan_input.")

    input = stan_input(data)
    replacements = Escape.replacements(data)
    alleles = sort(unique_alleles(Escape.hla_types(data),
        allele_depth = allele_depth(data)))

    sf = model_1() do exe
        model_path = exe[begin:end-1] * string(T)

        StanInterface.stan(
            model_path, input;
            stan_args = "adapt delta=0.85 algorithm=hmc engine=nuts max_depth=10",
            refresh = 1, iter = 1000, warmup = 300, stan_kwargs...
        )
    end

    !keep_all_parameters && reduce_size!(sf)

    return HLAModelResult(model, data, sf, replacements, alleles)
end

function reduce_size!(sf::StanInterface.Stanfit)
    for d in sf.result
        for key in keys(d)
            if !any(startswith.(key, ["lp__", "b_phy", "b0_hla", "beta_hla"]))
                delete!(d, key)
            end
        end
    end

    return sf
end


text:
# Submodule that constructs the system of Shearlets in 2D of size rows x cols

######################################################################
## ype of Filterconfigsuration to check the sizes
struct Filterswedgebandlow
  wedge::AbstractArray
  bandpass::AbstractArray
  lowpass::AbstractArray
end

#######################################################################
# Function that generates the whole Shearlet System filters (wedge, bandpass and lowpass) of size rows x cols
"""
  ...
  getwedgebandpasslowpassfilters2D(rows,cols,directionalFilter = filt_gen("directional_shearlet"), scalingFilter = filt_gen("scaling_shearlet"),gpu = false)
  generates the Wedge, Bandpass and LowspassFilter of size rows x cols
  ...
  """

function getwedgebandpasslowpassfilters2D(rows::Int, cols::Int, shearLevels,
                                          directionalFilter =
                                          filt_gen("directional_shearlet"),
                                          scalingFilter =
                                          filt_gen("scaling_shearlet"),
                                          waveletFilter =
                                          mirror(scalingFilter),
                                          scalingFilter2 = scalingFilter, gpu =
                                          false)

  FFTW.set_num_threads(Sys.CPU_THREADS)
  # Make shearLevels integer
  shearLevels = map(Int,shearLevels);
  # The number of scales
  NScales = length(shearLevels);

  # Initialize the bandpass and wedge filter
  # Bandpass filters partion the frequency plane into different scales
  if gpu
    bandpass = zeros(Float32,rows,cols,NScales);
    bandpass = AFArray(bandpass+im*zeros(bandpass));
  else
    bandpass = zeros(Float64,rows,cols,NScales);
    bandpass = bandpass+im*zeros(size(bandpass));
  end

  # Wedge filters partition the frequency plane into different directions
  wedge = Array{Any}(undef, 1, round(Int64,maximum((shearLevels) .+ 1)));

  # normalize the directional filter directional filter
  if gpu
    directionalFilter = AFArray(convert(Array{Float32},directionalFilter/sum(abs.(directionalFilter[:]))));
  else
    directionalFilter = directionalFilter/sum(abs.(directionalFilter[:]));
  end

  # Compute 1D high and lowpass filters at different scales

  # filterHigh{NScales} = g, and filterHigh{1} = g_J (analogous with filterLow)
  filterHigh = Array{Any}(undef, 1,NScales);
  filterLow = Array{Any}(undef, 1,NScales);
  filterLow2 = Array{Any}(undef, 1,round(Int64,maximum(shearLevels)+1));

  # Initialize wavelet highpass and lowpass filters
  if gpu
    filterHigh[size(filterHigh,2)] = AFArray(convert(Array{Float32},waveletFilter));
    filterLow[size(filterLow,2)] = AFArray(convert(Array{Float32},scalingFilter));
    filterLow2[size(filterLow2,2)] = AFArray(convert(Array{Float32},scalingFilter2));

    # Lets compute the filters in the other scales
    for j = (size(filterHigh)[2]-1):-1:1
      filterLow[j] = conv2(filterLow[size(filterLow,2)],upsample(filterLow[j+1],1,1,1))
      filterHigh[j] = conv(filterLow[size(filterLow,2)],upsample(filterHigh[j+1],1,1,1))
    end

    # Lets compute the filters in the other scales
    for j=(size(filterLow2)[2]-1):-1:1
      filterLow2[j] = conv(filterLow2[size(filterLow2,2)],upsample(filterLow2[j+1],1,1,1))
    end
    # Construct the bandpassfilter
    # Need to convert first to complex array since
    for j = 1:size(filterHigh,2)
      bandpass[:,:,j] = -fftshift(fft(ifftshift(padarray(filterHigh[j],[rows,cols],1))));
    end
  else
    filterHigh[size(filterHigh,2)] = waveletFilter;
    filterLow[size(filterLow,2)] = scalingFilter;
    filterLow2[size(filterLow2,2)] = scalingFilter2;

    # Lets compute the filters in the other scales
    for j = (size(filterHigh)[2]-1):-1:1
      filterLow[j] = conv(filterLow[size(filterLow,2)],upsample(filterLow[j+1],1,1))
      filterHigh[j] = conv(filterLow[size(filterLow,2)],upsample(filterHigh[j+1],1,1))
    end

    # Lets compute the filters in the other scales
    for j=(size(filterLow2)[2]-1):-1:1
      filterLow2[j] = conv(filterLow2[size(filterLow2,2)],upsample(filterLow2[j+1],1,1))
    end

    # Construct the bandpassfilter
    # Need to convert first to complex array since
    for j = 1:size(filterHigh,2)
      bandpass[:,:,j] = fftshift(fft(ifftshift(padarray(filterHigh[j],[rows,cols]))));
    end
  end

  ## construct wedge filters for achieving directional selectivity.
  # as the entries in the shearLevels array describe the number of differently
  # sheared atoms on a certain scale, a different set of wedge
  # filters has to be constructed for each value in shearLevels.
  if gpu
    for shearLevel = unique(shearLevels)
      #preallocate a total of floor(2^(shearLevel+1)+1) wedge filters, where
      #floor(2^(shearLevel+1)+1) is the number of different directions of
      #shearlet atoms associated with the horizontal (resp. vertical)
      #frequency cones.
      wedge[shearLevel+1] = AFArray(zeros(Complex{Float32},rows,cols,floor(2^(shearLevel+1)+1)));
      #upsample directional filter in y-direction. by upsampling the directional
      #filter in the time domain, we construct repeating wedges in the
      #frequency domain ( compare abs.(fftshift(fft2(ifftshift(directionalFilterUpsampled)))) and
      #abs.(fftshift(fft2(ifftshift(directionalFilter)))) ).
      directionalFilterUpsampled = upsample(directionalFilter,1,2^(shearLevel+1)-1,1);

      #remove high frequencies along the y-direction in the frequency domain.
      #by convolving the upsampled directional filter with a lowpass filter in y-direction, we remove all
      #but the central wedge in the frequency domain.
      wedgeHelp = conv2(directionalFilterUpsampled, filterLow2[size(filterLow2,2)-shearLevel][:,:])
      wedgeHelp = padarray(wedgeHelp,[rows,cols],1);
      #please note that wedgeHelp now corresponds to
      #conv(p_j,h_(J-j*alpha_j/2)') in the language of the paper. to see
      #this, consider the definition of p_j on page 14, the definition of w_j
      #on the same page an the definition of the digital sheralet filter on
      #page 15. furthermore, the g_j part of the 2D wavelet filter w_j is
      #invariant to shearings, hence it suffices to apply the digital shear
      #operator to wedgeHelp.;

      ## application of the digital shear operator (compare equation (22))

      #upsample wedge filter in x-direction. this operation corresponds to
      #the upsampling in equation (21) on page 15.
      wedgeUpsampled = upsample(wedgeHelp,2,2^shearLevel-1,1);
      #convolve wedge filter with lowpass filter, again following equation
      #(21) on page 14.
      lowpassHelp = padarray(filterLow2[size(filterLow2,2)-max(shearLevel-1,0)],size(wedgeUpsampled),1);
      if shearLevel >= 1
        wedgeUpsampled = fftshift(ifft(ifftshift(fftshift(fft(ifftshift(lowpassHelp))).*fftshift(fft(ifftshift(wedgeUpsampled))))));
      end
      lowpassHelpFlip = fliplr(lowpassHelp,1);
      #traverse all directions of the upper part of the left horizontal
      #frequency cone
      for k = -2^shearLevel:2^shearLevel
        #resample wedgeUpsampled as given in equation (22) on page 15.
        wedgeUpsampledSheared = dshear(wedgeUpsampled,k,2,1);
        #convolve again with flipped lowpass filter, as required by equation (22) on
        #page 15.
        if shearLevel >= 1
          wedgeUpsampledSheared = fftshift(ifft(ifftshift(fftshift(fft(ifftshift(lowpassHelpFlip))).*fftshift(fft(ifftshift(wedgeUpsampledSheared))))));
        end
        #obtain downsampled and renormalized and sheared wedge filter in the
        #frequency domain, according to equation (22) on page 15
        wedge[shearLevel+1][:,:,fix(2^shearLevel)+1-k] = fftshift(fft(ifftshift(2^shearLevel*wedgeUpsampledSheared[:,1:2^shearLevel:(2^shearLevel*cols-1)])));
      end
    end
  else
    for shearLevel = unique(shearLevels)
      #preallocate a total of floor(2^(shearLevel+1)+1) wedge filters, where
      #floor(2^(shearLevel+1)+1) is the number of different directions of
      #shearlet atoms associated with the horizontal (resp. vertical)
      #frequency cones.
      wedge[shearLevel+1] = zeros(rows,cols,floor(2^(shearLevel+1)+1))+zeros(rows,cols,floor(2^(shearLevel+1)+1))*im;

      #upsample directional filter in y-direction. by upsampling the directional
      #filter in the time domain, we construct repeating wedges in the
      #frequency domain ( compare abs.(fftshift(fft2(ifftshift(directionalFilterUpsampled)))) and
      #abs.(fftshift(fft2(ifftshift(directionalFilter)))) ).
      directionalFilterUpsampled = upsample(directionalFilter,1,2^(shearLevel+1)-1);

      #remove high frequencies along the y-direction in the frequency domain.
      #by convolving the upsampled directional filter with a lowpass filter in y-direction, we remove all
      #but the central wedge in the frequency domain.

      wedgeHelp = conv(directionalFilterUpsampled, filterLow2[size(filterLow2,2)-shearLevel][:,:])
      wedgeHelp = padarray(wedgeHelp,[rows,cols]);

      #please note that wedgeHelp now corresponds to
      #conv(p_j,h_(J-j*alpha_j/2)') in the language of the paper. to see
      #this, consider the definition of p_j on page 14, the definition of w_j
      #on the same page an the definition of the digital sheralet filter on
      #page 15. furthermore, the g_j part of the 2D wavelet filter w_j is
      #invariant to shearings, hence it suffices to apply the digital shear
      #operator to wedgeHelp.;

      ## application of the digital shear operator (compare equation (22))

      #upsample wedge filter in x-direction. this operation corresponds to
      #the upsampling in equation (21) on page 15.
      wedgeUpsampled = upsample(wedgeHelp,2,2^shearLevel-1);

      #convolve wedge filter with lowpass filter, again following equation
      #(21) on page 14.
      lowpassHelp = padarray(filterLow2[size(filterLow2,2)-max(shearLevel-1,0)],size(wedgeUpsampled));
      if shearLevel >= 1
        wedgeUpsampled = fftshift(ifft(ifftshift(fftshift(fft(ifftshift(lowpassHelp))).*fftshift(fft(ifftshift(wedgeUpsampled))))));
      end
      lowpassHelpFlip = fliplr(lowpassHelp);
      #traverse all directions of the upper part of the left horizontal
      #frequency cone
      for k = -2^shearLevel:2^shearLevel
        #resample wedgeUpsampled as given in equation (22) on page 15.
        wedgeUpsampledSheared = dshear(wedgeUpsampled,k,2);
        #convolve again with flipped lowpass filter, as required by equation (22) on
        #page 15.
        if shearLevel >= 1
          wedgeUpsampledSheared = fftshift(ifft(ifftshift(fftshift(fft(ifftshift(lowpassHelpFlip))).*fftshift(fft(ifftshift(wedgeUpsampledSheared))))));
        end
        #obtain downsampled and renormalized and sheared wedge filter in the
        #frequency domain, according to equation (22) on page 15
        wedge[shearLevel+1][:,:,fix(2^shearLevel)+1-k] = fftshift(fft(ifftshift(2^shearLevel*wedgeUpsampledSheared[:,1:2^shearLevel:(2^shearLevel*cols-1)])));
      end
    end
end

## compute low pass filter of shearlet system
if gpu
  lowpass = fftshift(fft(ifftshift(padarray(transpose(filterLow[1]')*filterLow[1]',[rows,cols],1))));
else
  lowpass = fftshift(fft(ifftshift(padarray(transpose(filterLow[1]')*filterLow[1]',[rows,cols]))));
end

# Generate the final array
return Filterswedgebandlow(wedge,bandpass,lowpass)
end #getwedgebandpasslowpassfilters2D

##############################################################
# Create a type for the Preparedfilters
struct Preparedfilters
  size::Tuple{Int,Int}
  shearLevels::Array{Int,1}
  cone1
  cone2
end

#######################################################################
# Function that prepare the filters
"""
  ...
  SlprepareFilters2D(rows, cols, nScales, shearLevels = ceil.((1:nScales)/2),
      directionalFilter = filt_gen("directional_shearlet"),
      scalingFilter = filt_gen("scaling_shearlet"),
      waveletFilter = mirror(scalingFilter),
      scalingFilter2 = scalingFilter, gpu = false) function that prepare the filters to generate
		   the shearlet system
  ...
  """
function preparefilters2D(rows, cols, nScales, shearLevels = ceil.((1:nScales)/2),
                          directionalFilter = filt_gen("directional_shearlet"),
                          scalingFilter = filt_gen("scaling_shearlet"),
                          waveletFilter = mirror(scalingFilter),
                          scalingFilter2 = scalingFilter, gpu = false)

  #Make sure the shearLevles are integer
  shearLevels = map(Int,shearLevels)
  # check filter sizes
  filters = checkfiltersizes(rows,cols,shearLevels,directionalFilter,scalingFilter,waveletFilter,scalingFilter2);

  # Save the new filters
  directionalFilter = filters.directionalFilter;
  scalingFilter = filters.scalingFilter;
  waveletFilter = filters.waveletFilter;
  scalingFilter2 = filters.scalingFilter2;

  # Define the cones
  cone1 = getwedgebandpasslowpassfilters2D(rows,cols,shearLevels,directionalFilter,scalingFilter,waveletFilter,scalingFilter2,gpu);
  if rows == cols
    cone2 = cone1;
  else
    cone2 = getwedgebandpasslowpassfilters2D(cols,rows,shearLevels,directionalFilter,scalingFilter,waveletFilter,scalingFilter2,gpu);
  end
  return Preparedfilters((rows,cols),shearLevels,cone1,cone2)
end # preparefilters2D

#######################################################################
# Function compute a index set describing a 2D shearlet system
"""
  ...
  getshearletidxs2D(shearLevels, full=0, restriction_type = [],restriction_value = []) Compute a index set describing a 2D shearletsystem, with restriction_type in
  ["cones","scales","shearings"], where shearLevels is a 1D array spcifying the level of shearing on each scale, and full is a boolean that orders to generate the full shearlet indices.
  ...
  """
function getshearletidxs2D(shearLevels, full=0, restriction_type = [],restriction_value = [])
	shearLevels = map(Int,shearLevels)
  # Set default values
  shearletIdxs = [];
  includeLowpass = 1;
  scales = 1:length(shearLevels);
  shearings = -2^(maximum(shearLevels)):2^(maximum(shearLevels));
  cones = 1:2;

  for i = 1:length(restriction_type)
    includeLowpass = 0;
    if restriction_type[i] == "scales"
      scales = restriction_value[i];
    elseif restriction_type[i] == "shearings"
      shearings = restriction_value[i];
    elseif restriction_type[i] == "cones"
      cones = restriction_value[i];
    end
  end

  for cone = intersect(1:2,cones)
    for scale = intersect(1:length(shearLevels),scales)
      for shearing = intersect(-2^shearLevels[scale]:2^shearLevels[scale],shearings)
        if convert(Bool,full) || cone == 1 || abs.(shearing)<2^(shearLevels[scale])
          shearletIdxs=[shearletIdxs;[cone,scale,shearing]];
        end
      end
    end
  end

  shearletIdxs = transpose(reshape(shearletIdxs,3,Int(length(shearletIdxs)/3)))
  if convert(Bool,includeLowpass) || convert(Bool,sum(0 .== scales)) || convert(Bool,sum(0 .== cones))
    shearletIdxs = [shearletIdxs;[0 0 0]];
  end
  return Int.(shearletIdxs)
end # getshearletidxs2D

####################################################################
# Type of shearletsystem in 2D
struct Shearletsystem2D{T<:Number, CT<:Union{Complex{T}, T}} # CT tells us
    # whether we should be using a complex or real fft
    shearlets::Array{Complex{T}, 3}
    size::Tuple{Int, Int}
    shearLevels::Array{Int, 1}
    full::Bool
    nShearlets::Int
    shearletIdxs::Array{Int, 2}
    dualFrameWeights::Array{T, 2}
    RMS::Array{T, 2}
    gpu::Bool
    padded::Bool
    padBy::Tuple{Int, Int}
end

function Base.show(io::IO, l::Shearletsystem2D{T,CT}) where {T,CT}
    print(io, "Shearletsystem2D{$(T),$(CT)}(input = $(l.size), levels = $(l.shearLevels),"*
          " nshearlets = $(l.nShearlets),gpu=$(l.gpu)"*(l.padded ?
                                                        ", padBy = $(l.padBy))" : ""))
end

#######################################################################
# helper methods to reduce code redundandancy

"""
    shearlets, dualFrameWeights = padShearlets(shearlets, dualFrameWeights, typeBecomes)

Add padding in the space domain and convert back. Also, if the type being
transformed is real, only store half the coefficients, as they're otherwise
redundant.
"""
function padShearlets(shearlets, dualFrameWeights, typeBecomes, padBy, upperFrameBound)
    shearlets = real.(fftshift(ifft(ifftshift(shearlets, (1,2)), (1, 2)),
                               (1, 2)))
    shearlets = cat([pad(shearlets[:,:,j], padBy) for
                     j=1:size(shearlets,3)]...; dims = 3)
    shearlets = ifftshift(shearlets, (1,2))
    if typeBecomes <: Real
        P = plan_rfft(zeros(typeBecomes, size(shearlets)[1:2]), (1,2))
        newSize = (div(size(shearlets,1), 2)+1, size(shearlets)[2:3]...)
    else
        P = plan_fft(zeros(typeBecomes, size(shearlets)[1:2]), (1,2))
        newSize = size(shearlets)
    end
    newShears = zeros(Complex{eltype(shearlets)}, newSize)
    totalMass = sum(abs.(shearlets[:,:,1]))
    # rowMean = sum((1:size(shearlets,1)) .* abs.(shearlets[:,:,1]), dims=(1,2))/totalMass
    # colMean = sum((1:size(shearlets,2))' .* abs.(shearlets[:,:,1]), dims=(1,2))/totalMass
    # println("newish location is $((rowMean, colMean))) out of $(size(shearlets)[1:2])")

    for j=1:size(shearlets, 3)
        newShears[:, :, j] = P * shearlets[:, :, j]
    end
    dualFrameWeights = sum(real.(abs.(newShears).^2), dims=3)[:,:]
    if upperFrameBound > 0
        totalMass = norm(dualFrameWeights, Inf)
        normalize = typeBecomes(upperFrameBound)/totalMass
        newShears = newShears .* normalize
        dualFrameWeights = sum(real.(abs.(newShears).^2), dims=3)[:,:]
    end
    return (newShears, dualFrameWeights)
end

"""
    shearlets, dualFrameWeights, RMS, rows, cols, nShearlets =
        generateShearlets(shearletIdxs, Preparedfilters,typeBecomes)

the core code for generating shearlets
"""
function generateShearlets(shearletIdxs, Preparedfilters, typeBecomes, gpu)
    rows = Preparedfilters.size[1];
    cols = Preparedfilters.size[2];
    nShearlets = size(shearletIdxs,1);
    if gpu
        shearlets = AFArray(zeros(Complex{Float32},rows,cols,nShearlets));
    else
        shearlets = zeros(typeBecomes, rows, cols, nShearlets)+im .*
            zeros(rows, cols, nShearlets);
    end
    # Compute shearlets
    for j = 1:nShearlets
        cone = shearletIdxs[j,1];
        scale = shearletIdxs[j,2];
        shearing = shearletIdxs[j,3];

        if cone == 0
            shearlets[:,:,j] = Preparedfilters.cone1.lowpass;
        elseif cone == 1
            #here, the fft of the digital shearlet filters described in
            #equation (23) on page 15 of "ShearLab 3D: Faithful Digital
            #Shearlet Transforms based on Compactly Supported Shearlets" is computed.
            #for more details on the construction of the wedge and bandpass
            #filters, please refer to getwedgebandpasslowpassfilters2D.
            shearlets[:,:,j] =
                Preparedfilters.cone1.wedge[Preparedfilters.shearLevels[scale]+1][:,:,-shearing+2^Preparedfilters.shearLevels[scale]+1].*conj(Preparedfilters.cone1.bandpass[:,:,scale]);
        else
            shearlets[:,:,j] =
                permutedims(Preparedfilters.cone2.wedge[Preparedfilters.shearLevels[scale]+1][:,:,shearing+2^Preparedfilters.shearLevels[scale]+1].*conj(Preparedfilters.cone2.bandpass[:,:,scale]),[2,1]);
        end
    end

	RMS =  transpose(sum(reshape(sum(real.(abs.(shearlets)).^2, dims=1), size(shearlets, 2), size(shearlets, 3)), dims=1));
    if gpu
        RMS = (sqrt.(RMS)/convert(T, sqrt.(rows*cols)));
        dualFrameWeights = sum(real(abs.(shearlets)).^2,dims=3);
    else
        RMS = (sqrt.(RMS)/sqrt.(rows*cols));
        dualFrameWeights = dropdims(sum(abs.(shearlets).^2,dims=3),dims=3); # need to stream to host before they fix abs of comples AFArray
    end
    return (shearlets, dualFrameWeights, RMS, rows, cols, nShearlets)
end

#######################################################################
# Function that generates the desired shearlet system
"""
    ...
    getshearletsystem2D(rows,cols,nScales,shearLevels=ceil.((1:nScales)/2),full = 0,directionalFilter = filt_gen("directional_shearlet"),quadratureMirrorFilter= filt_gen("scaling_shearlet"), gpu = false)

generates the desired shearlet system
  """
function getshearletsystem2D(rows, cols, nScales,
                             shearLevels=ceil.((1:nScales)/2),
                             full= false,
                             directionalFilter = filt_gen("directional_shearlet"),
                             quadratureMirrorFilter=
                             filt_gen("scaling_shearlet"), gpu=false;
                             typeBecomes=Float64, padded=true, tolerance =
                             1e-10, upperFrameBound = -1)

    # Set default value generates the desired shearlet systems
    shearLevels = Int.(shearLevels)

    #Generate prepared Filters and indices
    Preparedfilters = preparefilters2D(rows, cols, nScales, shearLevels,
                                       directionalFilter,
                                       quadratureMirrorFilter,
                                       mirror(quadratureMirrorFilter),
                                       quadratureMirrorFilter,  gpu);
    shearletIdxs = getshearletidxs2D(shearLevels, full);

    # Generate shearlets, RMS(rootmeansquare), dualFrameWeights
    shearlets, dualFrameWeights, RMS, rows, cols, nShearlets =
        generateShearlets(shearletIdxs, Preparedfilters, typeBecomes, gpu)

    # adjust sizes if we're padding, and/or store using half the coefficients
    # if the data is real
    if padded
        # this is probably inefficient, but it works
        padBy = getPadBy(shearlets, tolerance = tolerance)
        shearlets, dualFrameWeights = padShearlets(shearlets, dualFrameWeights,
                                                   typeBecomes, padBy, upperFrameBound)
    else
        padBy = (0,0)
    end

    if typeBecomes <: Real
        T = typeBecomes
        CT = typeBecomes
    else
        T = real(typeBecomes)
        CT = typeBecomes
    end
    # we were doing this a lot. Let's do it once
    shearlets = conj(shearlets)

    #return the system
    return Shearletsystem2D{T, CT}(shearlets, Preparedfilters.size,
                            Preparedfilters.shearLevels, full,
                            size(shearletIdxs, 1), shearletIdxs,
                            dualFrameWeights, RMS, gpu, padded, padBy)
end #getshearletsystem2D

# type for individual shearlets2D
struct Shearlets2D{T<:Number, CT <: Union{Complex{T}, Nothing}}
    shearlets::Array{Complex{T}, 3}
    RMS::Array{T, 2}
    dualFrameWeights::Array{T, 2}
    gpu::Bool
    padded::Bool
    padBy::Tuple{Int, Int}
end

#######################################################################
# Function that generates the desired shearlets
"""
   getshearlets2D(PreparedFilters,shearletIdxs,gpu = false)

generates the 2D shearlets in the frequency domain
"""
function getshearlets2D(Preparedfilters, shearletIdxs =
                        getshearletidxs2D(Preparedfilters.shearLevels), gpu =
                        false; typeBecomes=Float64, padded=true,
                        upperFrameBound = -1)
  # Generate shearlets, RMS(rootmeansquare), dualFrameWeights
    shearlets, dualFrameWeights, RMS, rows, cols, nShearlets =
        generateShearlets(shearletIdxs, Preparedfilters, typeBecomes, gpu)

    # adjust sizes if we're padding, and/or store using half the coefficients
    # if the data is real
    padBy = getPadBy(shearlets, tolerance = tolerance)
    if padded || typeBecomes <: Real
        # this is probably inefficient, but it works
        shearlets, dualFrameWeights = padShearlets(shearlets, dualFrameWeights,
                                                   typeBecomes, padBy, upperFrameBound)
    end

    if typeBecomes <: Real
        T = typeBecomes
        CT = typeBecomes
    else
        T = real(typeBecomes)
        CT = typeBecomes
    end

	return Shearlets2D{T, CT}(shearlets, RMS, dualFrameWeights, gpu, padded, padBy)
end #getshearlets2D


text:
LCR6xxx_Misc_Dict = Dict()
# gridmode
LCR6xxx_Misc_Dict["single"] = 1
LCR6xxx_Misc_Dict["double"] = 2
LCR6xxx_Misc_Dict["quad"]   = 3

# Reset
reset_instr(obj::T) where {T<:LCR6xxx} = write(obj,"*RST")
# IDN
get_idn(obj::T) where {T<:LCR6xxx} = query(obj,"*IDN?")
# Default setup
default_setup(obj::T) where {T<:LCR6xxx} = write(obj,"VBS app.SetToDefaultSetup")
# Screen Gridmode
function set_gridmode(obj::T;gridmode::String="single") where {T<:LCR6xxx}
	gridmode = LCR6xxx_Misc_Dict[gridmode];
	write(obj,"VBS app.Display.GridMode=$gridmode")
end
# Clear sweeps
clear_sweeps(obj::T) where {T<:LCR6xxx} = write(obj,"VBS app.ClearSweeps")


text:
const PREPROCESSING_EXTRAKEY = :input
const FAST = :fast
const PRECISE = :precise
const INTERCEPT_DEFAULT = true
const METHOD_DEFAULT = FAST
const TIME_DEFAULT = nothing
const PANEL_DEFAULT = nothing
const FIXED_EFFECT_DEFAULT = false
const INTERACTION_DEFAULT = nothing
const REMOVEOUTLIERS_DEFAULT = false
const REMOVEMISSINGS_DEFAULT = false
const SEASONALADJUSTMENT_DEFAULT = nothing
const DEFAULT_DATANAME = "dataname"


text:
## Metsa-Oy Forest Products Supply Chain
## Edward J. Xu <edxu96@outlook.com>
## Jan 24th, 2020

using JuMP, Gurobi

include("func.jl")
include("para.jl")
include("optim_1.jl")
include("optim_2.jl")
include("optim_3.jl")
include("result_2.jl")
include("result_3.jl")

function main()
  alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j, c_j,
    r_j, a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik, gamma_jk,
    delta_jk, gamma_paper_k, delta_paper_k, p_fuel, sigma, omega_i, omega_j,
    omega_paper, nu_saw, nu_plywood, nu_j, nu_paper, o_saw, o_plywood, o_j,
    o_paper, i_t1, i_t2, j_t2, pi_n, rho_nm = get_para()

  optim_1(alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j, c_j, r_j,
    a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik, gamma_jk,
    delta_jk, gamma_paper_k, delta_paper_k, p_fuel)

  optim_2(alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j, c_j,
    r_j, a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik, gamma_jk,
    delta_jk, gamma_paper_k, delta_paper_k, p_fuel, sigma, omega_i, omega_j,
    omega_paper, nu_saw, nu_plywood, nu_j, nu_paper, o_saw, o_plywood, o_j,
    o_paper, i_t1, i_t2, j_t2)

  optim_3(alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j,
    c_j, r_j, a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik,
    gamma_jk, delta_jk, gamma_paper_k, delta_paper_k, p_fuel, sigma, omega_i,
    omega_j, omega_paper, nu_saw, nu_plywood, nu_j, nu_paper, o_saw, o_plywood,
    o_j, o_paper, i_t1, i_t2, j_t2, pi_n, rho_nm)
end

main()


text:
# <editor-fold desc= top-level and general read-in functions for sets and parameter"

# XXX read-in all set files
function readSets!(files_dic::Dict{String,Array{String,1}},anyM::anyModel)

	# creates relevant sets, manually adds :mode and assigns short names
	set_arr = append!(map(x ->  Symbol(x[findfirst("set_",x)[1]+4:end-4]),files_dic["set"]), [:mode, :id])
	setLngShrt_dic = Dict(:timestep => :Ts, :region => :R, :carrier => :C, :technology => :Te, :mode => :M, :id => :id)

	setData_dic = Dict{Symbol,DataFrame}()
	anyM.sets = Dict{Symbol,Tree}()

	# loop over sets read-in data and create tree objects
	for setFile in files_dic["set"]
		setLong_sym = getindex(setFile[findfirst("set_",setFile)[1]+4:end-4] |> (y -> filter(x -> occursin(string(x),y),collectKeys(keys(setLngShrt_dic)))),1)
		setShort_sym = get!(setLngShrt_dic,setLong_sym,setLong_sym)
		if setShort_sym in keys(anyM.sets)
			push!(anyM.report,(3,"set read-in",string(setLong_sym),"multiple input files provided for set"))
		end
		setData_dic[setShort_sym] = convertReadIn(CSV.read(setFile;delim = anyM.options.csvDelim[1]),setFile,set_arr,setLngShrt_dic,anyM.report,anyM.lock)
		anyM.sets[setShort_sym] = createTree(setData_dic[setShort_sym],setLong_sym,anyM.report)
	    produceMessage(anyM.options,anyM.report, 3," - Read-in set file: ",setFile)
	end

	# manually adds mode set, creates a top node for all modes to allow for aggregation later
	if :mode in namesSym(setData_dic[:Te])
		modes_df = DataFrame(mode_1 = vcat(map(y -> split(replace(y," " => ""),";"),filter(x -> x != "",setData_dic[:Te][!,:mode]))...))
	else
		modes_df = DataFrame(mode_1 = String[])
	end
	anyM.sets[:M] = createTree(modes_df,:mode,anyM.report)

	# reports, if a required set was not defined or if non-unique carrier names were defined
	for set in filter(x -> !(x in (:mode, :id)), collectKeys(keys(setLngShrt_dic)))
		if !(setLngShrt_dic[set] in keys(anyM.sets))
			push!(anyM.report,(3,"set read-in",string(set),"no file provided to define set"))
		elseif set == :carrier || set == :technology
			# reports error if carrier names are non-unique
			strSet_arr = getfield.(values(anyM.sets[setLngShrt_dic[set]].nodes),:val)
			if length(strSet_arr) != length(unique(strSet_arr))
				push!(anyM.report,(3,"set read-in",set,"non-unique $set names provided"))
			end
		end
	end

	produceMessage(anyM.options,anyM.report, 2," - Read-in all set files")
	return setData_dic
end

# XXX read-in all parameter files, create model parts and assign parameter
function readParameters!(files_dic::Dict{String,Array{String,1}},setData_dic::Dict{Symbol,DataFrame},anyM::anyModel)
	# XXX read-in parameters (do not add to parts yet)
	paraTemp_dic = Dict{String, Dict{Symbol, DataFrame}}()

	# creates relevant sets, manually adds :mode and assigns short names
	set_arr = append!(map(x ->  Symbol(x[findfirst("set_",x)[1]+4:end-4]),files_dic["set"]), [:mode, :id])
	setLngShrt_dic = Dict(:timestep => :Ts, :region => :R, :carrier => :C, :technology => :Te, :mode => :M, :id => :id)

	# read-in parameter files and convert their content
	@threads for parFile in files_dic["par"]
		parData_df = convertReadIn(CSV.read(parFile;delim = anyM.options.csvDelim[1]),parFile,set_arr,setLngShrt_dic,anyM.report,anyM.lock,anyM.sets)
		if isempty(parData_df) || any(getindex.(anyM.report,1) .== 3) continue end
		para_obj =  writeParameter(parData_df, anyM.sets, setLngShrt_dic, parFile, anyM.report, anyM.lock)
		lock(anyM.lock)
		paraTemp_dic[parFile] = para_obj
		unlock(anyM.lock)
		produceMessage(anyM.options,anyM.report, 3," - Read-in parameter file: ",parFile)
	end
	produceMessage(anyM.options,anyM.report, 2," - Read-in all parameter files")
	return paraTemp_dic
end

# XXX read inputs folders for all 'csv' or 'jl' files starting with 'set', 'par', 'var' and 'eqn'
function readInputFolder(inputFolders::Array{String,1},files_dic::Dict{String,Array{String,1}} = Dict(b => String[] for b in ("set","par","var","eqn")))
	# loops over main folders provides in constructor
    for folder in inputFolders
	    for file in readdir(folder)
	        if occursin(".",file) fileType_str = file[findfirst(".",file)[1]+1:end] else fileType_str = "" end
	        fileGrp_str = file[1:3]
	        fileDir_str = string(folder,"/",file)
	        # loops over subfolders, if any exist
	        if (fileType_str == "csv" && fileGrp_str in ("set","par","var","eqn"))
	            files_dic[fileGrp_str] = push!(files_dic[fileGrp_str],fileDir_str)
	        elseif !isfile(fileDir_str)
	            files_dic = readInputFolder([fileDir_str],files_dic)
	        end
	    end
	end

    return files_dic
end

# XXX filters missing and adjusts data according to "all" statements
function convertReadIn(readIn_df::DataFrame,fileName_str::String,set_arr::Array{Symbol},setLngShrt_dic::Dict{Symbol,Symbol},report::Array{Tuple,1},lock_::ReentrantLock,sets::Dict{Symbol,Tree} = Dict{Symbol,Tree}())

	setNames_arr = filterSetColumns(readIn_df,set_arr)
    oprNames_arr = filterSetColumns(readIn_df,[:parameter,:variable,:value, :id])
	readInColAll_tup = tuple(namesSym(readIn_df)...)

	# drop irrelevant column that do not relate to a set or an operator or are completely empty
	select!(readIn_df, Not(setdiff(readInColAll_tup,vcat(setNames_arr[1],setNames_arr[2],oprNames_arr[1]))))
	emptyCol_arr = filter(x -> eltype(readIn_df[!,x]) == Missing,namesSym(readIn_df))
	setNames_arr[1] = setdiff(setNames_arr[1],emptyCol_arr)
	select!(readIn_df, Not(emptyCol_arr))

	# filter value columns
	readInCol_arr = namesSym(readIn_df)
	valCol_arr = filter(x -> occursin("value",string(x)),readInCol_arr)

	# XXX convert missing values and change array container type for editing later
	for j in 1:size(readIn_df,2)
		col = collect(readIn_df[!,j])

		if eltype(col) >: Int
			col = replace(string.(col),"missing" => "")
			if readInCol_arr[j] in valCol_arr
				replace!(col,"" => "NaN")
				col = parse.(Float64,col)
			end
			readIn_df[!,j] = col
		elseif eltype(col) >: Missing
			act_type = eltype(col) >: String ? String : Float64
			# convert remaining columns to strings and replace 'missing' with empty string
			col[findall(ismissing.(col))] .= act_type == String ? "" : NaN
			readIn_df[!,j] = convert(Array{act_type,1},col)
		else
			readIn_df[!,j] = col
		end
	end

	# XXX check types of columns
	if  !isempty(valCol_arr) && any(map(x -> eltype(readIn_df[!,x]),  findall(map(x -> x in valCol_arr, readInCol_arr))) .== String)
		lock(lock_)
		push!(report,(3,"parameter read-in",fileName_str,"detected strings in value column, file was not read-in"))
		unlock(lock_)
		return DataFrame()
	end

	for supCol in findall(eltype.(eachcol(readIn_df)) .!= String)
		if !occursin("value",string(readInCol_arr[supCol]))
			lock(lock_)
			push!(report,(3,"parameter read-in",fileName_str,"entries in $(readInCol_arr[supCol]) could be transferred to string (probably provided as floats), file was not read-in"))
			unlock(lock_)
			return DataFrame()
		end
	end

    # XXX rewrites rows with all commands into full format
    for col in setNames_arr[1]

		colVal_arr = readIn_df[!, col]
        # check column for keywords
        rowsAll_arr = map(x -> length(x) >= 3 && lowercase(x[1:3]) == "all",colVal_arr)

        if all(!,rowsAll_arr) continue end

        # determine relevant reference for "all", if parameter are read in
        if  !isempty(sets) # take reference from readin sets
            specSet_arr = split(String(col),"_")
            relSet_obj = sets[setLngShrt_dic[Symbol(specSet_arr[1])]]
            colValUni_arr = unique(map(x -> x.val,getNodesLvl(relSet_obj, parse(Int,specSet_arr[2]))))
        else # take reference from other column values, relevant when sets are currently read in
            colValUni_arr = sort(unique(filter(x -> !isempty(x),colVal_arr[(!).(rowsAll_arr)])))
        end

        # loop over rows with all
        for row in eachrow(readIn_df[rowsAll_arr,:])
            # append new rows to dataframe
            addRow_df = row
            allInfo_str = reduce(replace, ["all"=>"", "("=>"", ")"=>""], init=addRow_df[col])

            if occursin(":",allInfo_str)
                allVal_arr = split(allInfo_str,":")
                rplVal_arr = colValUni_arr[findall(x->x==allVal_arr[1], colValUni_arr)[1]:findall(x->x==allVal_arr[2], colValUni_arr)[1]]

				# reports if values within all expression could not be matched to sets
				if length(rplVal_arr) != length(allVal_arr)
					lock(lock_)
					push!(report,(2,"parameter read-in",fileName_str,"at least one value within all expression $(allInfo_str) could not be matched to an existing set"))
					unlock(lock_)
				end
            elseif occursin(",",allInfo_str)
                allVal_arr = split(allInfo_str,",")
                rplVal_arr = colValUni_arr[map(x -> in(x,allVal_arr),colValUni_arr)]

				# reports if values within all expression could not be matched to sets
				if length(rplVal_arr) != length(allVal_arr)
					lock(lock_)
					push!(report,(2,"parameter read-in",fileName_str,"at least one value within all expression $(allInfo_str) could not be matched to an existing set"))
					unlock(lock_)
				end
            else
                rplVal_arr = colValUni_arr
            end

            for addVal in rplVal_arr
                addRow_df[col] = addVal
                push!(readIn_df, [addRow_df[col] for col in readInCol_arr])
            end
        end
        #remove inital rows with all#
        delete!(readIn_df,findall(rowsAll_arr))
    end

	# XXX convert column names if sets are defined for multiple insances (e.g. two regions in case of trade related parameters)
	if split(fileName_str,"/")[end][1:3] == "par"
		# creates a dictionary that assigns everything after the set name seperated with a "_" to the respective set
		splitCol_arr = map(x -> split(String(x),"_"),setdiff(namesSym(readIn_df),oprNames_arr[1]))
		setCol_arr = unique(map(x -> Symbol(x[1]),splitCol_arr))
		grpCol_dic = Dict(x => map(z -> z[2:end],filter(y -> String(x) == y[1],splitCol_arr)) for x in setCol_arr)

		# loop over dictionary to check for irregular names and two sets having two columns assigned
		letters_arr = ("b","c","d","e","f","g","h")

		for set in keys(grpCol_dic)
			newCol_dic = Dict{Symbol,Symbol}()

			if any(map(x -> tryparse(Int,x),vcat(grpCol_dic[set]...)) .== nothing)
				lock(lock_)
				push!(report,(3,"parameter read-in",fileName_str,"column for set $(set) does not contain a number after _"))
				unlock(lock_)
				continue
			end

			# get the unique number of "_" that appear in columns assigned to the respective set
			# if the csv file has several columns of the same name JULIA adds a "_1" to the second column when reading in as a dataframe => so this checks if set is defined for multiple instances
			uniLen_arr = unique(map(x -> length(x),grpCol_dic[set]))
			if  length(uniLen_arr) == 2
				setNumbers_arr = map(x -> parse(Int,x[1]),grpCol_dic[set])
				switchBool_arr = fill(false,length(setNumbers_arr))
				for i in 2:length(setNumbers_arr)
					if setNumbers_arr[i-1] >= setNumbers_arr[i]
						switchBool_arr[i] = true
					end
				end
				switchInt_arr = findall(switchBool_arr)
				z = 0
				for (idx,k) in enumerate(grpCol_dic[set])
					if idx in switchInt_arr
						z = z+1
					end
					if z == 0
						newCol_dic[Symbol(join([set,k...],"_"))] = Symbol(set,"_",k[1])
					else
						newCol_dic[Symbol(join([set,k...],"_"))] = Symbol(set,"_",letters_arr[z],"_",k[1])
					end
				end

				DataFrames.rename!(readIn_df,newCol_dic)
			end
		end
	end

	return readIn_df
end

# </editor-fold>

# <editor-fold desc= creation of tree objects for sets"

# XXX creates tree object for set
function createTree(readIn_df::DataFrame, setLoad_sym::Symbol, report::Array{Tuple,1})

	setLoad_str = string(setLoad_sym)
	# create tree object and add the top node
	tree_obj = Tree()
	tree_obj.nodes[0] = Node(0,"none",0,1,Int[])

	# writes values of first column
	firstCol_sym = Symbol(setLoad_str,"_1")
	topNodes_arr =  filter(x -> !isempty(x),convert(Matrix,unique(readIn_df[!,namesSym(readIn_df) .== firstCol_sym])))

	for (idx, node) in enumerate(sort(topNodes_arr))
	    tree_obj.nodes[idx] = Node(idx,node,1,idx,Int[])
		tree_obj.srcTup[(node,)] = [idx]
		tree_obj.up[idx] = 0
	end
	tree_obj.nodes[0].down = collect(keys(tree_obj.up))

	# adds dictionary for occurrence of single strings
	for v in getNodesLvl(tree_obj, 1)
		a = v.val
		if haskey(tree_obj.srcStr,(a,1))
			push!(tree_obj.srcStr[(a,1)],v.idx)
		else
			tree_obj.srcStr[(a,1)] = [v.idx]
		end
	end

	# loop over subsequent columns and add respective tree levels
	height_int = maximum((map(x -> parse(Int,x[end]), filter(x-> (tryparse(Int,string(x[end])) != nothing) && x[1:minimum([length(x),length(setLoad_str)])] .== setLoad_str,[String(namesSym(readIn_df)[i]) for i = 1:size(readIn_df,2)]))))
	for i in 2:height_int
		createTreeLevel!(readIn_df, tree_obj, setLoad_str, i, report)
	end

	# adds max level
	tree_obj.height = height_int

    return tree_obj
end

# XXX adds nodex on level i to tree object
function createTreeLevel!(readIn_df::DataFrame, tree_obj::Tree, setLoad_str::String, i::Int, report::Array{Tuple,1})
	colNames_arr = filter(x -> occursin(setLoad_str,string(x)), namesSym(readIn_df))
	loLvl_Sym = Symbol(setLoad_str,"_",i)

	# removes upper columns with empty values only
	grpCol_arr = filter(x -> x in colNames_arr[1:(i-1)], colNames_arr)
	grpIn_df = unique(readIn_df[readIn_df[!,loLvl_Sym] .!= "",filter(x -> x in colNames_arr[1:i], colNames_arr)])
	grpRel_arr = setdiff(grpCol_arr,filter(x -> [""] == unique(grpIn_df[!,x]), grpCol_arr))
	select!(grpIn_df,Not(setdiff(grpCol_arr,grpRel_arr)))

	# provides the nodes of the lower level grouped by the upper nodes
	lowerNodes_gdf = groupby(grpIn_df,grpRel_arr)

	# checks for nodes wihtout any upper node assigned
	noUp_arr = findall(map(x -> all(x .== ""), [collect(parent(lowerNodes_gdf)[i, groupcols(lowerNodes_gdf)]) for i in lowerNodes_gdf.starts]))
	up_arr = setdiff(1:length(lowerNodes_gdf),noUp_arr)
	if !isempty(noUp_arr)
		noUpVal_arr = setdiff(union(map(x -> collect(x[!,i]),collect(lowerNodes_gdf[noUp_arr]))...),union(map(x -> collect(x[!,i]),collect(lowerNodes_gdf[up_arr]))...))
		for i in noUpVal_arr
			push!(report,(3,"set read-in",setLoad_str,"node named $(i) could not be assigned to an upper node"))
		end
	end

	# XXX assigns the upper nodes by id to strings of corresponding lower nodes
	startLvl_int = parse(Int,string(grpRel_arr[1])[end])
	upToLow_dic = Dict(lookupTupleTree(tuple(collect(lowerNode[1,grpRel_arr])...), tree_obj,startLvl_int)[1] =>  lowerNode[!,loLvl_Sym] for lowerNode in lowerNodes_gdf[up_arr])
	# XXX iterates over dict to write new nodes into tree
	createNodes!(upToLow_dic,tree_obj,i)

	# adds dictionary for occurrence of single strings
	for v in getNodesLvl(tree_obj, i)
		a = v.val
		if haskey(tree_obj.srcStr,(a,i))
			push!(tree_obj.srcStr[(a,i)],v.idx)
		else
			tree_obj.srcStr[(a,i)] = [v.idx]
		end
	end
end

# XXX create specific node on branch
function createNodes!(upToLow_dic::Dict{Int64,SubArray{String,1,Array{String,1},Tuple{Array{Int64,1}},false}},tree_obj::Tree,i::Int)
	upToLowSort_dic = Dict(map(x -> x => upToLow_dic[x] ,sort(collect(keys(upToLow_dic)))))
	for upperNodeId in sortSiblings(collect(keys(upToLowSort_dic)),tree_obj)
		numRow_int = length(tree_obj.nodes) -1
		exUp_int = length(tree_obj.nodes[upperNodeId].down)
		for (idx, lowerNode) in enumerate(sort(upToLowSort_dic[upperNodeId]))
			newIdx_int = numRow_int + idx
			tree_obj.nodes[newIdx_int] = Node(newIdx_int,lowerNode,i,idx+exUp_int,Int[])
			tree_obj.up[newIdx_int] = upperNodeId

			keyStr_tup = fill("",i)
			keyStr_tup[i] = lowerNode
			foreach(x -> keyStr_tup[x[2]] = tree_obj.nodes[x[1]].val, getAncestors(newIdx_int,tree_obj,:tup,1))

			tree_obj.srcTup[tuple(keyStr_tup...)] = [newIdx_int]
		end
		tree_obj.nodes[upperNodeId].down = union(tree_obj.nodes[upperNodeId].down,collect((numRow_int+1):(numRow_int+length(upToLowSort_dic[upperNodeId]))))
	end
end

# </editor-fold>

# <editor-fold desc= read-in of parameter data"

# XXX reads-in parameter data for respective sheet
function writeParameter(parData_df::DataFrame, sets::Dict{Symbol,Tree}, setLngShrt_dic::Dict{Symbol,Symbol}, fileName_str::String, report::Array{Tuple,1},lock_::ReentrantLock)

	setShrtLng_dic = Dict(value => key for (key, value) in setLngShrt_dic)
    set_arr = vcat(collect(setShrtLng_dic[key] for key in keys(sets))...,:id)
    setNames_arr = filterSetColumns(parData_df,set_arr)[1]
	para_dic = Dict{Symbol, DataFrame}()

	# creates array of all levels provided per set grouped by set
	setIndex_arr = map(setdiff(setNames_arr,set_arr)) do x
		splitSet = split(String(x),"_")
		return length(splitSet) == 1 ? (Symbol(splitSet[1]),1) : (Symbol(join(splitSet[1:end-1],"_")),parse(Int,splitSet[end]))
	end

	setIndex_arr = map(y -> (y,map(w -> w[2],filter(z -> z[1] == y,setIndex_arr))),unique(map(x -> x[1],setIndex_arr)))

	# creates array that is later edited to lookup and save set values, entries: set, blank for set values, levels, start level
	setIni_arr = [parEntry(z[1],initializeLookup(maximum(z[2])-minimum(z[2])+1),z[2],minimum(z[2])) for z in setIndex_arr]
    # creates special entry for sets with only one level, because they do not need to have a number at the end
    for i in intersect(set_arr,setNames_arr) push!(setIni_arr,parEntry(i,String[""],Int[],1)) end

	# throws error, if column level exceeds level of the respective set used
	for ele in setIndex_arr
		set = Symbol(split(string(ele[1]),"_")[1]) # extracts just the actual set name, if it has a letter in the end, because set is used multiple times
		if set != :id && sets[setLngShrt_dic[set]].height < maximum(ele[2])
			lock(lock_)
			push!(report,(2,"parameter read-in",fileName_str,"columns provided for $(ele[1]) exceed level of definition, parameter input ignored"))
			unlock(lock_)
			return para_dic
		end
	end

    # assigns relevant columns to sets
    setCol_dic = Dict{Symbol,Array}()
    for i in setIni_arr
		if isempty(i.lvl)
			setCol_dic[i.colSet] = [i.colSet]; i.lvl = [1]
		else
			setCol_dic[i.colSet] = [Symbol(i.colSet,"_",j) for j in i.lvl]
		end
	end

    # determines relevant parameter/value columns
    oprNames_arr = filterSetColumns(parData_df,[:parameter,:variable,:value],true)[1]
    oprLvl_arr = filter(x -> x != nothing,map(x -> tryparse(Int,split(x,"_")[end]),oprNames_arr))
    parVal_arr = isempty(oprLvl_arr) ? [[:parameter,:value]] : [[Symbol("parameter_",j),Symbol("value_",j)] for j in unique(oprLvl_arr)]

    # converts parameter columns to symbols
    for i in parVal_arr parData_df[!,i[1]] = map(x -> Symbol(x),parData_df[!,i[1]]) end

    # loop over rows to read respective parameter values
	convertParameter!(parData_df,sets,setIni_arr,parVal_arr,para_dic,setCol_dic,setLngShrt_dic,fileName_str,report,lock_)

    return para_dic
end

# XXX gets idx from set names and orders all data in dataframe for respective parameter
function convertParameter!(parData_df::DataFrame,sets::Dict{Symbol,Tree},setIni_arr::Array{parEntry,1},parVal_arr::Array{Array{Symbol,1},1},para_dic::Dict{Symbol,DataFrame},setCol_dic::Dict{Symbol,Array},setLngShrt_dic::Dict{Symbol,Symbol},fileName_str::String,report::Array{Tuple,1},lock_::ReentrantLock)
	setShrtLng_dic = Dict(value => key for (key, value) in setLngShrt_dic)
	for row in eachrow(parData_df)
		setId_dic = Dict{Symbol,Union{Int,Array{Int,1}}}()

		# XXX obtains node ids for row
		# overwrites default values for specific row in setIni_arr
		for i in setIni_arr, (index,j) in enumerate(i.lvl)
			i.entry[j+1-i.startLvl] = row[setCol_dic[i.colSet][index]]
		end

		# extracts specific set values and looks them up to obtain the respective node ids, in case of an id column just directly writes value
		relSets_arr = filter(x ->!all(("" .== x.entry) .| (false .== x.entry)),setIni_arr)
		for ele in relSets_arr
			if ele.colSet == :id
				setId_dic[:id] = parse(Int,ele.entry[1])
			else
				split_arr = split(string(ele.colSet),"_")
				setShort_sym = setLngShrt_dic[Symbol(split_arr[1])]
				saveDic_sym = length(split_arr) == 1 ? setShort_sym : Symbol(setShort_sym,"_",split_arr[2])
				setId_dic[saveDic_sym] = lookupTupleTree(tuple(ele.entry...),sets[setShort_sym],ele.startLvl)
			end
		end

		# goes to next iteration and writes report, if any set used is undefined
		if Int[] in values(setId_dic)
			undefinedDim_arr = map(filter(x -> setId_dic[x] == Int[],collect(keys(setId_dic)))) do x
				split_arr = split(String(x),"_")
				setName = setShrtLng_dic[Symbol(split_arr[1])]
				return length(split_arr) == 1 ? setName : Symbol(setName,"_",split_arr[2])
			end
			undefinedSets_arr = map(y -> join(map(z -> string(y.entry[z]," (lvl ",y.lvl[z],")") ,1:length(y.lvl))," > "),filter(x -> x.colSet in undefinedDim_arr,relSets_arr))

			lock(lock_)
			for undef in undefinedSets_arr
				push!(report,(2,"parameter read-in",fileName_str,"values provided for undefined set $(undef...)"))
			end
			unlock(lock_)

			continue
		end

		# creates all possible combinations of found values
		agnNodes_mat = []
		for i in Iterators.product(collect(values(setId_dic))...) push!(agnNodes_mat,collect(i)) end
		ordAgnNodes_ord = hcat(agnNodes_mat...)

		addEntry_df = DataFrame()

		for (index,y) in enumerate(keys(setId_dic)) addEntry_df[!,y] = ordAgnNodes_ord[index,:] end

		# XXX loop over parameter/value columns, prepares and writes
		for i in parVal_arr

			# extract parameter type and value
			par_sym = row[i[1]]
			if par_sym == Symbol()
				continue
			end

			# adds values to dataframe
			if isempty(addEntry_df)
				if isnan(row[i[2]]) continue end
				addEntry_df = DataFrame(val = row[i[2]])
			else
				if isnan(row[i[2]]) continue end
				addEntry_df[!,:val] .= row[i[2]]
			end

			# creates empty dataframe for parameter, if non-existent so far
			if !in(par_sym,keys(para_dic)) para_dic[par_sym] = DataFrame(val = Float64[]) end

			# adds 0 to dictionary for sets the parameter depends on, but that dont appear in the current file/row
			for missKey in setdiff(namesSym(para_dic[par_sym]),namesSym(addEntry_df)) addEntry_df[!,missKey] .= 0 end

			# adds new column to dataframe for respective parameter if required
			rows_int = nrow(para_dic[par_sym])
			for key in setdiff(namesSym(addEntry_df),namesSym(para_dic[par_sym]))
				para_dic[par_sym][!,key] = zeros(Int, rows_int)
			end

			select!(addEntry_df, namesSym(para_dic[par_sym]))
			append!(para_dic[par_sym],addEntry_df)
		end
	end
end

# XXX filters all columns of dataframe that are related to the sets
function filterSetColumns(input_df::DataFrame,input_arr::Array{Symbol},outStr_boo::Bool = false)
    colNames_arr = [String(namesSym(input_df)[i]) for i = 1:size(input_df,2)]

    # filters columns that relate to input array and further splits them based on "_" seperators
    inRelColNames_arr = collect(Iterators.flatten(map(y -> filter(x -> x[1:minimum([length(y),length(x)])] == y, colNames_arr),map(x-> string(x),input_arr))))
    # columns that relate to input because of numbering or identity
    ctrColNames_arr = vcat(filter(x -> isa(tryparse(Int,x[end:end]),Int),inRelColNames_arr),intersect(colNames_arr,map(x -> String(x),input_arr)))
    # remaining columns, solely used to filter mapping set files
    mapColNames_arr = setdiff(inRelColNames_arr,ctrColNames_arr)
    return_arr = [ctrColNames_arr, mapColNames_arr]

    # convert to symbol unless optional input is set to true
    if !outStr_boo return_arr = map(y -> map(x -> Symbol(x),y),return_arr) end
    return return_arr
end

# XXX initializes dictionary that saves lookups in tree
function initializeLookup(size_int::Int)
    Ini_arr = Array{String}(undef,size_int)
    Ini_arr .= ""
    return Ini_arr
end

# </editor-fold>


text:
using Test
using Statistics
using LinearAlgebra
using Random

using EnhancedGJK
using EnhancedGJK: projection_weights, projection_weights_reference
using CoordinateTransformations: IdentityTransformation, Translation
using StaticArrays: SVector
import GeometryTypes
const gt = GeometryTypes
using FileIO

const mesh_dir = joinpath(dirname(@__FILE__), "meshes")

@testset "reference distance" begin
    mesh = load(joinpath(mesh_dir, "base_link.obj"))
    for x in range(0.05, stop=1, length=10)
        for y in range(-1, stop=1, length=10)
            for z in range(-1, stop=1, length=10)
                point = SVector(x, y, z)
                gjk_dist = gjk(mesh, point)
                simple_dist = ReferenceDistance.signed_distance(mesh, point)
                @test isapprox(gjk_dist.signed_distance, simple_dist, atol=2e-4)
            end
        end
    end
end

@testset "reference interior distance" begin
    verts = [gt.Point(1., -1, 1), gt.Point(1., 1, 1), gt.Point(-1., 1, 1), gt.Point(-1., -1, 1),
             gt.Point(1., -1, -1), gt.Point(1., 1, -1), gt.Point(-1., 1, -1), gt.Point(-1., -1, -1)]
    ft = gt.Face{3, Int}
    faces = [ft(1, 2, 5), ft(6, 5, 2),
             ft(2, 3, 6), ft(7, 6, 3),
             ft(3, 4, 7), ft(8, 7, 4),
             ft(4, 1, 8), ft(5, 8, 1),
             ft(2, 1, 3), ft(4, 3, 1),
             ft(5, 6, 7), ft(7, 8, 5)]
    mesh = gt.HomogenousMesh(verts, faces)
    ReferenceDistance.signed_distance(mesh, SVector(0, 0, 0))

    for x in range(-1, stop=1, length=10)
        for y in range(-1, stop=1, length=10)
            for z in range(-1, stop=1, length=10)
                point = SVector(x, y, z)
                expected = -min(1 - abs(x), 1 - abs(y), 1 - abs(z))
                actual = ReferenceDistance.signed_distance(mesh, point)
                @test isapprox(expected, actual, atol=1e-15)
            end
        end
    end
end

@testset "table" begin
    width = 0.5
    thickness = 0.05
    surface_points = Vector{SVector{3, Float64}}()
    for z in [-thickness, thickness]
    for x in [-width, width]
        for y in [-width, width]
            push!(surface_points, SVector(x, y, z))
        end
    end
    end
    geometry = SVector{length(surface_points)}(surface_points)
    point = zeros(SVector{3, Float64})
    for x in range(-width, stop=width, length=21)
        for y in range(-width, stop=width, length=21)
            z = 0.1
            @test isapprox(gjk(geometry,
                point,
                IdentityTransformation(),
                Translation(SVector(x, y, z))).signed_distance, 0.05)
            z = 0.06
            @test isapprox(gjk(geometry,
                    point,
                    IdentityTransformation(),
                    Translation(SVector(x, y, z))).signed_distance,
                0.01,
                atol=1e-12)
            z = 0.05
            @test isapprox(gjk(geometry,
                    point,
                    IdentityTransformation(),
                    Translation(SVector(x, y, z))).signed_distance,
                0.0,
                atol=1e-12)
        end
    end
end

@testset "johnson distance subalgorithm" begin
    include("johnson_distance.jl")
end

@testset "simplex distance" begin
    simplex = SVector{3}(SVector{2, Float64}[[1., 0], [2., 0], [1., 1]])
    pt = SVector(0., 0)
    cache = CollisionCache(simplex, pt);
    result = gjk!(cache, IdentityTransformation(), IdentityTransformation())
    @test isapprox(result.signed_distance, 1.0)
    @test isapprox(result.closest_point_in_body.a, [1.0, 0.0])
    @test isapprox(result.closest_point_in_body.b, [0.0, 0.0])
end

@testset "mesh to mesh" begin
    mesh = load(joinpath(mesh_dir, "r_foot_chull.obj"))
    dx = 1.0
    foot_length = 0.172786 + 0.090933

    cache = CollisionCache(mesh, mesh)
    result = gjk!(cache, IdentityTransformation(), Translation(SVector(dx, 0, 0)))
    @test isapprox(result.signed_distance, dx - foot_length, atol=1e-3)

    cache = CollisionCache(mesh, mesh)
    result = gjk!(cache, Translation(SVector(dx, 0, 0)), IdentityTransformation())
    @test isapprox(result.signed_distance, dx - foot_length, atol=1e-3)

    cache = CollisionCache(mesh, mesh)
    expected_penetration = 0.01
    result = gjk!(cache, IdentityTransformation(), Translation(foot_length - expected_penetration, 0, 0))
    # TODO: penetration distance is inconsistent and inaccurate
    @test result.signed_distance < 0
end

@testset "neighbor mesh to mesh" begin
    mesh = NeighborMesh(load(joinpath(mesh_dir, "r_foot_chull.obj")))
    dx = 1.0
    foot_length = 0.172786 + 0.090933

    cache = CollisionCache(mesh, mesh)
    result = gjk!(cache, IdentityTransformation(), Translation(SVector(dx, 0, 0)))
    @test isapprox(result.signed_distance, dx - foot_length, atol=1e-3)

    cache = CollisionCache(mesh, mesh)
    result = gjk!(cache, Translation(SVector(dx, 0, 0)), IdentityTransformation())
    @test isapprox(result.signed_distance, dx - foot_length, atol=1e-3)

    cache = CollisionCache(mesh, mesh)
    expected_penetration = 0.01
    result = gjk!(cache, IdentityTransformation(), Translation(foot_length - expected_penetration, 0, 0))
    # TODO: penetration distance is inconsistent and inaccurate
    @test result.signed_distance < 0
end

@testset "geometry types" begin
    # Adapted from
    # https://github.com/JuliaGeometry/GeometryTypes.jl/blob/master/test/gjk.jl
    @testset "gjk examples" begin
        c1 = gt.Simplex(gt.Vec(-1.))
        c2 = gt.Simplex(gt.Vec(4.))
        cache = CollisionCache(c1, c2)
        result = gjk!(cache, IdentityTransformation(), IdentityTransformation())
        @test isapprox(result.signed_distance, 5.0)

        c1 = gt.Simplex(gt.Vec(-1.,0,0))
        c2 = gt.Simplex(gt.Vec(4.,0,0))
        cache = CollisionCache(c1, c2)
        result = gjk!(cache, IdentityTransformation(), IdentityTransformation())
        @test isapprox(result.signed_distance, 5.0)

        c1 = gt.FlexibleConvexHull([gt.Vec(0.,0), gt.Vec(0.,1), gt.Vec(1.,0),gt.Vec(1.,1)])
        c2 = gt.Simplex(gt.Vec(4.,0.5))
        cache = CollisionCache(c1, c2)
        result = gjk!(cache, IdentityTransformation(), IdentityTransformation())
        @test isapprox(result.signed_distance, 3.0)

        pt1 = gt.Vec(1,2,3.)
        pt2 = gt.Vec(3,4,5.)
        cache = CollisionCache(pt1, pt2)
        result = gjk!(cache, IdentityTransformation(), IdentityTransformation())
        @test isapprox(result.signed_distance, norm(pt1 - pt2))
    end

    @testset "gjk intersecting lines" begin
        c1 = gt.Simplex(gt.Vec(1,1.), gt.Vec(1, 2.))
        @test gjk(c1, c1).signed_distance == 0.

        c2 = gt.Simplex(gt.Vec(1,1.), gt.Vec(10, 2.))
        @test gjk(c1, c2).signed_distance == 0.
    end
end

@testset "benchmarks" begin
    include("../perf/runbenchmarks.jl")
end


text:
function compute(coins::Vector{Int}, n::Int)::Int
    combinations = vcat([1], repeat([0], n))
    for i ∈ 1:length(coins)
        for j ∈ coins[i] + 1:n + 1
            combinations[j] += combinations[j - coins[i]]
        end
    end
    return combinations[end]
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
import Base.bitstring

function describe(x::Valid{N, ES}, f = println) where {N,ES}
  isempty(x)    && return string("Valid{$N,$ES}(∅)")                               |> f
  isallreals(x) && return string("Valid{$N,$ES}(ℝp)")                              |> f

  #handle special characters.
  if (x.upper == maxpos(Vnum{N,ES}))
    (x.lower == zero(Vnum{N,ES})) && return string("Valid{$N,$ES}(ℝ(+))")          |> f
    (x.lower == minpos(Vnum{N,ES})) && return string("Valid{$N,$ES}(ℝ(+*))") |> f
    (x.lower == -maxpos(Vnum{N,ES})) && return string("Valid{$N,$ES}(ℝ)")         |> f
  end

  if (x.lower == -maxpos(Vnum{N,ES}))
    (x.upper == zero(Vnum{N,ES})) && return string("Valid{$N,$ES}(ℝ(-))")          |> f
    (x.upper == maxneg(Vnum{N,ES})) && return string("Valid{$N,$ES}(ℝ(-*))") |> f
  end

  if N <= 32
    lower_u = isulp(x.lower)
    upper_u = isulp(x.upper)

    lvalue = Float64(glb(x.lower))
    rvalue = Float64(lub(x.upper))

    if (lvalue == rvalue)
      string("Valid{$N,$ES}(", lvalue, " ex)") |> f
    else
      string("Valid{$N,$ES}(", lvalue, lower_u ? " op" : " ex", " → ", rvalue, upper_u ? " op)" : " ex)") |> f
    end
  else
  end
end

function describe(x::Vnum{N,ES}, f = println) where {N,ES}
  lvalue = Float64(glb(x))
  rvalue = Float64(lub(x))

  if (lvalue == rvalue)
    string("Vnum{$N,$ES}(", lvalue, " ex)") |> f
  else
    string("Vnum{$N,$ES}(", lvalue, " op", " → ", rvalue, " op)") |> f
  end
end

function  bitstring(x::Valid{N,ES}) where {N,ES}
  string( bitstring(x.lower), bitstring(x.upper))
end

function  bitstring(x::Valid{N,ES}, separator::AbstractString) where {N,ES}
  string( bitstring(x.lower, separator),"|", bitstring(x.upper, separator))
end

function show(io::IO, x::Valid{N,ES}) where {N,ES}
  print(io, x.lower, " → " , x.upper)
end

export describe

max_issues_repo_name: MohHizzani/SigmoidNumbers

max_forks_repo_name: MohHizzani/SigmoidNumbers


text:
function get_test_trajdata(roadway::Roadway)
    scene1 = Scene([
        Entity(VehicleState(VecSE2(0.0,0.0,0.0), roadway, 10.0), VehicleDef(), 1),
        Entity(VehicleState(VecSE2(3.0,0.0,0.0), roadway, 20.0), VehicleDef(), 2)
    ]
    )
    scene2 = Scene([
        Entity(VehicleState(VecSE2(1.0,0.0,0.0), roadway, 10.0), VehicleDef(), 1),
        Entity(VehicleState(VecSE2(5.0,0.0,0.0), roadway, 20.0), VehicleDef(), 2)
    ]
    )

    trajdata = [scene1, scene2]
    return trajdata
end

@testset "VehicleState" begin
    s = VehicleState(VecSE2(0.0,0.0,0.0), Frenet(NULL_ROADINDEX, 0.0, 0.0, 0.0), 10.0)
    @test isapprox(velf(s).s, 10.0)
    @test isapprox(velf(s).t,  0.0)
    show(IOBuffer(), s.posF)
    show(IOBuffer(), s)

    s = VehicleState(VecSE2(0.0,0.0,0.0), Frenet(NULL_ROADINDEX, 0.0, 0.0, 0.1), 10.0)
    @test isapprox(velf(s).s, 10.0*cos(0.1))
    @test isapprox(velf(s).t, 10.0*sin(0.1))
    @test isapprox(velg(s).x, 10.0)
    @test isapprox(velg(s).y, 0.0)
    @test isapprox(vel(s), 10.0)

    vehdef = VehicleDef(AgentClass.CAR, 5.0, 3.0)
    veh = Entity(s, vehdef, 1)
    @test isapprox(get_footpoint(veh), VecSE2(0.0,0.0,0.0))
    show(IOBuffer(), vehdef)
    show(IOBuffer(), veh)

    ri = RoadIndex(CurveIndex(1,0.1), LaneTag(1,2))
    @test isapprox(Frenet(ri, 0.0, 0.0, 0.1), Frenet(ri, 0.0, 0.0, 0.1))
    @test VehicleState(VecSE2(0.1,0.2,0.3), 1.0) == VehicleState(VecSE2(0.1,0.2,0.3), NULL_FRENET, 1.0)

    @test get_front(veh).x == vehdef.length/2
    @test get_rear(veh).x == -vehdef.length/2

    roadway = gen_straight_roadway(3, 100.0)
    veh = VehicleState(VecSE2(0.0, 0.0, 0.0), roadway, 0.0)
    veh = Entity(veh, VehicleDef(), 1)
    @test get_lane(roadway, veh).tag == LaneTag(1,1)
    @test get_lane(roadway, veh).tag == get_lane(roadway, veh.state).tag
    veh = VehicleState(VecSE2(0.0, 3.0, 0.0), roadway, 0.0)
    @test get_lane(roadway, veh).tag == LaneTag(1,2)
    veh = VehicleState(VecSE2(0.0, 6.0, 0.0), roadway, 0.0)
    @test get_lane(roadway, veh).tag == LaneTag(1,3)
end

@testset "scene" begin
    roadway = get_test_roadway()
    trajdata = get_test_trajdata(roadway)

    s1 = VehicleState(VecSE2(0.0, 0.0, 0.0), roadway, 0.0)
    s2 = VehicleState(VecSE2(5.0, 0.0, 0.0), roadway, 0.0)
    s3 = lerp(s1, s2, 0.5, roadway)
    @test s3.posG.x == 2.5
    @test s3.posF.s == 2.5

    vehstate = VehicleState(VecSE2(0.0, 0.0, 0.0), roadway, 0.0)
    vehstate1 = VehicleState(VecSE2(0.0, 0.0, 0.0), roadway[LaneTag(1,1)], roadway, 0.0)
    @test vehstate1 == vehstate
    veh = Entity(vehstate, VehicleDef(), 1)
    scene1 = Scene([veh])
    scene2 = Scene([veh])
    @test first(scene1.entities) == first(scene2.entities)
    @test scene1.n == scene2.n

    io = IOBuffer()
    show(io, scene1)
    close(io)

    veh2 = Entity(vehstate, BicycleModel(VehicleDef()), 1)
    veh3 = convert(Entity{VehicleState, VehicleDef, Int64}, veh2)
    @test veh3 == veh

    scene = Scene([veh, veh2, veh3])

    io = IOBuffer()
    show(io, scene)
    close(io)
    scene = Scene(typeof(veh))
    copyto!(scene, trajdata[1])
    @test length(scene) == 2
    for (i,veh) in enumerate(scene)
        @test scene[i].state == trajdata[1][i].state
        @test scene[i].def == trajdata[1][i].def
    end

    scene2 = Scene(deepcopy(scene.entities), 2)
    @test length(scene2) == 2
    for (i,veh) in enumerate(scene2)
        @test scene2[i].state == trajdata[1][i].state
        @test scene2[i].def == trajdata[1][i].def
    end

    @test get_by_id(scene, 1) == scene[1]

    empty!(scene2)
    @test length(scene2) == 0

    copyto!(scene2, scene)
    @test length(scene2) == 2
    for (i,veh) in enumerate(scene2)
        @test scene2[i].state == trajdata[1][i].state
        @test scene2[i].def == trajdata[1][i].def
    end

    delete!(scene2, scene2[1])
    @test length(scene2) == 1
    @test scene2[1].state == trajdata[1][2].state
    @test scene2[1].def == trajdata[1][2].def
    scene2[1] = deepcopy(scene[1])
    @test scene2[1].state == trajdata[1][1].state
    @test scene2[1].def == trajdata[1][1].def

    @test findfirst(1, scene) == 1
    @test findfirst(2, scene) == 2
    @test get_first_available_id(scene) == 3

    @test in(1, scene)
    @test in(2, scene)
    @test !in(3, scene)

    veh = scene[2]
    @test veh.state == trajdata[1][2].state
    @test veh.def == trajdata[1][2].def

    push!(scene, trajdata[1][1].state)
end


text: include("timevector.jl")


text: randn(1000) * 0.5 + 1

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Test
using MechGluecode
using MechanicalUnits
using MechanicalUnits: Unitfu.promote_to_derived
@import_expand(m, s, N) # Takes long if defined late)
using DifferentialEquations

@testset "Unitful time with unitless state" begin
    u0 = 30.0s
    tspan = (0.0, 10.0)s
    prob1 = ODEProblem((du,u,t,p) -> (du[1] = -0.2∙s⁻¹ * u[1]), [u0], tspan)
    prob2 = ODEProblem((u,t,p)    -> (-0.2∙s⁻¹ * u[1]), u0, tspan)
    prob3 = ODEProblem((u,t,p)    -> [-0.2∙s⁻¹* u[1]], [u0],tspan)
    for prob in [prob1, prob2, prob3]
        @test solve(prob, Tsit5()).retcode === :Success
    end
end
@testset "ODE quantity" begin
    f = (y, p, t) -> 0.5y / 3.0s  # The derivative, yeah
    u0 = 1.5N
    tspan = (0.0s, 1.0s)
    prob = ODEProblem(f, u0, tspan)
    sol = solve(prob)
    @test sol.t isa Vector{<:Time}
    @test sol.u isa Vector{<:Force}
end