text:
macro check_unique(ex, s = string(ex), msg="")
    quote
        $(esc(__source__))
        vals = $(esc(ex))
        u = unique(vals)
        if length(u) != length(vals)
            rep = Dict{Symbol, Int}()
            for v in vals
                n = get(rep, v, 0) + 1
                rep[v] = n
            end
            rep = filter(((_,n),)-> n > 1, pairs(rep))
            repk = keys(rep)
            throw(ArgumentError($(s)*" names not unique. Repeated names: "*string(repk)*" "*$(esc(msg))))
        end
    end
end

macro check_all_unique(s1, s2)
    quote
        vals = [getproperty($(esc(s1)), :x); getproperty($(esc(s2)), :x)]
        @check_unique vals "x"
        vals = [getproperty($(esc(s1)), :u); getproperty($(esc(s2)), :u)]
        @check_unique vals "u"
        vals = [getproperty($(esc(s1)), :y); getproperty($(esc(s2)), :y)]
        @check_unique vals "y"
    end
end

import ControlSystems as CS
import ControlSystems: nstates, blockdiag

const NameType = Union{Vector{Symbol}, NTuple{N, Symbol} where N}

"""
See `named_ss` for a convenient constructor.
"""
struct NamedStateSpace{T,S} <: AbstractStateSpace{T} where S <: AbstractStateSpace{T}
    sys::S
    x::Vector{Symbol}
    u::Vector{Symbol}
    y::Vector{Symbol}
end

function Base.promote_rule(::Type{U}, ::Type{NamedStateSpace{T, S}}) where
    {T, U<:AbstractStateSpace{T} , S<:AbstractStateSpace{T}}
    inner = promote_type(U,S)
    NamedStateSpace{T, inner}
end

function Base.promote_rule(::Type{NamedStateSpace{T, U}}, ::Type{NamedStateSpace{T, S}}) where
    {T, U<:AbstractStateSpace{T} , S<:AbstractStateSpace{T}}
    inner = promote_type(U,S)
    NamedStateSpace{T, inner}
end

function Base.convert(::Type{NamedStateSpace{T, S}}, s::U) where {T, S <: AbstractStateSpace, U <: AbstractStateSpace}
    s2 = S === U ? s : Base.convert(S, s)
    named_ss(s2, x = gensym("x"), u = gensym("u"), y = gensym("y"))
end

function Base.convert(::Type{NamedStateSpace{T, S}}, s::NamedStateSpace{T, U}) where {T, S <: AbstractStateSpace, U <: AbstractStateSpace}
    sys = Base.convert(S, s.sys)
    NamedStateSpace{T,typeof(sys)}(sys, s.x, s.u, s.y)
end

# Base.convert(::Type{RobustAndOptimalControl.NamedStateSpace{T, S}}, s::S) where {T, S<:RobustAndOptimalControl.NamedStateSpace} = s

const NamedIndex = Union{Symbol, Vector{Symbol}, Colon}

function Base.getproperty(G::NamedStateSpace, s::Symbol)
    s ∈ fieldnames(NamedStateSpace) && (return getfield(G,s))
    return getproperty(G.sys, s)
end

ControlSystems.numeric_type(G::NamedStateSpace) = ControlSystems.numeric_type(G.sys)

"""
    expand_symbol(s::Symbol, n::Int)

Takes a symbol and an integer and returns a vector of symbols with increasing numbers appended to the end. E.g.,
(:x, 3) -> [:x1, :x2, :x3]

The short-hand syntax `s^n` is also available, e.g., `:x^3 == expand_symbol(:x, 3)`.

Useful to create signal names for named systems.
"""
expand_symbol(s::Symbol, n::Int) = n == 1 ? [s] : [Symbol(string(s)*string(i)) for i in 1:n]
expand_symbol(v, n::Int) = v
Base.:^(s::Symbol, n::Int) = expand_symbol(s, n)

"""
    named_ss(sys::AbstractStateSpace{T}; x, u, y)

Create a `NamedStateSpace` system. This kind of system uses names rather than integer indices to refer to states, inputs and outputs

# Arguments:
- `sys`: A system to add names to.
- `x`: A list of symbols with names of the states.
- `u`: A list of symbols with names of the inputs.
- `y`: A list of symbols with names of the outputs.

Default names of signals if none are provided are `x,u,y`.

# Example
```julia
G1 = ss(1,1,1,0)
G2 = ss(1,1,1,0)
s1 = named_ss(G1, x = :x, u = :u1, y=:y1)
s2 = named_ss(G2, x = :z, u = :u2, y=:y2)

s1[:y1, :u1] # Index using symbols

fb = feedback(s1, s2, r = :r) #
```
"""
function named_ss(sys::AbstractStateSpace{T};
    x = [Symbol("x$i") for i in 1:sys.nx],
    u = [Symbol("u$i") for i in 1:size(sys,2)], # sinze is used instead of sys.nu for ExtendedStateSpace
    y = [Symbol("y$i") for i in 1:size(sys,1)],
    ) where T
    x = expand_symbol(x, sys.nx)
    u = expand_symbol(u, size(sys,2))
    y = expand_symbol(y, size(sys,1))
    length(x) == sys.nx ||
        throw(ArgumentError("Length of state names must match sys.nx ($(sys.nx)), got length $(length(x))"))
    length(u) == size(sys,2) ||
        throw(ArgumentError("Length of input names must match size(sys,2) ($(size(sys,2))), got length $(length(u))"))
    length(y) == size(sys,1) ||
        throw(ArgumentError("Length of output names must match size(sys,1) ($(size(sys,1))), got length $(length(y))"))

    @check_unique x "x"
    @check_unique u "u"
    @check_unique y "y"

    NamedStateSpace{T, typeof(sys)}(sys, x, u, y)
end

iterable(s::Symbol) = [s]
iterable(v) = v

function Base.getindex(sys::NamedStateSpace{T,S}, i::NamedIndex, j::NamedIndex) where {T,S}
    i,j = iterable.((i, j))
    ii = i isa Colon ? i : names2indices(i, sys.y)
    jj = j isa Colon ? j : names2indices(j, sys.u)

    return NamedStateSpace{T,S}(
        sys.sys[ii, jj],
        sys.x,
        sys.u[jj],
        sys.y[ii],
    ) |> sminreal
end

function Base.getindex(sys::NamedStateSpace{T,S}, inds...) where {T,S}
    if size(inds, 1) != 2
        error("Must specify 2 indices to index statespace model")
    end
    rows, cols = CS.index2range(inds...)
    return NamedStateSpace{T,S}(
        sys.sys[rows, cols],
        sys.x,
        sys.u[cols],
        sys.y[rows],
    ) |> sminreal
end

function Base.show(io::IO, G::NamedStateSpace)
    print(io, "Named")
    show(io, G.sys)
    print(io, "\nWith state  names: "); println(io, join(G.x, ' '))
    print(io, "     input  names: "); println(io, join(G.u, ' '))
    print(io, "     output names: "); println(io, join(G.y, ' '))
end

function Base.:-(s1::NamedStateSpace{T,S}) where {T <: CS.TimeEvolution, S}
    return NamedStateSpace{T,S}(
        -s1.sys,
        s1.x,
        s1.u,
        s1.y,
    )
end

function Base.:+(s1::NamedStateSpace{T,S}, s2::NamedStateSpace{T,S}) where {T <: CS.TimeEvolution, S}
    return NamedStateSpace{T,S}(
        s1.sys+s2.sys,
        [s1.x; s2.x],
        s1.u,
        s1.y,
    )
end

function Base.:*(s1::NamedStateSpace{T}, s2::NamedStateSpace{T}) where {T <: CS.TimeEvolution}
    @check_all_unique s1 s2
    if s1.u != s2.y
        connection_map = join(["$y -> $u" for (u,y) in zip(s1.u, s2.y) if u != y], '\n')
        @warn "Connected signals have different names\n $connection_map" maxlog=2
    end
    sys = s1.sys*s2.sys
    S = typeof(sys)
    return NamedStateSpace{T,S}(
        sys,
        [s1.x; s2.x],
        s2.u,
        s1.y,
    )
end

function Base.:*(s1::Number, s2::NamedStateSpace{T, S}) where {T <: CS.TimeEvolution, S}
    return NamedStateSpace{T,S}(
        s1*s2.sys,
        s2.x,
        s2.u,
        [Symbol(string(y)*"_scaled") for y in s2.y]
    )
end

function Base.:/(s::NamedStateSpace{T, S}, n::Number) where {T <: CS.TimeEvolution, S}
    (1/n)*s
end
##

function Base.hcat(systems::NamedStateSpace{T,S}...) where {T,S}
    x = reduce(vcat, getproperty.(systems, :x))
    u = reduce(vcat, getproperty.(systems, :u))
    @check_unique x
    @check_unique u
    return NamedStateSpace{T,S}(
        hcat(getproperty.(systems, :sys)...),
        x,
        u,
        systems[1].y,
    )
end

function Base.vcat(systems::NamedStateSpace{T,S}...) where {T,S}
    x = reduce(vcat, getproperty.(systems, :x))
    y = reduce(vcat, getproperty.(systems, :y))
    @check_unique x
    @check_unique y
    return NamedStateSpace{T,S}(
        vcat(getproperty.(systems, :sys)...),
        x,
        systems[1].u,
        y,
    )
end

"""
    measure(s::NamedStateSpace, names)

Return a system with specified states as measurement outputs.
"""
function measure(s::NamedStateSpace, names)
    inds = names2indices(names, s.x)
    A,B = ssdata(s)
    ny = length(inds)
    C = zeros(eltype(A), ny, s.nx)
    for i = 1:ny
        C[i, inds[i]] = 1
    end
    s2 = ss(A,B,C,0, s.timeevol)
    sminreal(named_ss(s2; s.x, s.u, y=names))
end

"""
    feedback(s1::NamedStateSpace, s2::NamedStateSpace;
            w1 = [],  u1 = (:), z1 = (:), y1 = (:),
            w2 = (:), u2 = (:), z2 = (:), y2 = [], kwargs...)

Feedback between two named systems. The lists of signals to connect are expected to be lists of symbols with signal names, `(:)` indicating all signals, or `[]` indicating no signals.

All signal sets `w1,u1,z1,y1,u2,y2,w2,z2` have the same meaning as for the advanced feedback methods for regular systems.
The added signal set `r` is used to optionally provide a new name for the input of the feedback loop.

To simplify creating complicated feedback interconnections, see `connect`.
"""
function ControlSystems.feedback(s1::NamedStateSpace{T}, s2::NamedStateSpace{T};
    u1=:, w1=[],z1=:,y1=:,u2=:,y2=:,w2=:,z2=[], kwargs...) where {T <: CS.TimeEvolution}
    W1 = names2indices(w1, s1.u)
    U1 = names2indices(u1, s1.u)
    Z1 = names2indices(z1, s1.y)
    Y1 = names2indices(y1, s1.y)

    W2 = names2indices(w2, s2.u)
    U2 = names2indices(u2, s2.u)
    Z2 = names2indices(z2, s2.y)
    Y2 = names2indices(y2, s2.y)

    # TODO: add feedthrough if user requests to have some inputs as outputs

    sys = feedback(s1.sys, s2.sys; W1, W2, U1, U2, Z1, Z2, Y1, Y2, kwargs...)
    fbname = gensym(:feedback)
    x1  = [s1.x; s2.x]
    @check_unique x1
    x1 = [Symbol(string(x1)*string(fbname)) for x1 in x1] # add unique name postfix
    @assert sys.nu == length(W1) + length(W2)
    @assert sys.ny == length(Z1) + length(Z2)
    @assert sys.nx == length(x1)
    nsys = NamedStateSpace{T,typeof(sys)}(sys, x1, s1.u[[W1; W2]], s1.y[[Z1; Z2]])
    sminreal(nsys)
end

"""
    connect(systems, connections; w1, z1 = (:), verbose = true, kwargs...)

Create complicated feedback interconnection.

Addition and subtraction nodes are achieved by creating a linear combination node, i.e., a system with a `D` matrix only.

# Arguments:
- `systems`: A vector of named systems to be connected
- `connections`: a vector of pairs indicating output => input mappings.
    - `u1`: input mappings  (alternative input argument)
    - `y1`: output mappings (alternative input argument)
- `w1`: external signals
- `z1`: outputs (can overlap with `y1`)
- `verbose`: Issue warnings for signals that have no connection

Example:
The following complicated feedback interconnection

```
                 yF
              ┌────────────────────────────────┐
              │                                │
    ┌───────┐ │  ┌───────┐ yR   ┌─────────┐    │    ┌───────┐
uF  │       │ │  │       ├──────►         │ yC │  uP│       │    yP
────►   F   ├─┴──►   R   │      │    C    ├────+────►   P   ├────┬────►
    │       │    │       │   ┌──►         │         │       │    │
    └───────┘    └───────┘   │  └─────────┘         └───────┘    │
                             │                                   │
                             └───────────────────────────────────┘
```
can be created by
```
F = named_ss(ssrand(1, 1, 2, proper=true), x=:xF, u=:uF, y=:yF)
R = named_ss(ssrand(1, 1, 2, proper=true), x=:xR, u=:uR, y=:yR)
C = named_ss(ssrand(1, 1, 2, proper=true), x=:xC, u=:uC, y=:yC)
P = named_ss(ssrand(1, 1, 3, proper=true), x=:xP, u=:uP, y=:yP)

addP = sumblock("uP = yF + yC") # Sum node before P
addC = sumblock("uC = yR - yP") # Sum node before C

connections = [
    :yP => :yP # Output to input
    :uP => :uP
    :yC => :yC
    :yF => :yF
    :yF => :uR
    :uC => :uC
    :yR => :yR
]
w1 = [:uF] # External inputs

G = connect([F, R, C, P, addP, addC], connections; w1)
```

If an external input is to be connected to multiple points, use a `splitter` to split up the signal into a set of unique names which are then used in the connections.
"""
function connect(systems; u1::Vector{Symbol}, y1::Vector{Symbol}, w1::Vector{Symbol}, z1 = (:), verbose = true, kwargs...)
    full = append(systems...)
    @assert length(y1) == length(u1)
    @check_unique u1 u1 "Connected inputs not unique. If you want to connect several signals to the same input, use a summation node, e.g., named_ss(ss([1  1]), u=[:u1, :u2], y=:usum)"
    @check_unique full.u "system inputs"
    @check_unique full.y "system outputs"

    if verbose
        leftover_inputs = setdiff(full.u, [u1; w1])
        isempty(leftover_inputs) || @warn("The following inputs were unconnected $leftover_inputs")
        leftover_outputs = setdiff(full.y, z1 == (:) ? y1 : [y1; z1])
        isempty(leftover_outputs) || @warn("The following outputs were unconnected $leftover_outputs")
    end

    z2 = []
    w2 = []

    # Connections
    y2 = (:)
    u2 = (:)

    fb = named_ss(ss(I(length(y1)), full.timeevol))
    G = feedback(full, fb; z1, z2, w1, w2, u1, u2, y1, y2, pos_feedback=true, kwargs...)
end

function connect(systems, pairs::AbstractVector{<:Pair}; kwargs...)
    connect(systems; u1 = last.(pairs), y1 = first.(pairs), kwargs...)
end

function splitter(u::Symbol, n::Int, timeevol = Continuous())
    named_ss(ss(ones(n), timeevol), u = [u], y = u^n)
end

# function sumblock(ex::Expr; Ts=0, n=1)
#     @assert ex.head == :(=)
#     sumname = ex.args[1]::Symbol
#     rhs = ex.args[2]::Expr
#     op, s1, s2 = rhs.args
#     timeevol = Ts <= 0 ? ControlSystems.Continuous() : ControlSystems.Discrete(Ts)
#     s = op == :(-) ? -1 : 1
#     named_ss(ss([I(n) s*I(n)], timeevol), u=[s1^n; s2^n], y=sumname)
# end

"""
    sumblock(ex::String; Ts = 0, n = 1)

Create a summation node that sums (or subtracts) vectors of length `n`.

# Arguments:
- `Ts`: Sample time
- `n`: The length of the input and output vectors. Set `n=1` for scalars.

# Examples:
```
julia> sumblock("uP = vf + yL")
NamedStateSpace{Continuous, Int64}
D =
 1  1

With state  names:
     input  names: vf yL
     output names: uP

julia> sumblock("x_diff = xr - xh"; n=3)
NamedStateSpace{Continuous, Int64}
D =
 1  0  0  -1   0   0
 0  1  0   0  -1   0
 0  0  1   0   0  -1

With state  names:
     input  names: xr1 xr2 xr3 xh1 xh2 xh3
     output names: x_diff1 x_diff2 x_diff3

julia> sumblock("a = b + c - d")
NamedStateSpace{Continuous, Int64}
D =
 1  1  -1

With state  names:
     input  names: b c d
     output names: a
```
"""
function sumblock(ex::String; Ts=0, n=1)
    timeevol = Ts <= 0 ? ControlSystems.Continuous() : ControlSystems.Discrete(Ts)
    sumname = Symbol(strip(only(match(r"(.+?)\s?=", ex).captures)))^n
    rhs = split(ex, '=', keepempty=false)[2] |> strip
    rhs = replace(rhs, r"([\+\-])" => s" \1 ") # insert whitespace to make sure split works below.
    s = 1
    mats = []
    names = []
    for sym in split(rhs, ' ', keepempty=false)
        if sym == "+"
            s = 1
        elseif sym == "-"
            s = -1
        else
            push!(mats, s*I(n))
            push!(names, Symbol(sym)^n)
        end
    end
    names = reduce(vcat, names)
    D = reduce(hcat, mats)
    named_ss(ss(D, timeevol), u=names, y=sumname)
end

function ControlSystems.sminreal(s::NamedStateSpace)
    local sys
    try
        sys = sminreal(s.sys)
    catch
        return s
    end
    _, _, _, inds = CS.struct_ctrb_obsv(s.sys) # we do this one more time to get the inds. This implies repeated calculations, but will allow inner systems of exotic types that have a special method for sminreal to keep their type.
    named_ss(sys; x=s.x[inds], s.u, s.y)
end

names2indices(::Colon, allnames) = 1:length(allnames)

function names2indices(names, allnames)
    inds = [findfirst(==(n), allnames) for n in names]
    for i in eachindex(inds)
        inds[i] === nothing && error("The indexed NamedSystem has no signal named $(names[i]), available names are $(allnames)")
    end
    inds
end
function names2indices(name::Symbol, allnames)
    i = findfirst(==(name), allnames)
    i === nothing && error("The indexed NamedSystem has no signal named $name, available names are $(allnames)")
    i:i # return a vector rather than scalar for slices of matrices to not drop dim
end

function ExtendedStateSpace(P::NamedStateSpace; z=[], y=[], w=[], u=[])
    zi = names2indices(z, P.y)
    yi = names2indices(y, P.y)
    wi = names2indices(w, P.u)
    ui = names2indices(u, P.u)
    ss(P.A, P.B[:, wi], P.B[:, ui], P.C[zi, :], P.C[yi, :],
        P.D[zi, wi], P.D[zi, ui], P.D[yi, wi], P.D[yi, ui], P.timeevol)
end

function named_ss(sys::ExtendedStateSpace{T};
    x = [Symbol("x$i") for i in 1:sys.nx],
    u = [Symbol("u$i") for i in 1:sys.nu],
    y = [Symbol("y$i") for i in 1:sys.ny],
    w = [Symbol("w$i") for i in 1:sys.nw],
    z = [Symbol("z$i") for i in 1:sys.nz],
    ) where T
    x = expand_symbol(x, sys.nx)
    u = expand_symbol(u, sys.nu)
    y = expand_symbol(y, sys.ny)
    w = expand_symbol(w, sys.nw)
    z = expand_symbol(z, sys.nz)
    length(x) == sys.nx ||
        throw(ArgumentError("Length of state names must match sys.nx ($(sys.nx))"))
    length(u) == sys.nu ||
        throw(ArgumentError("Length of input names must match sys.nu ($(sys.nu))"))
    length(y) == sys.ny ||
        throw(ArgumentError("Length of output names must match sys.ny ($(sys.ny))"))
    length(w) == sys.nw ||
        throw(ArgumentError("Length of disturbance names must match sys.nw ($(sys.nw))"))
    length(z) == sys.nz ||
        throw(ArgumentError("Length of performance names must match sys.nz ($(sys.nz))"))

    sys2 = ss(sys)
    NamedStateSpace{T, typeof(sys2)}(sys2, x, [w; u], [z; y])
end

CS.@recipe function plot(res::CS.SimResult{<:Any, <:Any, <:Any, <:Any, <:NamedStateSpace})
    s = res.sys
    label --> permutedims(["From $n" for n in s.u])
    yguide --> permutedims(["$n" for n in s.y])
    res2 = CS.SimResult(ntuple(i->getfield(res, i), 4)..., s.sys) # new SimResult without NamedSS to avoid infinite recursion
    CS.@series begin
        res2
    end
end

function CS.bodeplot(s::NamedStateSpace, args...;
    title  = permutedims(["From $n" for n in s.u]),
    kwargs...)
    bodeplot(s.sys, args...; kwargs...)
    CS.Plots.plot!(;
        title,
        ylabel = permutedims(["$n" for n in s.y]),
    )
end

function CS.c2d(s::NamedStateSpace, Ts::Real, args...;
    kwargs...)
    named_ss(c2d(s.sys, Ts, args...; kwargs...); s.x, s.u, s.y)
end

function CS.append(systems::NamedStateSpace...)
    systype = promote_type([typeof(s.sys) for s in systems]...)
    timeevol = common_timeevol(systems...)
    A = blockdiag([s.A for s in systems]...)
    B = blockdiag([_remove_empty_cols(s.B) for s in systems]...)
    C = blockdiag([s.C for s in systems]...)
    D = blockdiag([_remove_empty_cols(s.D) for s in systems]...)

    # TODO: If some name below is repeated, perhaps the B matrix can be reduced in width?
    # This could be optional, i.e., unique_names=true
    x = reduce(vcat, getproperty.(systems, :x))
    y = reduce(vcat, getproperty.(systems, :y))
    u = reduce(vcat, getproperty.(systems, :u))

    return named_ss(systype(A, B, C, D, timeevol); x, y, u)
end


text:
export Tree
export data, parent, children, isroot, is_parents_last_child, leaves

mutable struct Tree{T}
    data::T
    parent::Union{Nothing, Tree{T}}
    children::Union{Nothing, Vector{Tree{T}}}
end

data(node::Tree) = node.data
parent(node::Tree) = node.parent
children(node::Tree) = node.children

Tree(data::T) where {T} = Tree{T}(data, nothing, nothing)
function Tree(data::T, parent::Tree{T}) where{T}
    this = Tree{T}(data, parent, nothing)
    if isnothing(children(parent))
        parent.children = [this]
    else
        push!(children(parent), this)
    end
    return this
end

isroot(node::Tree) = parent(node) === nothing
is_parents_last_child(node::Tree) = children(parent(node))[end] === node
function leaves(node::Tree{T}) where {T}
    leaf_nodes = Tree{T}[]
    node_children = children(node)
    if !isnothing(node_children)
        for child ∈ node_children
            get_leaves!(child, leaf_nodes)
        end
    else
        push!(leaf_nodes, node)
    end
    return leaf_nodes
end
function get_leaves!(node::Tree{T}, leaf_nodes::Vector{Tree{T}}) where {T}
    node_children = children(node)
    if !isnothing(node_children)
        for child ∈ node_children
            get_leaves!(child, leaf_nodes)
        end
    else
        push!(leaf_nodes, node)
    end
    return nothing
end

function Base.show(io::IO, node::Tree)
    println(io, data(node))
    node_children = children(node)
    if !isnothing(node_children)
        nchildren = length(node_children)
        if 7 < nchildren
            show(io, node_children[1], "")
            show(io, node_children[2], "")
            println(io, "│  ⋮")
            println(io, "│  ⋮ (", nchildren-4, " additional children)")
            println(io, "│  ⋮")
            show(io, node_children[end-1], "")
            show(io, node_children[end],   "")
        else
            for child ∈ node_children
                show(io, child, "")
            end
        end
    end
end
function Base.show(io::IO, node::Tree, predecessor_string::String)
    next_predecessor_string = ""
    if !isroot(node)
        if is_parents_last_child(node)
            print(io, predecessor_string * "└─ ")
            next_predecessor_string = predecessor_string * "   "
        else
            print(io, predecessor_string * "├─ ")
            next_predecessor_string = predecessor_string * "│  "
        end
    end
    println(io, data(node))
    node_children = children(node)
    if !isnothing(node_children)
        nchildren = length(node_children)
        if 7 < nchildren
            show(io, node_children[1], next_predecessor_string)
            show(io, node_children[2], next_predecessor_string)
            println(io, next_predecessor_string * "│  ⋮")
            println(io, next_predecessor_string * "│  ⋮ (", nchildren-4,
                    " additional children)")
            println(io, next_predecessor_string * "│  ⋮")
            show(io, node_children[end-1], next_predecessor_string)
            show(io, node_children[end], next_predecessor_string)
        else
            for child ∈ node_children
                show(io, child, next_predecessor_string)
            end
        end
    end
end

ext: jl

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
enum Color {
  RED,
  GREEN,
  BLUE
};

int main(){
  Color color = Color.RED;
  if(color == Color.RED){
    printString("red");
  } else{
    printString("not red");
  }
  return 0;
}

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCellDG
include("../useful_routines.jl")

function test_transformation_strain_simple_tension()
    L = 1.0
    W = 1.0
    lambda1, mu1 = 1.0, 2.0
    K1 = (lambda1 + 2mu1 / 3)
    lambda2, mu2 = 2.0, 4.0
    theta0 = -0.067
    e22 = K1 * theta0 / (lambda1 + 2mu1)
    penaltyfactor = 1e2
    nelmts = 2
    dx = 1.0 / nelmts
    penalty = penaltyfactor / dx * (lambda1 + mu1)
    eta = 1
    stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)
    transfstress =
        CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)

    polyorder = 1
    numqp = required_quadrature_order(polyorder) + 2

    elasticitybasis = LagrangeTensorProductBasis(2, polyorder)
    levelsetbasis = HermiteTensorProductBasis(2)

    dim, nf = size(interpolation_points(levelsetbasis))
    refpoints = interpolation_points(elasticitybasis)

    cgmesh = CutCellDG.CGMesh([0.0, 0.0], [L, W], [nelmts, nelmts], nf)
    mesh = CutCellDG.DGMesh([0.0, 0.0], [L, W], [nelmts, nelmts], refpoints)

    levelset =
        CutCellDG.LevelSet(x -> 1.0, cgmesh, levelsetbasis)

    cutmesh = CutCellDG.CutMesh(mesh, levelset)
    cellquads = CutCellDG.CellQuadratures(cutmesh, levelset, numqp)
    interfacequads = CutCellDG.InterfaceQuadratures(cutmesh, levelset, numqp)
    facequads = CutCellDG.FaceQuadratures(cutmesh, levelset, numqp)

    sysmatrix = CutCellDG.SystemMatrix()
    sysrhs = CutCellDG.SystemRHS()

    CutCellDG.assemble_displacement_bilinear_forms!(
        sysmatrix,
        elasticitybasis,
        cellquads,
        stiffness,
        cutmesh,
    )
    CutCellDG.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        elasticitybasis,
        cellquads,
        cutmesh,
    )
    CutCellDG.assemble_interelement_condition!(
        sysmatrix,
        elasticitybasis,
        facequads,
        stiffness,
        cutmesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_interelement_transformation_linear_form!(
        sysrhs,
        transfstress,
        elasticitybasis,
        facequads,
        cutmesh,
    )

    CutCellDG.assemble_penalty_displacement_component_bc!(
        sysmatrix,
        sysrhs,
        x -> 0.0,
        elasticitybasis,
        facequads,
        stiffness,
        cutmesh,
        x -> x[1] ≈ 0.0,
        [1.0, 0.0],
        penalty,
    )
    CutCellDG.assemble_penalty_displacement_component_transformation_linear_form!(
        sysrhs,
        transfstress,
        elasticitybasis,
        facequads,
        cutmesh,
        x -> x[1] ≈ 0.0,
        [1.0, 0.0],
    )

    CutCellDG.assemble_penalty_displacement_component_bc!(
        sysmatrix,
        sysrhs,
        x -> 0.0,
        elasticitybasis,
        facequads,
        stiffness,
        cutmesh,
        x -> x[2] ≈ 0.0,
        [0.0, 1.0],
        penalty,
    )
    CutCellDG.assemble_penalty_displacement_component_transformation_linear_form!(
        sysrhs,
        transfstress,
        elasticitybasis,
        facequads,
        cutmesh,
        x -> x[2] ≈ 0.0,
        [0.0, 1.0],
    )

    CutCellDG.assemble_penalty_displacement_component_bc!(
        sysmatrix,
        sysrhs,
        x -> 0.0,
        elasticitybasis,
        facequads,
        stiffness,
        cutmesh,
        x -> x[1] ≈ L,
        [1.0, 0.0],
        penalty,
    )
    CutCellDG.assemble_penalty_displacement_component_transformation_linear_form!(
        sysrhs,
        transfstress,
        elasticitybasis,
        facequads,
        cutmesh,
        x -> x[1] ≈ L,
        [1.0, 0.0],
    )

    matrix = CutCellDG.sparse_displacement_operator(sysmatrix, cutmesh)
    rhs = CutCellDG.displacement_rhs_vector(sysrhs, cutmesh)

    sol = matrix \ rhs
    disp = reshape(sol, 2, :)

    nodalcoordinates = CutCellDG.nodal_coordinates(mesh)

    testdisp = copy(nodalcoordinates)
    testdisp[1, :] .= 0.0
    testdisp[2, :] .*= e22

    @test allapprox(disp, testdisp, 1e2eps())
end

test_transformation_strain_simple_tension()


text:
"""
    NonTraditionalBetaPlane{FT} <: AbstractRotation

A Coriolis implementation that accounts for the latitudinal variation of both
the locally vertical and the locally horizontal components of the rotation vector.
The "traditional" approximation in ocean models accounts for only the locally
vertical component of the rotation vector (see [`BetaPlane`](@ref)).

This implementation is based off of section 5 of Dellar (2011). It conserve energy,
angular momentum, and potential vorticity.

References
==========
Dellar, P. (2011). Variations on a beta-plane: Derivation of non-traditional
    beta-plane equations from Hamilton's principle on a sphere. Journal of
    Fluid Mechanics, 674, 174-195. doi:10.1017/S0022112010006464
"""
struct NonTraditionalBetaPlane{FT} <: AbstractRotation
    fz :: FT
    fy :: FT
    β  :: FT
    γ  :: FT
    R  :: FT
end

"""
    NonTraditionalBetaPlane(FT=Float64;
        fz=nothing, fy=nothing, β=nothing, γ=nothing,
        rotation_rate=Ω_Earth, latitude=nothing, radius=R_Earth)

The user may directly specify `fz`, `fy`, `β`, `γ`, and `radius` or the three
parameters `rotation_rate`, `latitude`, and `radius` that specify the rotation rate
and radius of a planet, and the central latitude (where y = 0) at which the
non-traditional `β`-plane approximation is to be made.

By default, the `rotation_rate` and planet `radius` is assumed to be Earth's.
"""
function NonTraditionalBetaPlane(FT=Float64;
        fz=nothing, fy=nothing, β=nothing, γ=nothing,
        rotation_rate=Ω_Earth, latitude=nothing, radius=R_Earth)

    Ω, φ, R =  rotation_rate, latitude, radius

    use_f = !all(isnothing.((fz, fy, β, γ))) && isnothing(latitude)
    use_planet_parameters = !isnothing(latitude) && all(isnothing.((fz, fy, β, γ)))

    if !xor(use_f, use_planet_parameters)
        throw(ArgumentError("Either the keywords fz, fy, β, γ, and radius must be specified, " *
                            "*or* all of rotation_rate, latitude, and radius."))
    end

    if use_planet_parameters
        fz =  2Ω*sind(φ)
        fy =  2Ω*cosd(φ)
        β  =  2Ω*cosd(φ)/R
        γ  = -4Ω*sind(φ)/R
    end

    return NonTraditionalBetaPlane{FT}(fz, fy, β, γ, R)
end

@inline two_Ωʸ(P, y, z) = P.fy * (1 -  z/P.R) + P.γ * y
@inline two_Ωᶻ(P, y, z) = P.fz * (1 + 2z/P.R) + P.β * y

# This function is eventually interpolated to fcc to contribute to x_f_cross_U.
@inline two_Ωʸw_minus_two_Ωᶻv(i, j, k, grid, coriolis, U) =
    @inbounds (  two_Ωʸ(coriolis, grid.yC[j], grid.zC[k]) * ℑzᵃᵃᶜ(i, j, k, grid, U.w)
               - two_Ωᶻ(coriolis, grid.yC[j], grid.zC[k]) * ℑyᵃᶜᵃ(i, j, k, grid, U.v))

@inline x_f_cross_U(i, j, k, grid, coriolis::NonTraditionalBetaPlane, U) =
    ℑxᶠᵃᵃ(i, j, k, grid, two_Ωʸw_minus_two_Ωᶻv, coriolis, U)

@inline y_f_cross_U(i, j, k, grid, coriolis::NonTraditionalBetaPlane, U) =
    @inbounds  two_Ωᶻ(coriolis, grid.yF[k], grid.zC[k]) * ℑxyᶜᶠᵃ(i, j, k, grid, U.u)

@inline z_f_cross_U(i, j, k, grid, coriolis::NonTraditionalBetaPlane, U) =
    @inbounds -two_Ωʸ(coriolis, grid.yC[j], grid.zF[k]) * ℑxzᶜᵃᶠ(i, j, k, grid, U.u)

Base.show(io::IO, β_plane::NonTraditionalBetaPlane{FT}) where FT =
    print(io, "NonTraditionalBetaPlane{$FT}: ",
          @sprintf("fz = %.2e, fy = %.2e, β = %.2e, γ = %.2e, R = %.2e",
                   β_plane.fz, β_plane.fy, β_plane.β, β_plane.γ, β_plane.R))


text:
# distritbutions
@rule typeof(-)(:out, Marginalisation) (m_in1::UnivariateNormalDistributionsFamily, m_in2::UnivariateNormalDistributionsFamily) = begin
    min1, vin1 = mean_var(m_in1)
    min2, vin2 = mean_var(m_in2)
    return NormalMeanVariance(min1 - min2, vin1 + vin2)
end

@rule typeof(-)(:out, Marginalisation) (m_in1::MultivariateNormalDistributionsFamily, m_in2::MultivariateNormalDistributionsFamily) = begin
    min1, vin1 = mean_cov(m_in1)
    min2, vin2 = mean_cov(m_in2)
    return MvNormalMeanCovariance(min1 - min2, vin1 + vin2)
end

# PointMass
@rule typeof(-)(:out, Marginalisation) (m_in1::PointMass, m_in2::PointMass) = PointMass(mean(m_in1) - mean(m_in2))

@rule typeof(-)(:out, Marginalisation) (m_in1::NormalMeanPrecision, m_in2::PointMass) = begin
    min1, win1 = mean_precision(m_in1)
    return NormalMeanPrecision(min1 - mean(m_in2), win1)
end

@rule typeof(-)(:out, Marginalisation) (m_in1::PointMass, m_in2::NormalMeanPrecision) = begin
    min2, win2 = mean_precision(m_in2)
    return NormalMeanPrecision(mean(m_in1) - min2, win2)
end

@rule typeof(-)(:out, Marginalisation) (m_in1::MvNormalMeanPrecision, m_in2::PointMass) = begin
    min1, win1 = mean_precision(m_in1)
    return MvNormalMeanPrecision(min1 - mean(m_in2), win1)
end

@rule typeof(-)(:out, Marginalisation) (m_in1::PointMass, m_in2::MvNormalMeanPrecision) = begin
    min2, win2 = mean_precision(m_in2)
    return MvNormalMeanPrecision(mean(m_in1) - min2, win2)
end

@rule typeof(-)(:out, Marginalisation) (m_in1::UnivariateNormalDistributionsFamily, m_in2::PointMass) = begin
    min1, vin1 = mean_var(m_in1)
    return NormalMeanVariance(min1 - mean(m_in2), vin1)
end

@rule typeof(-)(:out, Marginalisation) (m_in1::PointMass, m_in2::UnivariateNormalDistributionsFamily) = begin
    min2, vin2 = mean_var(m_in2)
    return NormalMeanVariance(mean(m_in1) - min2, vin2)
end

@rule typeof(-)(:out, Marginalisation) (m_in1::MultivariateNormalDistributionsFamily, m_in2::PointMass) = begin
    min1, vin1 = mean_cov(m_in1)
    return MvNormalMeanCovariance(min1 - mean(m_in2), vin1)
end

@rule typeof(-)(:out, Marginalisation) (m_in1::PointMass, m_in2::MultivariateNormalDistributionsFamily) = begin
    min2, vin2 = mean_cov(m_in2)
    return MvNormalMeanCovariance(mean(m_in1) - min2, vin2)
end


text:
using ModelingToolkitStandardLibrary.Thermal, ModelingToolkit, OrdinaryDiffEq, Test
@parameters t

# Test HeatCapacitor, TemperatureSensor, RelativeTemperatureSensor, FixedTemperature
@testset "Heat systems" begin
    T, C, G = 10.0, 10.0, 10.0
    @variables final_T(t)
    @named mass1         = HeatCapacitor(C=C)
    @named mass2         = HeatCapacitor(C=C)
    @named th_conductor  = ThermalConductor(G=G)
    @named reltem_sensor = RelativeTemperatureSensor()
    @named T_sensor1     = TemperatureSensor()
    @named T_sensor2     = TemperatureSensor()
    @named tem_src       = FixedTemperature(T=T)

    @info "Building a single-body system..."
    eqs = [
        connect(mass1.port, th_conductor.port_a)
        connect(th_conductor.port_b, reltem_sensor.port_a)
        connect(reltem_sensor.port_b, tem_src.port)
    ]
    @named h1 = ODESystem(eqs, t, systems=[mass1, reltem_sensor, tem_src, th_conductor])
    sys = structural_simplify(h1)

    u0 = [
        mass1.T => 2.0
    ]
    prob = ODEProblem(sys, u0, (0, 2.0))
    sol = solve(prob, Rodas4())

    # Check if Relative temperature sensor reads the temperature of heat capacitor
    # when connected to a thermal conductor and a fixed temperature source
    @test sol[reltem_sensor.T] + sol[tem_src.port.T] == sol[mass1.T] + sol[th_conductor.dT]

    @info "Building a two-body system..."
    eqs = [
        connect(T_sensor1.port, mass1.port, th_conductor.port_a)
        connect(th_conductor.port_b, mass2.port, T_sensor2.port)
        final_T ~ (mass1.C * mass1.T + mass2.C * mass2.T) /
        (mass1.C + mass2.C)
    ]
    @named h2 = ODESystem(eqs, t, [final_T], [],
                        systems=[mass1, mass2, T_sensor1, T_sensor2, th_conductor])
    sys = structural_simplify(h2)

    u0 = [
        mass1.T => 1.0
        mass2.T => 10.0
        final_T => 12
    ]
    prob = ODEProblem(sys, u0, (0, 3.0))
    sol = solve(prob, Rodas4())

    m1, m2 = sol.u[end]
    @test m1 ≈ m2 atol=1e-1
    mass_T = reduce(hcat, sol.u)
    @test sol[T_sensor1.T] == mass_T[1, :]
    @test sol[T_sensor2.T] == mass_T[2, :]
end

# Test HeatFlowSensor, FixedHeatFlow, ThermalResistor, ThermalConductor
@testset "Heat flow system" begin
    C, G, R = 10, 10, 10
    @named flow_src     = FixedHeatFlow(Q_flow=50, alpha=100)
    @named mass1        = HeatCapacitor(C=C)
    @named hf_sensor1   = HeatFlowSensor()
    @named hf_sensor2   = HeatFlowSensor()
    @named th_conductor = ThermalConductor(G=G)
    @named th_resistor  = ThermalResistor(R=R)
    @named th_ground    = FixedTemperature(T=0)

    @info "Building a heat-flow system..."
    eqs = [
        connect(mass1.port, th_resistor.port_a, th_conductor.port_a)
        connect(th_conductor.port_b, flow_src.port, hf_sensor1.port_a, hf_sensor2.port_a)
        connect(th_resistor.port_b, hf_sensor1.port_b, hf_sensor2.port_b, th_ground.port)
    ]
    @named h2 = ODESystem(eqs, t,
                          systems=[mass1, hf_sensor1, hf_sensor2,
                                   th_resistor, flow_src, th_ground, th_conductor])
    sys = structural_simplify(h2)

    u0 = [
        mass1.T            => 10.0
        th_resistor.Q_flow => 1.0
    ]
    prob = ODEProblem(sys, u0, (0, 3.0))
    sol  = solve(prob, Rodas4())

    @test sol[th_conductor.dT] .* G == sol[th_conductor.Q_flow]
    @test sol[th_conductor.Q_flow] ≈ sol[hf_sensor1.Q_flow] + sol[flow_src.port.Q_flow]

    @test sol[mass1.T] == sol[th_resistor.port_a.T]
    @test sol[th_resistor.dT] ./ R ≈ sol[th_resistor.Q_flow]

end

# Tests ConvectiveResistor and includes FixedTemperature and ThermalResistor
@testset "Piston cylinder wall" begin
    Tᵧ, Tᵪ = 1000, 10 # ᵧ -> gas and ᵪ -> coolant
    Rᵧ, Rᵪ = 50e-4, 10e-4 # R = 1/h; h is convection co-efficient
    R_wall = 1.5e-4
    @named coolant     = ConvectiveResistor(R=Rᵪ)
    @named gas         = ConvectiveResistor(R=Rᵧ)
    @named wall        = ThermalResistor(R=R_wall)
    @named gas_tem     = FixedTemperature(T=Tᵧ)
    @named coolant_tem = FixedTemperature(T=Tᵪ)

    @info "Building a piston-cylinder..."
    eqs = [
        connect(gas_tem.port, gas.solidport)
        connect(gas.fluidport, wall.port_a)
        connect(wall.port_b, coolant.fluidport)
        connect(coolant.solidport, coolant_tem.port)
    ]
    @named piston = ODESystem(eqs, t, systems=[gas_tem, wall, gas, coolant, coolant_tem])
    sys = structural_simplify(piston)

    u0 = [
        coolant.dT  => 5.0
        wall.Q_flow => 10.0
    ]
    prob = ODEProblem(sys, u0, (0, 3.0))
    sol = solve(prob, Rodas4())

    # Heat-flow-rate is equal in magnitude
    # and opposite in direction
    @test sol[gas.Q_flow] + sol[coolant.Q_flow] == zeros(length(sol))
end

# Test ConvectiveConductor, BodyRadiation
@testset "Radiator system" begin
    Tᵧ, Tᵪ = 1000, 10 # ᵧ -> gas and ᵪ -> coolant
    R_wall = 10
    G = 0.04
    σ   = 5.6703744191844294e-8 # Stefan-Boltzmann constant

    @named base        = ThermalResistor(R=R_wall)
    @named gas_tem     = FixedTemperature(T=Tᵧ)
    @named coolant_tem = FixedTemperature(T=Tᵪ)
    @named radiator    = BodyRadiation(G=G)
    @named dissipator  = ConvectiveConductor(G=10)
    @named mass        = HeatCapacitor(C=10)

    @info "Building a radiator..."
    eqs = [
        connect(gas_tem.port, radiator.port_a, base.port_a, dissipator.solid, mass.port)
        connect(coolant_tem.port, base.port_b, radiator.port_b, dissipator.fluid)
    ]
    @named rad = ODESystem(eqs, t, systems=[base, gas_tem, radiator, dissipator, coolant_tem, mass])
    sys = structural_simplify(rad)

    u0 = [
        base.Q_flow => 10
        dissipator.Q_flow => 10
        mass.T => Tᵧ
    ]
    prob = ODEProblem(sys, u0, (0, 3.0))
    sol = solve(prob, Rodas4())

    @test sol[dissipator.dT] == sol[radiator.port_a.T] - sol[radiator.port_b.T]
    rad_Q_flow = G*σ*(Tᵧ^4 - Tᵪ^4)
    @test sol[radiator.Q_flow] == fill(rad_Q_flow, length(sol[radiator.Q_flow]))
end

@testset "Thermal Collector" begin
    @named flow_src    = FixedHeatFlow(Q_flow=50, alpha=100)
    @named hf_sensor   = HeatFlowSensor()
    @named th_ground   = FixedTemperature(T=0)
    @named collector   = ThermalCollector(m=2)
    @named th_resistor = ThermalResistor(R=10)
    @named tem_src     = FixedTemperature(T=10)
    @named mass        = HeatCapacitor(C=10)

    @info "Building a heat collector..."
    eqs = [
        connect(flow_src.port, collector.port_a1, th_resistor.port_a)
        connect(tem_src.port, collector.port_a2)
        connect(hf_sensor.port_a, collector.port_b)
        connect(hf_sensor.port_b, mass.port, th_resistor.port_b)
        connect(mass.port, th_ground.port)
        ]
    @named coll = ODESystem(eqs, t,
                            systems=[hf_sensor,flow_src, tem_src,
                            collector, th_resistor, mass])
    sys = structural_simplify(coll)

    u0 = [
        th_resistor.Q_flow => 1.0,
        mass.T => 0.0,
    ]
    prob = ODEProblem(sys, u0, (0, 3.0))
    sol  = solve(prob, Rodas4())

    @test sol[collector.port_b.Q_flow] + sol[collector.port_a1.Q_flow] + sol[collector.port_a2.Q_flow] ==
        zeros(length(sol[collector.port_b.Q_flow]))
    @test sol[collector.port_b.T] == sol[collector.port_a1.T] == sol[collector.port_a2.T]
end


text:
mutable struct TestMDP <: MDP{Bool, Bool} end
mutable struct TestPOMDP <: POMDP{Bool, Bool, Bool} end

let
    # @test actions(TestMDP()) == (true, false)
    # @test actions(TestPOMDP()) == (true, false)

    a = [1,2,3]
    @test support(a) == a
    @test support((1,2,3)) == (1,2,3)

    # @test states(TestMDP()) == (true, false)
    # @test states(TestPOMDP()) == (true, false)

    # @test observations(TestPOMDP()) == (true, false)
    # @test n_observations(TestPOMDP()) == 2

    # @test stateindex(TestMDP(), 1) == 1
    # @test actionindex(TestMDP(), 2) == 2
    # @test obsindex(TestMDP(), 3) == 3
end


text:
using HDF5
using Test

@testset "Dataspaces" begin
    hsize_t = HDF5.hsize_t
    # Reference objects without using high-level API
    ds_null   = HDF5.Dataspace(HDF5.h5s_create(HDF5.H5S_NULL))
    ds_scalar = HDF5.Dataspace(HDF5.h5s_create(HDF5.H5S_SCALAR))
    ds_zerosz = HDF5.Dataspace(HDF5.h5s_create_simple(1, hsize_t[0], hsize_t[0]))
    ds_vector = HDF5.Dataspace(HDF5.h5s_create_simple(1, hsize_t[5], hsize_t[5]))
    ds_matrix = HDF5.Dataspace(HDF5.h5s_create_simple(2, hsize_t[7, 5], hsize_t[7, 5]))
    ds_maxdim = HDF5.Dataspace(HDF5.h5s_create_simple(2, hsize_t[7, 5], hsize_t[20, 20]))
    ds_unlim  = HDF5.Dataspace(HDF5.h5s_create_simple(1, hsize_t[1], [HDF5.H5S_UNLIMITED]))

    # Testing basic property accessors of dataspaces

    @test isvalid(ds_scalar)

    @test ndims(ds_null)   === 0
    @test ndims(ds_scalar) === 0
    @test ndims(ds_zerosz) === 1
    @test ndims(ds_vector) === 1
    @test ndims(ds_matrix) === 2

    # Test that properties of existing datasets can be extracted.
    # Note: Julia reverses the order of dimensions when using the high-level API versus
    #       the dimensions used above to create the reference objects.
    @test size(ds_null)   === ()
    @test size(ds_scalar) === ()
    @test size(ds_zerosz) === (0,)
    @test size(ds_vector) === (5,)
    @test size(ds_matrix) === (5, 7)
    @test size(ds_maxdim) === (5, 7)

    @test size(ds_null, 5)   === 1
    @test size(ds_scalar, 5) === 1
    @test size(ds_zerosz, 5) === 1
    @test size(ds_vector, 5) === 1
    @test size(ds_matrix, 5) === 1
    @test size(ds_maxdim, 5) === 1
    @test_throws ArgumentError("invalid dimension d; must be positive integer") size(ds_null, 0)
    @test_throws ArgumentError("invalid dimension d; must be positive integer") size(ds_scalar, -1)

    @test length(ds_null)   === 0
    @test length(ds_scalar) === 1
    @test length(ds_zerosz) === 0
    @test length(ds_vector) === 5
    @test length(ds_matrix) === 35
    @test length(ds_maxdim) === 35

    @test isempty(ds_null)
    @test !isempty(ds_scalar)
    @test isempty(ds_zerosz)
    @test !isempty(ds_vector)

    @test HDF5.isnull(ds_null)
    @test !HDF5.isnull(ds_scalar)
    @test !HDF5.isnull(ds_zerosz)
    @test !HDF5.isnull(ds_vector)

    @test HDF5.get_extent_dims(ds_null)   === ((), ())
    @test HDF5.get_extent_dims(ds_scalar) === ((), ())
    @test HDF5.get_extent_dims(ds_zerosz) === ((0,), (0,))
    @test HDF5.get_extent_dims(ds_vector) === ((5,), (5,))
    @test HDF5.get_extent_dims(ds_matrix) === ((5, 7), (5, 7))
    @test HDF5.get_extent_dims(ds_maxdim) === ((5, 7), (20, 20))
    @test HDF5.get_extent_dims(ds_unlim)  === ((1,), (-1,))

    # Can create new copies
    ds_tmp  = copy(ds_maxdim)
    ds_tmp2 = HDF5.Dataspace(ds_tmp.id) # copy of ID, but new Julia object
    @test ds_tmp.id == ds_tmp2.id != ds_maxdim.id
    # Equality and hashing
    @test ds_tmp  == ds_maxdim
    @test ds_tmp !== ds_maxdim
    @test hash(ds_tmp) != hash(ds_maxdim)
    @test ds_tmp  == ds_tmp2
    @test ds_tmp !== ds_tmp2
    @test hash(ds_tmp) == hash(ds_tmp2)

    # Behavior of closing dataspace objects
    close(ds_tmp)
    @test ds_tmp.id == -1
    @test !isvalid(ds_tmp)
    @test !isvalid(ds_tmp2)

    # Validity checking in high-level operations
    @test_throws ErrorException("File or object has been closed") copy(ds_tmp)
    @test_throws ErrorException("File or object has been closed") ndims(ds_tmp)
    @test_throws ErrorException("File or object has been closed") size(ds_tmp)
    @test_throws ErrorException("File or object has been closed") size(ds_tmp, 1)
    @test_throws ErrorException("File or object has been closed") length(ds_tmp)
    @test_throws ErrorException("File or object has been closed") ds_tmp == ds_tmp2
    @test close(ds_tmp) === nothing # no error

    # Test ability to create explicitly-sized dataspaces

    @test dataspace(()) == ds_scalar
    @test dataspace((5,)) == ds_vector
    @test dataspace((5, 7)) == ds_matrix != ds_maxdim
    @test dataspace((5, 7), max_dims = (20, 20)) == ds_maxdim != ds_matrix
    @test dataspace((1,), max_dims = (-1,)) == ds_unlim
    # for ≥ 2 numbers, same as single tuple argument
    @test dataspace(5, 7) == ds_matrix
    @test dataspace(5, 7, 1) == dataspace((5, 7, 1))

    # Test dataspaces derived from data

    @test dataspace(nothing) == ds_null
    @test dataspace(HDF5.EmptyArray{Bool}()) == ds_null

    @test dataspace(fill(1.0)) == ds_scalar
    @test dataspace(1) == ds_scalar
    @test dataspace(1 + 1im) == ds_scalar
    @test dataspace("string") == ds_scalar

    @test dataspace(zeros(0)) == ds_zerosz
    @test dataspace(zeros(0, 0)) != ds_zerosz
    @test dataspace(zeros(5, 7)) == ds_matrix
    @test dataspace(HDF5.VLen([[1]])) == dataspace((1,))
    @test dataspace(HDF5.VLen([[1], [2]])) == dataspace((2,))

    # Constructing dataspace for/from HDF5 dataset or attribute

    mktemp() do path, io
        close(io)
        h5open(path, "w") do hid
            dset = create_dataset(hid, "dset", datatype(Int), ds_matrix)
            attr = create_attribute(dset, "attr", datatype(Bool), ds_vector)
            @test dataspace(dset)  == ds_matrix
            @test dataspace(dset) !== ds_matrix
            @test dataspace(attr)  == ds_vector
            @test dataspace(attr) !== ds_vector
            close(dset)
            close(attr)
            @test_throws ErrorException("File or object has been closed") dataspace(dset)
            @test_throws ErrorException("File or object has been closed") dataspace(attr)
        end
    end

    # Test mid-level routines: set/get_extent_dims

    dspace_norm = dataspace((100, 4))
    @test HDF5.get_extent_dims(dspace_norm)[1] == HDF5.get_extent_dims(dspace_norm)[2] == (100, 4)
    HDF5.set_extent_dims(dspace_norm, (8, 2))
    @test HDF5.get_extent_dims(dspace_norm)[1] == HDF5.get_extent_dims(dspace_norm)[2] == (8, 2)

    dspace_maxd = dataspace((100, 4), max_dims = (256, 5))
    @test HDF5.get_extent_dims(dspace_maxd)[1] == (100, 4)
    @test HDF5.get_extent_dims(dspace_maxd)[2] == (256, 5)
    HDF5.set_extent_dims(dspace_maxd, (8, 2))
    @test HDF5.get_extent_dims(dspace_maxd)[1] == (8, 2)
    HDF5.set_extent_dims(dspace_maxd, (3, 1), (4,2))
    @test HDF5.get_extent_dims(dspace_maxd)[1] == (3, 1)
    @test HDF5.get_extent_dims(dspace_maxd)[2] == (4, 2)
    HDF5.set_extent_dims(dspace_maxd, (3, 1), (-1, -1)) # unlimited max size
    @test HDF5.get_extent_dims(dspace_maxd)[1] == (3, 1)
    @test HDF5.get_extent_dims(dspace_maxd)[2] == (-1, -1)
end


text:
module PolynomialsBenchs

include("MonomialBasesBenchs.jl")

end # module


text:
"""
Generate the `src/services/{service}.jl` file.
"""
function _generate_high_level_wrapper(services::AbstractArray{<:AbstractDict}, repo_name::String, auth::GitHub.OAuth2)
    for service in services
        service_name = service["path"]
        @info "Generating high level wrapper for $service_name"

        service_blob = blob(repo_name, service["sha"]; auth=auth)
        service = JSON.parse(String(base64decode(service_blob.content)))

        service_name = lowercase(service["metadata"]["serviceId"])
        service_name = replace(service_name, ' ' => '_')
        operations = service["operations"]
        shapes = service["shapes"]

        protocol = service["metadata"]["protocol"]

        operations = sort!(_generate_high_level_definitions(service_name, protocol, operations, shapes))

        service_path = joinpath(@__DIR__, "..", "services", "$service_name.jl")
        open(service_path, "w") do f
            println(f, """
                    # This file is auto-generated by AWSMetadata.jl
                    using AWS
                    using AWS.AWSServices: $service_name
                    using AWS.Compat
                    using AWS.UUIDs
                    """)
            join(f, operations, "\n")
        end
    end
end

"""
Generate high-level definitions for the `service`.
All high-level definitions and documentation to be written into `services/{Service}.jl`
"""
function _generate_high_level_definitions(
    service_name::String,
    protocol::String,
    operations::AbstractDict,
    shapes::AbstractDict
)
    operation_definitions = String[]

    for (_, operation) in operations
        name = operation["name"]
        method = operation["http"]["method"]
        request_uri = operation["http"]["requestUri"]

        documentation = ""

        if haskey(operation, "documentation")
            documentation = _clean_documentation(operation["documentation"])
        end

        required_parameters = Dict{String, Any}()
        optional_parameters = Dict{String, Any}()

        if haskey(operation, "input")
            required_parameters, optional_parameters = _get_function_parameters(operation["input"]["shape"], shapes)
        end

        operation_definition = _generate_high_level_definition(
            service_name,
            protocol,
            name,
            method,
            request_uri,
            required_parameters,
            optional_parameters,
            documentation
        )

        push!(operation_definitions, operation_definition)
    end

    return operation_definitions
end

"""
Generate the high-level definition for a services function.
"""
function _generate_high_level_definition(
    service_name::String,
    protocol::String,
    name::String,
    method::String,
    request_uri::String,
    required_parameters::AbstractDict,
    optional_parameters::AbstractDict,
    documentation::String
)
    """
    Generate function definition for a service request given required, header and idempotent parameters.
    """
    function _generate_rest_operation_defintion(required_params::AbstractDict, optional_params::AbstractDict, function_name::String, service_name::String, method::String, request_uri::String)
        request_uri = replace(request_uri, '{' => "\$(")  # Replace { with $(
        request_uri = replace(request_uri, '}' => ')')  # Replace } with )
        request_uri = replace(request_uri, '+' => "")  # Remove + from the request URI

        # Pre Julia-1.3 workaround
        req_keys = [replace(key, "-"=>"_") for key in collect(keys(required_params))]

        required_params = filter(p -> (p[2]["location"] != "uri"), required_params)
        header_params = filter(p -> (p[2]["location"] == "header"), required_params)
        required_params = setdiff(required_params, header_params)

        idempotent_params = filter(p -> (p[2]["idempotent"]), optional_params)

        req_kv = ["\"$(p[1])\"=>$(replace(p[1], "-" => "_"))" for p in required_params]
        header_kv = ["\"$(p[1])\"=>$(replace(p[1], "-" => "_"))" for p in header_params]
        idempotent_kv = ["\"$(p[1])\"=>string(uuid4())" for p in idempotent_params]

        required_keys = !isempty(req_keys)
        headers = !isempty(header_params)
        idempotent = !isempty(idempotent_params)

        req_str = !isempty(req_kv) ? "Dict{String, Any}($(join(req_kv, ", "))" : ""
        params_str = (!isempty(req_kv) || idempotent) ? "$(join(vcat(req_kv, idempotent_kv), ", "))" : ""
        headers_str = headers ? "\"headers\"=>Dict{String, Any}($(join(header_kv, ", ")))" : ""
        params_headers_str = "Dict{String, Any}($(join([s for s in (params_str, headers_str) if !isempty(s)], ", ")))"

        formatted_function_name = _format_name(function_name)

        if required_keys && (idempotent || headers)
            return """
                $formatted_function_name($(join(req_keys, ", ")); aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", $params_headers_str; aws_config=aws_config)
                $formatted_function_name($(join(req_keys, ", ")), params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", Dict{String, Any}(mergewith(_merge, $params_headers_str, params)); aws_config=aws_config)
                """
        elseif !required_keys && (idempotent || headers)
            return """
                $formatted_function_name(; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", $params_headers_str; aws_config=aws_config)
                $formatted_function_name(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", Dict{String, Any}(mergewith(_merge, $params_headers_str, params)); aws_config=aws_config)
                """
        elseif required_keys && !isempty(req_kv)
            return """
                $formatted_function_name($(join(req_keys, ", ")); aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", $req_str); aws_config=aws_config)
                $formatted_function_name($(join(req_keys, ", ")), params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", Dict{String, Any}(mergewith(_merge, $req_str), params)); aws_config=aws_config)
                """
        elseif required_keys
            return """
                $formatted_function_name($(join(req_keys, ", ")); aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\"; aws_config=aws_config)
                $formatted_function_name($(join(req_keys, ", ")), params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", params; aws_config=aws_config)
                """
        else
            return """
                $formatted_function_name(; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\"; aws_config=aws_config)
                $formatted_function_name(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", params; aws_config=aws_config)
                """
        end
    end

    """
    Generate a JSON/Query high level definition.
    """
    function _generate_json_query_opeation_definition(required_params::AbstractDict, optional_params::AbstractDict, function_name::String, service_name::String)
        req_keys = [replace(key, '-' => '_') for key in collect(keys(required_params))]

        idempotent_params = filter(p -> (p[2]["idempotent"]), optional_params)

        req_kv = ["\"$(p[1])\"=>$(replace(p[1], "-" => "_"))" for p in required_params]
        idempotent_kv = ["\"$(p[1])\"=>string(uuid4())" for p in idempotent_params]

        required = !isempty(req_kv)
        idempotent = !isempty(idempotent_kv)

        formatted_function_name = _format_name(function_name)

        if required && idempotent
            return """
                $formatted_function_name($(join(req_keys, ", ")); aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}($(join(req_kv, ", ")), $(join(idempotent_kv, ", "))); aws_config=aws_config)
                $formatted_function_name($(join(req_keys, ", ")), params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}(mergewith(_merge, Dict{String, Any}($(join(req_kv, ", ")), $(join(idempotent_kv, ", "))), params)); aws_config=aws_config)
                """
        elseif required
            return """
                $formatted_function_name($(join(req_keys, ", ")); aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}($(join(req_kv, ", "))); aws_config=aws_config)
                $formatted_function_name($(join(req_keys, ", ")), params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}(mergewith(_merge, Dict{String, Any}($(join(req_kv, ", "))), params)); aws_config=aws_config)
                """
        elseif idempotent
            return """
                $formatted_function_name(; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}($(join(idempotent_kv, ", "))); aws_config=aws_config)
                $formatted_function_name(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}(mergewith(_merge, Dict{String, Any}($(join(idempotent_kv, ", "))), params)); aws_config=aws_config)
                """
        else
            return """
                $formatted_function_name(; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$function_name\"; aws_config=aws_config)
                $formatted_function_name(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$function_name\", params; aws_config=aws_config)
                """
        end
    end

    """
    Generate the docstring for the `function_name`.
    """
    function _generate_docstring(function_name, documentation, required_parameters, optional_parameters)
        function_name = _format_name(function_name)
        args = join((_format_name(key) for (key, val) in required_parameters), ", ")
        maybejoin = isempty(args) ? "" : ", "
        operation_definition = """
            $(repeat('"', 3))
                $function_name($(args))
                $function_name($(args)$(maybejoin)params::Dict{String,<:Any})

            $(_wraplines(documentation))\n
            """

        # Add in the required parameters if applicable
        if !isempty(required_parameters)
            operation_definition *= "# Arguments\n"

            for (required_key, required_value) in required_parameters
                key = _format_name(required_key)
                operation_definition *= _wraplines("- `$key`: $(required_value["documentation"])"; delim="\n  ")
                operation_definition *= "\n"
            end

            operation_definition *= "\n"
        end

        # Add in the optional parameters if applicable
        if !isempty(optional_parameters)
            operation_definition *= """
                # Optional Parameters
                Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
                """

            for (optional_key, optional_value) in optional_parameters
                operation_definition *= _wraplines("- `\"$optional_key\"`: $(optional_value["documentation"])"; delim="\n  ")
                operation_definition *= "\n"
            end
        end

        operation_definition *= repeat('"', 3)
    end

    doc_string = _generate_docstring(name, documentation, required_parameters, optional_parameters)

    if protocol in ("json", "query", "ec2")
        function_string = _generate_json_query_opeation_definition(required_parameters, optional_parameters, name, service_name)
    elseif protocol in ("rest-json", "rest-xml")
        function_string = _generate_rest_operation_defintion(required_parameters, optional_parameters, name, service_name, method, request_uri)
    else
        throw(ProtocolNotDefined("$function_name is using a new protocol; $protocol which is not supported."))
    end

    return string(doc_string, '\n', function_string)
end


text:
using PyPlot, PyCall

np = pyimport("numpy")

origin = "lower"

delta = 0.025

x = y = np.arange(-3.0, 3.01, delta)
X, Y = np.meshgrid(x, y)
Z1 = @. np.exp(-X^2 - Y^2)
Z2 = @. np.exp(-(X - 1)^2 - (Y - 1)^2)
Z = (Z1 .- Z2) .* 2

nr, nc = size(Z)

# mask a circle in the middle:
interior = @. sqrt(X^2 + Y^2) < 0.5
Z[interior] .= NaN

# We are using automatic selection of contour levels;
# this is usually not such a good idea, because they don't
# occur on nice boundaries, but we do it here for purposes
# of illustration.

fig1, ax2 = plt.subplots(constrained_layout=true)
CS = ax2.contourf(X, Y, Z, 10, cmap=plt.cm.bone, origin=origin)

# Note that in the following, we explicitly pass in a subset of
# the contour levels used for the filled contours.  Alternatively,
# We could pass in additional levels to provide extra resolution,
# or leave out the levels kwarg to use all of the original levels.

CS2 = ax2.contour(CS, levels=CS.levels[1:2:end], colors="r", origin=origin)

ax2.set_title("1 masked region")
ax2.set_xlabel("word length anomaly")
ax2.set_ylabel("sentence length anomaly")

# Make a colorbar for the ContourSet returned by the contourf call.
cbar = fig1.colorbar(CS)
cbar.ax.set_ylabel("verbosity coefficient")
# Add the contour line levels to the colorbar
cbar.add_lines(CS2)

plt.show()


text:
type LocalQuantities{T <: AbstractFloat} <: AbstractLocalQuantities
    density::T
    momentum::DenseVector{T}
    velocity::DenseVector{T}
    feq::DenseVector{T}
end

function LocalQuantities{T}(fᵢ::DenseVector{T}, lattice::Lattice{T})
    const ρ = density(fᵢ)
    const μ = momentum(fᵢ, lattice.celerities)
    const ν = velocity(μ, ρ)
    const feq = equilibrium(ρ, μ, lattice.celerities, lattice.weights)
    LocalQuantities{T}(ρ, μ, ν, feq)
end

# Pre-allocates a local quantity aggregator
function LocalQuantities(lattice::Lattice)
    const T = eltype(lattice.weights)

    LocalQuantities{T}(
        convert(T, 0), zeros(T, ndims(lattice)),
        zeros(T, ndims(lattice)), zeros(T, length(lattice))
    )
end

# Recomputes local-quantities in-place
function LocalQuantities!{T}(out::LocalQuantities{T}, fᵢ::DenseVector{T}, lattice::Lattice{T})
    out.density = density(fᵢ)
    momentum!(out.momentum, fᵢ, lattice.celerities)
    velocity!(out.velocity, out.momentum, out.density)
    equilibrium!(out.feq, out.density, out.momentum, lattice.celerities, lattice.weights)
    out
end


text:
using ClickHouse

@testset "server capabilites" begin

    @test ClickHouse.has_temporary_tables(50264)
    @test ClickHouse.has_temporary_tables(50274)
    @test !ClickHouse.has_temporary_tables(50263)
end


text:
#TODO: OpenCV functions
mutable struct Distortion
    camera_matrix::Matrix
    distortion_coefs::Matrix
    rectification_matrix::Matrix
    projection_matrix::Matrix
    mapx
    mapy
    rmapx
    rmapy
end

function Distortion()

    camera_matrix = [
        305.5718893575089 0 303.0797142544728;
        0 308.8338858195428 231.8845403702499;
        0 0 1
    ]

    # distortion parameters - (k1, k2, t1, t2, k3)
    distortion_coefs = [
        -0.2, 0.0305,
        0.0005859930422629722, -0.0006697840226199427, 0
    ]

    distortion_coefs = reshape(distortion_coefs, (1, 5))

    # R - Rectification matrix - stereo cameras only, so identity
    rectification_matrix = eye(3)

    # P - Projection Matrix - specifies the intrinsic (camera) matrix
    #  of the processed (rectified) image
    projection_matrix = [
        220.2460277141687, 0, 301.8668918355899, 0;
        0, 238.6758484095299, 227.0880056118307, 0;
        0, 0, 1, 0
    ]

    # Initialize mappings

    # Used for rectification
    mapx = nothing
    mapy = nothing

    # Used for distortion
    rmapx = nothing
    rmapy = nothing

    Distortion(camera_matrix, distortion_coefs, rectification_matrix,
               projection_matrix, mapx, mapy, rmapx, rmapy)
end

function distort(dist::Distorion, observation)
    ##
    #Distort observation using parameters in constructor
    ##

    if isnothing(dist.mapx)
        # Not initialized - initialize all the transformations we'll need
        dist.mapx = zeros(size(observation))
        dist.mapy = zeros(size(observation))

        H, W, _ = size(observation)

        # Initialize dist.mapx and dist.mapy (updated)
        dist.mapx, dist.mapy = cv2.initUndistortRectifyMap(dist.camera_matrix,
                                                           dist.distortion_coefs, dist.rectification_matrix,
                                                           dist.projection_matrix, (W, H), cv2.CV_32FC1)

        # Invert the transformations for the distortion
        dist.rmapx, dist.rmapy = _invert_map(dist, dist.mapx, dist.mapy)
    end

    return cv2.remap(observation, self.rmapx, self.rmapy, interpolation=cv2.INTER_NEAREST)
end

function _undistort(dist::Distortion, observation)
    ##
    #Undistorts a distorted image using camera parameters
    ##

    # If mapx is None, then distort was never called
    @assert !isnothing(dist.mapx) "You cannot call undistort on a rectified image"

    return cv2.remap(observation, self.mapx, self.mapy, cv2.INTER_NEAREST)
end

function _invert_map(self, mapx, mapy):
    ##
    #Utility function for simulating distortion
    #Source: https://github.com/duckietown/Software/blob/master18/catkin_ws
    #... /src/10-lane-control/ground_projection/include/ground_projection/
    #... ground_projection_geometry.py
    ##

    H, W, _ = size(mapx)
    rmapx = similar(mapx)
    fill!(rmapx, NaN)
    rmapy = similar(mapx)
    fill!(rmapy, NaN)

    for y in 1:H, x in 1:W
        tx = mapx[y, x]
        ty = mapy[y, x]

        tx = Int(round(tx))
        ty = int(round(ty))

        if (1 ≤ tx ≤ W) && (1 ≤ ty ≤ H)
            rmapx[ty, tx] = x
            rmapy[ty, tx] = y
        end
    end
    _fill_holes(dist, rmapx, rmapy)
    return rmapx, rmapy
end

function _fill_holes(dist::Distortion, rmapx, rmapy)
    ##
    #Utility function for simulating distortion
    #Source: https://github.com/duckietown/Software/blob/master18/catkin_ws
    #... /src/10-lane-control/ground_projection/include/ground_projection/
    #... ground_projection_geometry.py
    ##
    H, W, _ = size(rmapx)

    R = 2
    F = R * 2 + 1

    norm(_) = hypot(_[1], _[2])

    deltas0 = [(i - R - 1, j - R - 1) for i in 1:F, j in 1:F]
    deltas0 = [x for x in deltas0 if norm(x) ≤ R]
    sort(deltas0, by=norm)

    holes = Set()

    for i = 1:H, j in 1:W
        if isnan(rmapx[i, j])
            union!(holes, [(i, j)])
        end
    end

    while !isempty(holes)
        nholes = length(holes)
        nholes_filled = 0

        for (i, j) in holes
            # there is nan
            nholes += 1
            for di, dj in deltas0
                u = i + di
                v = j + dj
                if (1 ≤ u ≤ H) && (1 ≤ v ≤ W)
                    if !isnan(rmapx[u, v])
                        rmapx[i, j] = rmapx[u, v]
                        rmapy[i, j] = rmapy[u, v]
                        nholes_filled += 1
                        pop!(holes, (i, j))
                        break
                    end
                end
            end
        end

        nholes_filled == 0 && break
    end
end


text:
#__precompile__()
module GUIFrameworkPlayground

    using SDL2, Cairo, ColorTypes, Graphics, DataStructures
    using GLVisualize

    import Cairo: rectangle
    import Graphics: width, height
    import GLVisualize: GLWindow, GLFW, ModernGL, GLAbstraction, GeometryTypes

    abstract type AbstractGUIBackend end

    window(b::AbstractGUIBackend) = error("This method must be implemented.")
    init!(b::AbstractGUIBackend) = error("This method must be implemented.")

    abstract type AbstractGUIWindow end

    info("loading events")
    include("events.jl")

    include("backends/SDL2GUIBackend.jl")
    include("backends/OpenGLGUIBackend.jl")

end # module

##

##


text:
using BattleSim
using Test

include("testutil.jl")

@testset "All tests" begin

@testset "Arrow depletion" begin
    archer = Archer("Robin Hood", 28)

    arrows = archer.arrows
    resupply!(archer)
    @test arrows == archer.arrows

    shoot!(archer)

    @test arrows > archer.arrows

    arrows = archer.arrows
    resupply!(archer)

    @test archer.arrows > arrows
end

@testset "Archer attack" begin
    archer = Archer("Robin Hood", 28)
    piky = Pikeman("Swiss Guard", 32)

    archer_health = archer.health
    piky_health = piky.health

    for _ in 1:4
        attack!(archer, piky)
    end

    @test archer_health == archer.health
    @test piky_health > piky.health
end

@testset "Archer vs Knight" begin

    # Since this is not deterministic we have to rely on the fact that
    # if this is coded correctly knights should statistically speaking
    # win more often.

    wins, losses =  wins_and_losses(Knight, Archer)
    @test wins > losses

    wins, losses =  wins_and_losses(Archer, Knight)
    @test wins < losses
end

@testset "Pikman vs Knight" begin
    wins, losses =  wins_and_losses(Knight, Pikeman)
    @test wins < losses

    wins, losses =  wins_and_losses(Pikeman, Knight)
    @test wins > losses
end

@testset "Archer vs Pikeman" begin
    wins, losses =  wins_and_losses(Archer, Pikeman)
    @test wins > losses

    wins, losses =  wins_and_losses(Pikeman, Archer)
    @test wins < losses
end

end

max_stars_repo_name: ordovician/BattleSim

max_issues_repo_name: ordovician/BattleSim

max_forks_repo_name: ordovician/BattleSim


text:
# Generated functions for tensor-product polynomials where nodes are in an array of arrays

@generated function chebyshev_weights(f::AbstractArray{T,N},nodes::Array{Array{T,1},1},order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T,N,S}

  chebyshev_polynomials = :( poly = Array{Array{T,2},1}(undef,N);
                             for k = 1:N;
                               orderk = order[k];
                               xk = nodes[k];

                               # normalize nodes

                               if domain[1,k] == domain[2,k];
                                 fill!(xk,(domain[1,k].+domain[2,k])/2);
                               else;
                                 xk = 2*(xk.-domain[2,k])/(domain[1,k].-domain[2,k]).-one(T);
                               end;

                               polynomial = Array{T}(undef,length(xk),orderk+1);
                               for i = 1:orderk+1;
                                 for j = 1:length(xk);
                                   if i == 1;
                                     polynomial[j,i] = one(T);
                                   elseif i == 2;
                                     polynomial[j,i] = xk[j];
                                   else;
                                     polynomial[j,i] = 2*xk[j]*polynomial[j,i-1]-polynomial[j,i-2];
                                   end;
                                 end;
                               end;
                               poly[i] = polynomial;#push!(poly,polynomial);
                             end )

  i_vars = Array{Symbol}(undef,N)
  s_vars = Array{Symbol}(undef,N)
  for i = 1:N
    i_vars[i] = Symbol("i$i")
    s_vars[i] = Symbol("s$i")
  end

  # Construct denominator term

  inner_prod_denominator = :( poly[1][s1,i1] )
  for i = 2:N
    inner_prod_denominator = :( poly[$i][$(s_vars[i]),$(i_vars[i])]*$inner_prod_denominator )
  end
  denominator_term = :( $inner_prod_denominator*$inner_prod_denominator )

  # Construct numerator term

  inner_prod_numerator = :( poly[1][s1,i1]*f[$(s_vars...)] )
  for i = 2:N
      inner_prod_numerator = :( poly[$i][$(s_vars[i]),$(i_vars[i])]*$inner_prod_numerator )
  end
  numerator_term = inner_prod_numerator

  # Construct weights term

  weights_term = :( weights[$(i_vars...)] = numerator/denominator )

  # Construct the inner loops that compute numerator and denominator

  inner = :( numerator += $numerator_term;
             denominator += $denominator_term )
  outer  = inner

  for i = N:-1:1#1:N
    outer = :(
      for $(s_vars[i]) = 1:size(f,$i)
        $outer
      end
    )
  end

  # Construct the outer loops that compute the weights

  new_inner = :( numerator = zero(T);
                 denominator = zero(T);
                 $outer;
                 $weights_term )
  new_outer = new_inner

  for i = N:-1:1#1:N
    new_outer = :(
      for $(i_vars[i]) = 1:(order[$i]+1)
        $new_outer
      end
    )
  end

  # Put it all together to compute the weights array

  final = :( $chebyshev_polynomials;
             weights = Array{T}(undef,(order+1)...);
             $new_outer;
             return weights )

  return final

end

@generated function chebyshev_weights(f::AbstractArray{T,N},poly::Array{Array{T,2},1},order::Array{S,1}) where {T,N,S}

  i_vars = Array{Symbol}(undef,N)
  s_vars = Array{Symbol}(undef,N)
  for i = 1:N
    i_vars[i] = Symbol("i$i")
    s_vars[i] = Symbol("s$i")
  end

  # Construct denominator term

  inner_prod_denominator = :( poly[1][s1,i1] )
  for i = 2:N
    inner_prod_denominator = :( poly[$i][$(s_vars[i]),$(i_vars[i])]*$inner_prod_denominator )
  end
  denominator_term = :( $inner_prod_denominator*$inner_prod_denominator )

  # Construct numerator term

  inner_prod_numerator = :( poly[1][s1,i1]*f[$(s_vars...)] )
  for i = 2:N
      inner_prod_numerator = :( poly[$i][$(s_vars[i]),$(i_vars[i])]*$inner_prod_numerator )
  end
  numerator_term = inner_prod_numerator

  # Construct weights term

  weights_term = :( weights[$(i_vars...)] = numerator/denominator )

  # Construct the inner loops that compute numerator and denominator

  inner = :( numerator += $numerator_term;
             denominator += $denominator_term )
  outer  = inner

  for i = N:-1:1#1:N
    outer = :(
      for $(s_vars[i]) = 1:size(f,$i)
        $outer
      end
    )
  end

  # Construct the outer loops that compute the weights

  new_inner = :( numerator = zero(T);
                 denominator = zero(T);
                 $outer;
                 $weights_term )
  new_outer = new_inner

  for i = N:-1:1#1:N
    new_outer = :(
      for $(i_vars[i]) = 1:(order[$i]+1)
        $new_outer
      end
    )
  end

    # Put it all together to compute the weights array

  final = :( weights = Array{T}(undef,(order.+1)...);
             $new_outer;
             return weights )

  return final

end

# Generated functions for tensor-product polynomials where nodes are in a tuple

@generated function chebyshev_weights(f::AbstractArray{T,N},nodes::NTuple{N,Array{T,1}},order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T,N,S}

  chebyshev_polynomials = :( poly = Array{Array{T,2},1}(undef,N);
                             for k = 1:N;
                               orderk = order[k];
                               xk = nodes[k];

                               # normalize nodes

                               if domain[1,k] == domain[2,k];
                                 fill!(xk,(domain[1,k].+domain[2,k])/2);
                               else;
                                 xk = 2*(xk.-domain[2,k])/(domain[1,k].-domain[2,k]).-one(T);
                               end;

                               polynomial = Array{T}(undef,length(xk),orderk+1);
                               for i = 1:orderk+1;
                                 for j = 1:length(xk);
                                   if i == 1;
                                     polynomial[j,i] = one(T);
                                   elseif i == 2;
                                     polynomial[j,i] = xk[j];
                                   else;
                                     polynomial[j,i] = 2*xk[j]*polynomial[j,i-1]-polynomial[j,i-2];
                                   end;
                                 end;
                               end;
                               poly[k] = polynomial;#push!(poly,polynomial);
                             end )

  i_vars = Array{Symbol}(undef,N)
  s_vars = Array{Symbol}(undef,N)
  for i = 1:N
    i_vars[i] = Symbol("i$i")
    s_vars[i] = Symbol("s$i")
  end

  # Construct denominator term

  inner_prod_denominator = :( poly[1][s1,i1] )
  for i = 2:N
    inner_prod_denominator = :( poly[$i][$(s_vars[i]),$(i_vars[i])]*$inner_prod_denominator )
  end
  denominator_term = :( $inner_prod_denominator*$inner_prod_denominator )

  # Construct numerator term

  inner_prod_numerator = :( poly[1][s1,i1]*f[$(s_vars...)] )
  for i = 2:N
      inner_prod_numerator = :( poly[$i][$(s_vars[i]),$(i_vars[i])]*$inner_prod_numerator )
  end
  numerator_term = inner_prod_numerator

  # Construct weights term

  weights_term = :( weights[$(i_vars...)] = numerator/denominator )

  # Construct the inner loops that compute numerator and denominator

  inner = :( numerator += $numerator_term;
             denominator += $denominator_term )
  outer  = inner

  for i = N:-1:1#1:N
    outer = :(
      for $(s_vars[i]) = 1:size(f,$i)
        $outer
      end
    )
  end

  # Construct the outer loops that compute the weights

  new_inner = :( numerator = zero(T);
                 denominator = zero(T);
                 $outer;
                 $weights_term )
  new_outer = new_inner

  for i = N:-1:1#1:N
    new_outer = :(
      for $(i_vars[i]) = 1:(order[$i]+1)
        $new_outer
      end
    )
  end

  # Put it all together to compute the weights array

  final = :( $chebyshev_polynomials;
             weights = Array{T}(undef,(order.+1)...);
             $new_outer;
             return weights )

  return final

end

@generated function chebyshev_weights(f::AbstractArray{T,N},poly::NTuple{N,Array{T,2}},order::Array{S,1}) where {T,N,S}

  i_vars = Array{Symbol}(undef,N)
  s_vars = Array{Symbol}(undef,N)
  for i = 1:N
    i_vars[i] = Symbol("i$i")
    s_vars[i] = Symbol("s$i")
  end

  # Construct denominator term

  inner_prod_denominator = :( poly[1][s1,i1] )
  for i = 2:N
    inner_prod_denominator = :( poly[$i][$(s_vars[i]),$(i_vars[i])]*$inner_prod_denominator )
  end
  denominator_term = :( $inner_prod_denominator*$inner_prod_denominator )

  # Construct numerator term

  inner_prod_numerator = :( poly[1][s1,i1]*f[$(s_vars...)] )
  for i = 2:N
      inner_prod_numerator = :( poly[$i][$(s_vars[i]),$(i_vars[i])]*$inner_prod_numerator )
  end
  numerator_term = inner_prod_numerator

  # Construct weights term

  weights_term = :( weights[$(i_vars...)] = numerator/denominator )

  # Construct the inner loops that compute numerator and denominator

  inner = :( numerator += $numerator_term;
             denominator += $denominator_term )
  outer  = inner

  for i = N:-1:1#1:N
    outer = :(
      for $(s_vars[i]) = 1:size(f,$i)
        $outer
      end
    )
  end

  # Construct the outer loops that compute the weights

  new_inner = :( numerator = zero(T);
                 denominator = zero(T);
                 $outer;
                 $weights_term )
  new_outer = new_inner

  for i = N:-1:1#1:N
    new_outer = :(
      for $(i_vars[i]) = 1:(order[$i]+1)
        $new_outer
      end
    )
  end

  # Put it all together to compute the weights array

  final = :( weights = Array{T}(undef,(order.+1)...);
             $new_outer;
             return weights )

  return final

end

# Generated functions for complete polynomials where nodes are in an array of arrays

@generated function chebyshev_weights(f::AbstractArray{T,N},nodes::Array{Array{T,1},1},order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T,N,S}

  chebyshev_polynomials = :( poly = Array{Array{T,2},1}(undef,N);
                             for k = 1:N;
                               xk = nodes[k];

                               # normalize nodes

                               if domain[1,k] == domain[2,k];
                                 fill!(xk,(domain[1,k].+domain[2,k])/2);
                               else;
                                 xk = 2*(xk.-domain[2,k])/(domain[1,k].-domain[2,k]).-one(T);
                               end;

                               polynomial = Array{T}(undef,length(xk),order+1);
                               for i = 1:order+1;
                                 for j = 1:length(xk);
                                   if i == 1;
                                     polynomial[j,i] = one(T);
                                   elseif i == 2;
                                     polynomial[j,i] = xk[j];
                                   else;
                                     polynomial[j,i] = 2*xk[j]*polynomial[j,i-1]-polynomial[j,i-2];
                                   end;
                                 end;
                               end;
                               poly[k] = polynomial;#push!(poly,polynomial);
                             end )

  i_vars = Array{Symbol}(undef,N)
  s_vars = Array{Symbol}(undef,N)
  for i = 1:N
    i_vars[i] = Symbol("i$i")
    s_vars[i] = Symbol("s$i")
  end

  # Construct denominator term

  inner_prod_denominator = :( poly[1][s1,i1] )
  for i = 2:N
    inner_prod_denominator = :( poly[$i][$(s_vars[i]),$(i_vars[i])]*$inner_prod_denominator )
  end
  denominator_term = :( $inner_prod_denominator*$inner_prod_denominator )

  # Construct numerator term

  inner_prod_numerator = :( poly[1][s1,i1]*f[$(s_vars...)] )
  for i = 2:N
      inner_prod_numerator = :( poly[$i][$(s_vars[i]),$(i_vars[i])]*$inner_prod_numerator )
  end
  numerator_term = inner_prod_numerator

  # Construct weights term

  weights_term = :( weights[$(i_vars...)] = numerator/denominator )

  # Construct the inner loops that compute numerator and denominator

  inner = :( numerator += $numerator_term;
             denominator += $denominator_term )
  outer  = inner

  for i = N:-1:1#1:N
    outer = :(
      for $(s_vars[i]) = 1:size(f,$i)
        $outer
      end
    )
  end

  # Construct the outer loops that compute the weights

  new_inner = :( numerator = zero(T);
                 denominator = zero(T);
                 if sum(tuple($(i_vars...))) <= order+N;
                   $outer;
                   $weights_term;
                 end )
  new_outer = new_inner

  for i = N:-1:1#1:N
    new_outer = :(
      for $(i_vars[i]) = 1:(order+1)
        $new_outer
      end
    )
  end

  # Initialize the weight Array

#=

  initial_weight = string("weights = zeros(order+1,",)
  for i = 2:N
    initial_weight = string(initial_weight,"order+1,")
  end
  initial_weight = parse(string(initial_weight,")"))

=#

  # Put it all together to compute the weights array

  final = :( $chebyshev_polynomials;
             weights = Array{T}(undef,(order.+1)...);
#             $initial_weight;
             $new_outer;
             return weights )

  return final

end

@generated function chebyshev_weights(f::AbstractArray{T,N},poly::Array{Array{T,2},1},order::S) where {T,N,S}

  i_vars = Array{Symbol}(undef,N)
  s_vars = Array{Symbol}(undef,N)
  for i = 1:N
    i_vars[i] = Symbol("i$i")
    s_vars[i] = Symbol("s$i")
  end

  # Construct denominator term

  inner_prod_denominator = :( poly[1][s1,i1] )
  for i = 2:N
    inner_prod_denominator = :( poly[$i][$(s_vars[i]),$(i_vars[i])]*$inner_prod_denominator )
  end
  denominator_term = :( $inner_prod_denominator*$inner_prod_denominator )

  # Construct numerator term

  inner_prod_numerator = :( poly[1][s1,i1]*f[$(s_vars...)] )
  for i = 2:N
      inner_prod_numerator = :( poly[$i][$(s_vars[i]),$(i_vars[i])]*$inner_prod_numerator )
  end
  numerator_term = inner_prod_numerator

  # Construct weights term

  weights_term = :( weights[$(i_vars...)] = numerator/denominator )

  # Construct the inner loops that compute numerator and denominator

  inner = :( numerator += $numerator_term;
             denominator += $denominator_term )
  outer  = inner

  for i = N:-1:1#1:N
    outer = :(
      for $(s_vars[i]) = 1:size(f,$i)
        $outer
      end
    )
  end

  # Construct the outer loops that compute the weights

  new_inner = :( numerator = zero(T);
                 denominator = zero(T);
                 if sum(tuple($(i_vars...))) <= order+N;
                   $outer;
                   $weights_term;
                 end )
  new_outer = new_inner

  for i = N:-1:1#1:N
    new_outer = :(
      for $(i_vars[i]) = 1:(order+1)
        $new_outer
      end
    )
  end

  # Initialize the weight Array

#=

  initial_weight = string("weights = zeros(order+1,",)
  for i = 2:N
    initial_weight = string(initial_weight,"order+1,")
  end
  initial_weight = parse(string(initial_weight,")"))

=#

  # Put it all together to compute the weights array

  final = :( weights = Array{T}(undef,(order.+1)...);
#             $initial_weight;
             $new_outer;
             return weights )

  return final

end

# Generated functions for complete polynomials where nodes are in a tuple

@generated function chebyshev_weights(f::AbstractArray{T,N},nodes::NTuple{N,Array{T,1}},order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T,N,S}

  chebyshev_polynomials = :( poly = Array{Array{T,2},1}(undef,N);
                             for k = 1:N;
                               xk = nodes[k];

                               # normalize nodes

                               if domain[1,k] == domain[2,k];
                                 fill!(xk,(domain[1,k].+domain[2,k])/2);
                               else;
                                 xk = 2*(xk.-domain[2,k])/(domain[1,k].-domain[2,k]).-one(T);
                               end;

                               polynomial = Array{T}(undef,length(xk),order+1);
                               for i = 1:order+1;
                                 for j = 1:length(xk);
                                   if i == 1;
                                     polynomial[j,i] = one(T);
                                   elseif i == 2;
                                     polynomial[j,i] = xk[j];
                                   else;
                                     polynomial[j,i] = 2*xk[j]*polynomial[j,i-1]-polynomial[j,i-2];
                                   end;
                                 end;
                               end;
                               poly[k] = polynomial;#push!(poly,polynomial);
                             end )

  i_vars = Array{Symbol}(undef,N)
  s_vars = Array{Symbol}(undef,N)
  for i = 1:N
    i_vars[i] = Symbol("i$i")
    s_vars[i] = Symbol("s$i")
  end

  # Construct denominator term

  inner_prod_denominator = :( poly[1][s1,i1] )
  for i = 2:N
    inner_prod_denominator = :( poly[$i][$(s_vars[i]),$(i_vars[i])]*$inner_prod_denominator )
  end
  denominator_term = :( $inner_prod_denominator*$inner_prod_denominator )

  # Construct numerator term

  inner_prod_numerator = :( poly[1][s1,i1]*f[$(s_vars...)] )
  for i = 2:N
      inner_prod_numerator = :( poly[$i][$(s_vars[i]),$(i_vars[i])]*$inner_prod_numerator )
  end
  numerator_term = inner_prod_numerator

  # Construct weights term

  weights_term = :( weights[$(i_vars...)] = numerator/denominator )

  # Construct the inner loops that compute numerator and denominator

  inner = :( numerator += $numerator_term;
             denominator += $denominator_term )
  outer  = inner

  for i = N:-1:1#1:N
    outer = :(
      for $(s_vars[i]) = 1:size(f,$i)
        $outer
      end
    )
  end

  # Construct the outer loops that compute the weights

  new_inner = :( numerator = zero(T);
                 denominator = zero(T);
                 if sum(tuple($(i_vars...))) <= order+N;
                   $outer;
                   $weights_term;
                 end )
  new_outer = new_inner

  for i = N:-1:1#1:N
    new_outer = :(
      for $(i_vars[i]) = 1:(order+1)
        $new_outer
      end
    )
  end

  # Initialize the weight Array

  initial_weight = string("weights = zeros(order+1,",)
  for i = 2:N
    initial_weight = string(initial_weight,"order+1,")
  end
  initial_weight = Meta.parse(string(initial_weight,")"))

  # Put it all together to compute the weights array

  final = :( $chebyshev_polynomials;
#             weights = Array{T}(undef,(order.+1)...);
             $initial_weight;
             $new_outer;
             return weights )

  return final

end

@generated function chebyshev_weights(f::AbstractArray{T,N},poly::NTuple{N,Array{T,2}},order::S) where {T,N,S}

  i_vars = Array{Symbol}(undef,N)
  s_vars = Array{Symbol}(undef,N)
  for i = 1:N
    i_vars[i] = Symbol("i$i")
    s_vars[i] = Symbol("s$i")
  end

  # Construct denominator term

  inner_prod_denominator = :( poly[1][s1,i1] )
  for i = 2:N
    inner_prod_denominator = :( poly[$i][$(s_vars[i]),$(i_vars[i])]*$inner_prod_denominator )
  end
  denominator_term = :( $inner_prod_denominator*$inner_prod_denominator )

  # Construct numerator term

  inner_prod_numerator = :( poly[1][s1,i1]*f[$(s_vars...)] )
  for i = 2:N
      inner_prod_numerator = :( poly[$i][$(s_vars[i]),$(i_vars[i])]*$inner_prod_numerator )
  end
  numerator_term = inner_prod_numerator

  # Construct weights term

  weights_term = :( weights[$(i_vars...)] = numerator/denominator )

  # Construct the inner loops that compute numerator and denominator

  inner = :( numerator += $numerator_term;
             denominator += $denominator_term )
  outer  = inner

  for i = N:-1:1#1:N
    outer = :(
      for $(s_vars[i]) = 1:size(f,$i)
        $outer
      end
    )
  end

  # Construct the outer loops that compute the weights

  new_inner = :( numerator = zero(T);
                 denominator = zero(T);
                 if sum(tuple($(i_vars...))) <= order+N;
                   $outer;
                   $weights_term;
                 end )
  new_outer = new_inner

  for i = N:-1:1#1:N
    new_outer = :(
      for $(i_vars[i]) = 1:(order+1)
        $new_outer
      end
    )
  end

  # Initialize the weight Array

  initial_weight = string("weights = zeros(order+1,",)
  for i = 2:N
    initial_weight = string(initial_weight,"order+1,")
  end
  initial_weight = Meta.parse(string(initial_weight,")"))

  # Put it all together to compute the weights array

  final = :( #weights = Array{T}(undef,(order.+1)...);
             $initial_weight;
             $new_outer;
             return weights )

  return final

end

# Functions to ensure backward compatibility with an older API

# Tensor product polynomials

function chebyshev_weights(f::AbstractArray{T,1},nodes_1::Array{T,1},order::Array{S,1},domain=[ones(1);-ones(1)]) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights(f,(nodes_1,),order,reshape(domain,2,1))

  return weights

end

function chebyshev_weights(f::AbstractArray{T,2},nodes_1::Array{T,1},nodes_2::Array{T,1},order::Array{S,1},domain=[ones(1);-ones(1)]) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights(f,(nodes_1,nodes_2),order,domain)

  return weights

end

function chebyshev_weights(f::AbstractArray{T,3},nodes_1::Array{T,1},nodes_2::Array{T,1},nodes_3::Array{T,1},order::Array{S,1},domain=[ones(1);-ones(1)]) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights(f,(nodes_1,nodes_2,nodes_3),order,domain)

  return weights

end

function chebyshev_weights(f::AbstractArray{T,4},nodes_1::Array{T,1},nodes_2::Array{T,1},nodes_3::Array{T,1},nodes_4::Array{T,1},order::Array{S,1},domain=[ones(1);-ones(1)]) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights(f,(nodes_1,nodes_2,nodes_3,nodes_4),order,domain)

  return weights

end

function chebyshev_weights(f::AbstractArray{T,5},nodes_1::Array{T,1},nodes_2::Array{T,1},nodes_3::Array{T,1},nodes_4::Array{T,1},nodes_5::Array{T,1},order::Array{S,1},domain=[ones(1);-ones(1)]) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights(f,(nodes_1,nodes_2,nodes_3,nodes_4,nodes_5),order,domain)

  return weights

end

function chebyshev_weights(f::AbstractArray{T,6},nodes_1::Array{T,1},nodes_2::Array{T,1},nodes_3::Array{T,1},nodes_4::Array{T,1},nodes_5::Array{T,1},nodes_6::Array{T,1},order::Array{S,1},domain=[ones(1);-ones(1)]) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights(f,(nodes_1,nodes_2,nodes_3,nodes_4,nodes_5,nodes_6),order,domain)

  return weights

end

# Complete polynomials

function chebyshev_weights(f::AbstractArray{T,1},nodes_1::Array{T,1},order::S,domain=[ones(1);-ones(1)]) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights(f,(nodes_1,),order,reshape(domain,2,1))

  return weights

end

function chebyshev_weights(f::AbstractArray{T,2},nodes_1::Array{T,1},nodes_2::Array{T,1},order::S,domain=[ones(1);-ones(1)]) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights(f,(nodes_1,nodes_2),order,domain)

  return weights

end

function chebyshev_weights(f::AbstractArray{T,3},nodes_1::Array{T,1},nodes_2::Array{T,1},nodes_3::Array{T,1},order::S,domain=[ones(1);-ones(1)]) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights(f,(nodes_1,nodes_2,nodes_3),order,domain)

  return weights

end

function chebyshev_weights(f::AbstractArray{T,4},nodes_1::Array{T,1},nodes_2::Array{T,1},nodes_3::Array{T,1},nodes_4::Array{T,1},order::S,domain=[ones(1);-ones(1)]) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights(f,(nodes_1,nodes_2,nodes_3,nodes_4),order,domain)

  return weights

end

function chebyshev_weights(f::AbstractArray{T,5},nodes_1::Array{T,1},nodes_2::Array{T,1},nodes_3::Array{T,1},nodes_4::Array{T,1},nodes_5::Array{T,1},order::S,domain=[ones(1);-ones(1)]) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights(f,(nodes_1,nodes_2,nodes_3,nodes_4,nodes_5),order,domain)

  return weights

end

function chebyshev_weights(f::AbstractArray{T,6},nodes_1::Array{T,1},nodes_2::Array{T,1},nodes_3::Array{T,1},nodes_4::Array{T,1},nodes_5::Array{T,1},nodes_6::Array{T,1},order::S,domain=[ones(1);-ones(1)]) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights(f,(nodes_1,nodes_2,nodes_3,nodes_4,nodes_5,nodes_6),order,domain)

  return weights

end

max_issues_repo_name: RJDennis/ChebyshevApprox

max_forks_repo_name: RJDennis/ChebyshevApprox


text:
module OrbitsWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AttitudeProvider = @jimport org.orekit.attitudes.AttitudeProvider
const CR3BPDifferentialCorrection = @jimport org.orekit.orbits.CR3BPDifferentialCorrection
const CR3BPSystem = @jimport org.orekit.bodies.CR3BPSystem
const CartesianOrbit = @jimport org.orekit.orbits.CartesianOrbit
const CircularOrbit = @jimport org.orekit.orbits.CircularOrbit
const Class = @jimport java.lang.Class
const Collection = @jimport java.util.Collection
const Enum = @jimport java.lang.Enum
const EquinoctialOrbit = @jimport org.orekit.orbits.EquinoctialOrbit
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const FieldCartesianOrbit = @jimport org.orekit.orbits.FieldCartesianOrbit
const FieldCircularOrbit = @jimport org.orekit.orbits.FieldCircularOrbit
const FieldEquinoctialOrbit = @jimport org.orekit.orbits.FieldEquinoctialOrbit
const FieldKeplerianOrbit = @jimport org.orekit.orbits.FieldKeplerianOrbit
const FieldOrbit = @jimport org.orekit.orbits.FieldOrbit
const FieldPVCoordinates = @jimport org.orekit.utils.FieldPVCoordinates
const FieldTimeInterpolable = @jimport org.orekit.time.FieldTimeInterpolable
const FieldTimeShiftable = @jimport org.orekit.time.FieldTimeShiftable
const Frame = @jimport org.orekit.frames.Frame
const HaloOrbit = @jimport org.orekit.orbits.HaloOrbit
const JString = @jimport java.lang.String
const KeplerianOrbit = @jimport org.orekit.orbits.KeplerianOrbit
const LagrangianPoints = @jimport org.orekit.utils.LagrangianPoints
const LibrationOrbit = @jimport org.orekit.orbits.LibrationOrbit
const LibrationOrbitFamily = @jimport org.orekit.orbits.LibrationOrbitFamily
const LibrationOrbitType = @jimport org.orekit.orbits.LibrationOrbitType
const LyapunovOrbit = @jimport org.orekit.orbits.LyapunovOrbit
const Object = @jimport java.lang.Object
const Optional = @jimport java.util.Optional
const Orbit = @jimport org.orekit.orbits.Orbit
const OrbitType = @jimport org.orekit.orbits.OrbitType
const PVCoordinates = @jimport org.orekit.utils.PVCoordinates
const ParameterDriversList = @jimport org.orekit.utils.ParameterDriversList
const PositionAngle = @jimport org.orekit.orbits.PositionAngle
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const RichardsonExpansion = @jimport org.orekit.orbits.RichardsonExpansion
const STMEquations = @jimport org.orekit.propagation.numerical.cr3bp.STMEquations
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const Stream = @jimport java.util.stream.Stream
const TimeInterpolable = @jimport org.orekit.time.TimeInterpolable
const TimeScale = @jimport org.orekit.time.TimeScale
const TimeStampedFieldPVCoordinates = @jimport org.orekit.utils.TimeStampedFieldPVCoordinates
const TimeStampedPVCoordinates = @jimport org.orekit.utils.TimeStampedPVCoordinates

export CR3BPDifferentialCorrection
export CartesianOrbit
export CircularOrbit
export EquinoctialOrbit
export FieldCartesianOrbit
export FieldCircularOrbit
export FieldEquinoctialOrbit
export FieldKeplerianOrbit
export FieldOrbit
export HaloOrbit
export KeplerianOrbit
export LibrationOrbit
export LibrationOrbitFamily
export LibrationOrbitType
export LyapunovOrbit
export Orbit
export OrbitType
export PositionAngle
export RichardsonExpansion
export add_kepler_contribution
export apply_differential_correction
export compare_to
export compute
export compute_halo_first_guess
export compute_lyapunov
export compute_lyapunov_first_guess
export convert_type
export describe_constable
export eccentric_to_mean
export eccentric_to_true
export elliptic_eccentric_to_mean
export elliptic_eccentric_to_true
export equals
export get_a
export get_a_dot
export get_alpha
export get_alpha_dot
export get_alpha_e
export get_alpha_e_dot
export get_alpha_m
export get_alpha_m_dot
export get_alpha_v
export get_alpha_v_dot
export get_anomaly
export get_anomaly_dot
export get_circular_ex
export get_circular_ex_dot
export get_circular_ey
export get_circular_ey_dot
export get_class
export get_cr3bp_system
export get_date
export get_declaring_class
export get_drivers
export get_e
export get_e_dot
export get_eccentric_anomaly
export get_eccentric_anomaly_dot
export get_equinoctial_ex
export get_equinoctial_ex_dot
export get_equinoctial_ey
export get_equinoctial_ey_dot
export get_frame
export get_halo_orbital_period
export get_hx
export get_hx_dot
export get_hy
export get_hy_dot
export get_i
export get_i_dot
export get_initial_pv
export get_jacobian_wrt_cartesian
export get_jacobian_wrt_parameters
export get_keplerian_mean_motion
export get_keplerian_period
export get_l
export get_l_dot
export get_lagrangian_point
export get_le
export get_le_dot
export get_lm
export get_lm_dot
export get_lv
export get_lv_dot
export get_lyapunov_orbital_period
export get_manifolds
export get_mean_anomaly
export get_mean_anomaly_dot
export get_mu
export get_orbital_period
export get_perigee_argument
export get_perigee_argument_dot
export get_pv_coordinates
export get_right_ascension_of_ascending_node
export get_right_ascension_of_ascending_node_dot
export get_true_anomaly
export get_true_anomaly_dot
export get_type
export has_derivatives
export hash_code
export hyperbolic_eccentric_to_mean
export hyperbolic_eccentric_to_true
export interpolate
export map_array_to_orbit
export map_orbit_to_array
export mean_to_eccentric
export mean_to_elliptic_eccentric
export mean_to_hyperbolic_eccentric
export name
export notify
export notify_all
export ordinal
export shifted_by
export to_orbit
export to_string
export true_to_eccentric
export true_to_elliptic_eccentric
export true_to_hyperbolic_eccentric
export value_of
export values
export wait

include("c_r3_bp_differential_correction.jl")
include("cartesian_orbit.jl")
include("circular_orbit.jl")
include("equinoctial_orbit.jl")
include("field_cartesian_orbit.jl")
include("field_circular_orbit.jl")
include("field_equinoctial_orbit.jl")
include("field_keplerian_orbit.jl")
include("field_orbit.jl")
include("halo_orbit.jl")
include("keplerian_orbit.jl")
include("libration_orbit.jl")
include("libration_orbit_family.jl")
include("libration_orbit_type.jl")
include("lyapunov_orbit.jl")
include("orbit.jl")
include("orbit_type.jl")
include("position_angle.jl")
include("richardson_expansion.jl")

end