text:
import Base.BaseDocs: @kw_str

"""
**欢迎来到 Julia $(string(VERSION)).** 完整的中文手册可以在这里找到

    https://docs.juliacn.com/

更多中文资料和教程，也请关注 Julia 中文社区

    https://cn.julialang.org

新手请参考中文 discourse 上的新手指引

    https://discourse.juliacn.com/t/topic/159

输入 `?`， 然后输入你想要查看帮助文档的函数或者宏名称就可以查看它们的文档。例如 `?cos`, 或者 `?@time` 然后按回车键即可。

在 REPL 中输入 `ENV["REPL_LOCALE"]=""` 将恢复英文模式。再次回到中文模式请输入 `ENV["REPL_LOCALE"]="zh_CN"`。
"""
kw"help", kw"?", kw"julia", kw""

"""
    using

`using Foo` 将会加载一个名为 `Foo` 的模块（module）或者一个包，然后其 [`export`](@ref) 的名称将可以直接使用。不论是否被 `export`，名称都可以通过点来访问（例如，输入 `Foo.foo` 来访问到 `foo`）。查看[手册中关于模块的部分](@ref modules)以获取更多细节。
"""
kw"using"

"""
    import

`import Foo` 将会加载一个名为 `Foo` 的模块（module）或者一个包。`Foo` 模块中的名称可以通过点来访问到（例如，输入 `Foo.foo` 可以获取到 `foo`）。查看[手册中关于模块的部分](@ref modules)以获取更多细节。
"""
kw"import"

"""
    export

`export` 被用来在模块中告诉Julia哪些函数或者名字可以由用户使用。例如 `export foo` 将在 [`using`](@ref) 这个 module 的时候使得 `foo`可以直接被访问到。查看[手册中关于模块的部分](@ref modules)以获取更多细节。
"""
kw"export"

"""
    abstract type

`abstract type` 声明来一个不能实例化的类型，它将仅仅作为类型图中的一个节点存在，从而能够描述一系列相互关联的具体类型（concrete type）：这些具体类型都是抽象类型的子节点。抽象类型在概念上使得 Julia 的类型系统不仅仅是一系列对象的集合。例如：

```julia
abstract type Number end
abstract type Real <: Number end
```

[`Number`](@ref) 没有父节点（父类型）, 而 [`Real`](@ref) 是 `Number` 的一个抽象子类型。
"""
kw"abstract type"

"""
    module

`module` 会声明一个 `Module` 类型的实例用于描述一个独立的变量名空间。在一个模块（module）里，你可以控制来自于其它模块的名字是否可见（通过载入，import），你也可以决定你的名字有哪些是可以公开的（通过暴露，export）。模块使得你在在创建上层定义时无需担心命名冲突。查看[手册中关于模块的部分](@ref modules)以获取更多细节。

# 例子
```julia
module Foo
import Base.show
export MyType, foo

struct MyType
    x
end

bar(x) = 2x
foo(a::MyType) = bar(a.x) + 1
show(io::IO, a::MyType) = print(io, "MyType \$(a.x)")
end
```
"""
kw"module"

"""
    baremodule

`baremodule` 将声明一个不包含 `using Base` 或者 `eval` 定义的模块。但是它将仍然载入 `Core` 模块。
"""
kw"baremodule"

"""
    primitive type

`primitive type` 声明了一个其数据仅仅由一系列二进制数表示的具体类型。比较常见的例子是整数类型和浮点类型。下面是一些内置的原始类型（primitive type）：

```julia
primitive type Char 32 end
primitive type Bool <: Integer 8 end
```

名称后面的数字表达了这个类型存储所需的比特数目。目前这个数字要求是 8 bit 的倍数。[`Bool`](@ref) 类型的声明展示了一个原始类型如何选择成为另一个类型的子类型。
"""
kw"primitive type"

"""
    macro

`macro` 定义了一种会将生成的代码包含在最终程序体中的方法，这称之为宏。一个宏将一系列输入映射到一个表达式，然后所返回的表达式将会被直接进行编译而不需要在运行时调用 `eval` 函数。宏的输入可以包括表达式、字面量和符号。例如：

# 例子

```jldoctest
julia> macro sayhello(name)
           return :( println("Hello, ", \$name, "!") )
       end
@sayhello (macro with 1 method)

julia> @sayhello "小明"
Hello, 小明!
```
"""
kw"macro"

"""
    local

`local`将会定义一个新的局部变量。

查看[手册：变量作用域](@ref scope-of-variables)以获取更详细的信息。

# 例子
```jldoctest
julia> function foo(n)
           x = 0
           for i = 1:n
               local x # introduce a loop-local x
               x = i
           end
           x
       end
foo (generic function with 1 method)

julia> foo(10)
0
```
"""
kw"local"

"""
    global

`global x` 将会使得当前作用域和当前作用所包含的作用域里的 `x` 指向名为 `x` 的全局变量。查看[手册：变量作用域](@ref scope-of-variables)以获取更多信息。

# 例子
```jldoctest
julia> z = 3
3

julia> function foo()
           global z = 6 # use the z variable defined outside foo
       end
foo (generic function with 1 method)

julia> foo()
6

julia> z
6
```
"""
kw"global"

"""
    let

`let` 会在每次被运行时声明一个新的变量绑定。这个新的变量绑定将拥有一个新的地址。这里的不同只有当变量通过闭包生存在它们的作用域外时才会显现。`let` 语法接受逗号分割的一系列赋值语句和变量名：

```julia
let var1 = value1, var2, var3 = value3
    code
end
```

这些赋值语句是按照顺序求值的，等号右边的表达式将会首先求值，然后才绑定给左边的变量。因此这使得 `let x = x` 这样的表达式有意义，因为这两个 `x` 变量将具有不同的地址。
"""
kw"let"

"""
    quote

`quote` 会将其包含的代码扩变成一个多重的表达式对象，而无需显示调用 `Expr` 的构造器。这称之为引用，比如说

```julia
ex = quote
    x = 1
    y = 2
    x + y
end
```

和其它引用方式不同的是，`:( ... )`形式的引用（被包含时）将会在表达式树里引入一个在操作表达式树时必须要考虑的 `QuoteNode` 元素。而在其它场景下，`:( ... )`和 `quote .. end` 代码块是被同等对待的。
"""
kw"quote"

"""
    '

厄米算符（共轭转置），参见 [`adjoint`](@ref)

# 例子
```jldoctest
julia> A = [1.0 -2.0im; 4.0im 2.0]
2×2 Array{Complex{Float64},2}:
 1.0+0.0im  -0.0-2.0im
 0.0+4.0im   2.0+0.0im

julia> A'
2×2 Array{Complex{Float64},2}:
  1.0-0.0im  0.0-4.0im
 -0.0+2.0im  2.0-0.0im
```
"""
kw"'"

"""
    const

`const` 被用来声明常数全局变量。在大部分（尤其是性能敏感的代码）全局变量应当被声明为常数。

```julia
const x = 5
```

可以使用单个 `const` 声明多个常数变量。

```julia
const y, z = 7, 11
```

注意 `const` 只会作用于一个 `=` 操作，因此 `const x = y = 1` 声明了 `x` 是常数，而 `y` 不是。在另一方面，`const x = const y = 1`声明了 `x` 和 `y` 都是常数。

注意「常数性质」并不会强制容器内部变成常数，所以如果 `x` 是一个数组或者字典（举例来讲）你仍然可以给它们添加或者删除元素。

严格来讲，你甚至可以重新定义 `const`（常数）变量，尽管这将会让编译器产生一个警告。唯一严格的要求是这个变量的**类型**不能改变，这也是为什么常数变量会比一般的全局变量更快的原因。
"""
kw"const"

"""
    function

函数由 `function` 关键词定义：

```julia
function add(a, b)
    return a + b
end
```

或者是更短的形式：

```julia
add(a, b) = a + b
```

[`return`](@ref) 关键词的使用方法和其它语言完全一样，但是常常是不使用的。一个没有显示声明 `return` 的函数将返回函数体最后一个表达式。
"""
kw"function"

"""
    return

`return` 可以用来在函数体中立即退出并返回给定值，例如

```julia
function compare(a, b)
    a == b && return "equal to"
    a < b ? "less than" : "greater than"
end
```

通常，你可以在函数体的任意位置放置 `return` 语句，包括在多层嵌套的循环和条件表达式中，但要注意 `do` 块。例如：

```julia
function test1(xs)
    for x in xs
        iseven(x) && return 2x
    end
end

function test2(xs)
    map(xs) do x
        iseven(x) && return 2x
        x
    end
end
```

在第一个例子中，return 一碰到偶数就跳出包含它的函数，因此 `test1([5,6,7])` 返回 `12`。

你可能希望第二个例子的行为与此相同，但实际上，这里的 `return` 只会跳出（在 `do` 块中的）*内部*函数并把值返回给 `map`。于是，`test2([5,6,7])` 返回 `[5,12,7]`。
"""
kw"return"

# 仿照 https://docs.juliacn.com/latest/manual/control-flow/#man-conditional-evaluation-1
"""
    if/elseif/else

`if`/`elseif`/`else` 执行条件表达式（Conditional evaluation）可以根据布尔表达式的值，让部分代码被执行或者不被执行。下面是对 `if`-`elseif`-`else` 条件语法的分析：

```julia
if x < y
    println("x is less than y")
elseif x > y
    println("x is greater than y")
else
    println("x is equal to y")
end
```

如果表达式 `x < y` 是 `true`，那么对应的代码块会被执行；否则判断条件表达式 `x > y`，如果它是 `true`，则执行对应的代码块；如果没有表达式是 true，则执行 `else` 代码块。`elseif` 和 `else` 代码块是可选的，并且可以使用任意多个 `elseif` 代码块。
"""
kw"if", kw"elseif", kw"else"

"""
    for

`for` 循环通过迭代一系列值来重复计算循环体。

# 例子
```jldoctest
julia> for i in [1, 4, 0]
           println(i)
       end
1
4
0
```
"""
kw"for"

# 仿照 https://docs.juliacn.com/latest/manual/control-flow/#man-loops-1
"""
    while

`while` 循环会重复执行条件表达式，并在该表达式为 `true` 时继续执行 while 循环的主体部分。当 while 循环第一次执行时，如果条件表达式为 false，那么主体代码就一次也不会被执行。

# 例子
jldoctest
julia> i = 1
1

julia> while i < 5
           println(i)
           global i += 1
       end
1
2
3
4
```
"""
kw"while"

"""
    end

`end` 标记一个表达式块的结束，例如 [`module`](@ref)、[`struct`](@ref)、[`mutable struct`](@ref)、[`begin`](@ref)、[`let`](@ref)、[`for`](@ref) 等。`end` 在索引数组时也可以用来表示维度的最后一个索引。

# 例子
```jldoctest
julia> A = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> A[end, :]
2-element Array{Int64,1}:
 3
 4
```
"""
kw"end"

# 仿照 https://docs.juliacn.com/latest/manual/control-flow/#try/catch-%E8%AF%AD%E5%8F%A5-1
"""
    try/catch

`try/catch` 语句可以用来捕获 `Exception`，并进行异常处理。例如，一个自定义的平方根函数可以通过 `Exception` 来实现自动按需调用求解实数或者复数平方根的方法：

```julia
f(x) = try
    sqrt(x)
catch
    sqrt(complex(x, 0))
end
```

`try/catch` 语句允许保存 `Exception` 到一个变量中，例如 `catch y`。

`try/catch` 组件的强大之处在于能够将高度嵌套的计算立刻解耦成更高层次地调用函数。
"""
kw"try", kw"catch"

# 仿照 https://docs.juliacn.com/latest/manual/control-flow/#finally-%E5%AD%90%E5%8F%A5-1
"""
    finally

无论代码块是如何退出的，都让代码块在退出时运行某段代码。这里是一个确保一个打开的文件被关闭的例子：

```julia
f = open("file")
try
    operate_on_file(f)
finally
    close(f)
end
```

当控制流离开 `try` 代码块（例如，遇到 `return`，或者正常结束），`close(f)` 就会被执行。如果 `try` 代码块由于异常退出，这个异常会继续传递。`catch` 代码块可以和 `try` 还有 `finally` 配合使用。这时 `finally` 代码块会在 `catch` 处理错误之后才运行。
"""
kw"finally"

"""
    break

立即跳出当前循环。

# 例子
```jldoctest
julia> i = 0
0

julia> while true
           global i += 1
           i > 5 && break
           println(i)
       end
1
2
3
4
5
```
"""
kw"break"

"""
    continue

跳过当前循环迭代的剩余部分。

# 例子
```jldoctest
julia> for i = 1:6
           iseven(i) && continue
           println(i)
       end
1
3
5
```
"""
kw"continue"

"""
    do

创建一个匿名函数。例如：

```julia
map(1:10) do x
    2x
end
```

等价于 `map(x->2x, 1:10)`。

像这样便可使用多个参数：

```julia
map(1:10, 11:20) do x, y
    x + y
end
```
"""
kw"do"

"""
    ...

「splat」运算符 `...` 表示参数序列。`...` 可以在函数定义中用来表示该函数接受任意数量的参数。`...` 也可以用来将函数作用于参数序列。

# 例子
```jldoctest
julia> add(xs...) = reduce(+, xs)
add (generic function with 1 method)

julia> add(1, 2, 3, 4, 5)
15

julia> add([1, 2, 3]...)
6

julia> add(7, 1:100..., 1000:1100...)
111107
```
"""
kw"..."

"""
    ;

`;` 在 Julia 中具有与许多类 C 语言相似的作用，用于分隔前一个语句的结尾。`;` 在换行中不是必要的，但可以用于在单行中分隔语句或者将多个表达式连接为单个表达式。`;` 也用于抑制 REPL 和类似界面中的输出打印。

# 例子
```julia
julia> function foo()
           x = "Hello, "; x *= "World!"
           return x
       end
foo (generic function with 1 method)

julia> bar() = (x = "Hello, Mars!"; return x)
bar (generic function with 1 method)

julia> foo();

julia> bar()
"Hello, Mars!"
```
"""
kw";"

"""
    x && y

短路布尔 AND。
"""
kw"&&"

"""
    x || y

短路布尔 OR。
"""
kw"||"

"""
    ccall((function_name, library), returntype, (argtype1, ...), argvalue1, ...)
    ccall(function_name, returntype, (argtype1, ...), argvalue1, ...)
    ccall(function_pointer, returntype, (argtype1, ...), argvalue1, ...)

调用由 C 导出的共享库里的函数，该函数由元组 `(function_name, library)` 指定，其中的每个组件都是字符串或符号。若不指定库，还可以使用 `function_name` 的符号或字符串，它会在当前进程中解析。另外，`ccall` 也可用于调用函数指针 `function_pointer`，比如 `dlsym` 返回的函数指针。

请注意参数类型元组必须是字面上的元组，而不是元组类型的变量或表达式。

通过自动插入对 `unsafe_convert(argtype, cconvert(argtype, argvalue))` 的调用，每个传给 `ccall` 的 `argvalue` 将被类型转换为对应的 `argtype`。（有关的详细信息，请参阅 [`unsafe_convert`](@ref Base.unsafe_convert) 和 [`cconvert`](@ref Base.cconvert) 的文档。）在大多数情况下，这只会简单地调用 `convert(argtype, argvalue)`。
"""
kw"ccall"

"""
    begin

`begin...end` 表示一个代码块。

```julia
begin
    println("Hello, ")
    println("World!")
end
```

通常，`begin` 不会是必需的，因为诸如 [`function`](@ref) and [`let`](@ref) 之类的关键字会隐式地开始代码块。另请参阅 [`;`](@ref)。
"""
kw"begin"

"""
    struct

struct 是 Julia 中最常用的数据类型，由名称和一组字段指定。

```julia
struct Point
    x
    y
end
```

可对字段施加类型限制，该限制也可被参数化：

```julia
    struct Point{X}
        x::X
        y::Float64
    end
```

struct 可以通过 `<:` 语法声明一个抽象超类型：
A struct can also declare an abstract super type via `<:` syntax:

```julia
struct Point <: AbstractPoint
    x
    y
end
```

`struct` 默认是不可变的；这些类型的实例在构造后不能被修改。如需修改实例，请使用 [`mutable struct`](@ref) 来声明一个可以修改其实例的类型。

有关更多细节，比如怎么定义构造函数，请参阅手册的 [复合类型](@ref) 章节。
"""
kw"struct"

"""
    mutable struct

`mutable struct` 类似于 [`struct`](@ref)，但另外允许在构造后设置类型的字段。有关详细信息，请参阅 [复合类型](@ref)。
"""
kw"mutable struct"

"""
    new

仅在内部构造函数中可用的特殊函数，用来创建该类型的对象。有关更多信息，请参阅手册的 [内部构造方法](@ref) 章节。
"""
kw"new"

"""
    where

`where` 关键字创建一个类型，该类型是其他类型在一些变量上所有值的迭代并集。例如 `Vector{T} where T<:Real` 包含所有元素类型是某种 `Real` 的 [`Vector`](@ref)。

如果省略，变量上界默认为 `Any`：

```julia
Vector{T} where T    # short for `where T<:Any`
```

变量也可以具有下界：

```julia
Vector{T} where T>:Int
Vector{T} where Int<:T<:Real
```

嵌套的 `where` 也有简洁的语法。例如，这行代码：

```julia
Pair{T, S} where S<:Array{T} where T<:Number
```

可以缩写为：

```julia
Pair{T, S} where {T<:Number, S<:Array{T}}
```

这种形式常见于方法签名：

请注意，在这种形式中，最外层变量列在最前面。这与使用语法 `T{p1, p2, ...}` 将类型「作用」于参数值时所替换变量的次序相匹配。
"""
kw"where"

"""
    ans

一个引用最后一次计算结果的变量，在交互式提示符中会自动设置。
"""
kw"ans"

"""
    Union{}

`Union{}`，即空的类型 [`Union`](@ref)，是没有值的类型。也就是说，它具有决定性性质：对于任何 `x`，`isa(x, Union{}) == false`。`Base.Bottom` 被定义为其别名，`Union{}` 的类型是 `Core.TypeofBottom`。

# 例子
```jldoctest
julia> isa(nothing, Union{})
false
```
"""
kw"Union{}", Base.Bottom

"""
    ::

`::` 运算符用于类型声明，在程序中可被附加到表达式和变量后。详见手册的 [类型声明](@ref) 章节

在类型声明外，`::` 用于断言程序中的表达式和变量具有给定类型。

# 例子
```jldoctest
julia> (1+2)::AbstractFloat
ERROR: TypeError: typeassert: expected AbstractFloat, got Int64

julia> (1+2)::Int
3
```
"""
kw"::"

"""
Julia 的基础库。
"""
kw"Base"


text:
using Test
using StaticArrays
using BenchmarkTools

using iLQGames:
    Unicycle4D,
    NPlayerUnicycleCost,
    SystemTrajectory,
    AffineStrategy,
    LinearizationStyle,
    FeedbackLinearization,
    TrivialLinearization,
    iLQSolver,
    GeneralGame,
    NPlayer2DDoubleIntegratorCost,
    generate_nplayer_navigation_game,
    dynamics,
    n_states,
    n_controls,
    n_players,
    player_costs,
    horizon,
    solve!,
    feedbacklin,
    player_id,
    xindex,
    uindex,
    goalcost,
    ξ_from,
    subsystems,
    transform_to_feedbacklin,
    feedbacklin,
    x_from,
    λ_issingular,
    sparse_feedbacklin,
    next_x

using iLQGames.TestUtils

# generate a game
T_horizon = 10.
ΔT = 0.1

x01 = @SVector [-3., 0., 0., 0.]
x02 = @SVector [0.,  3., 0., 0.]
x0 = vcat(x01, x02)
# goal states (goal position of other player with opposite orientation)
xg1 = @SVector [3., 0., 0., 0.]
xg2 = @SVector [0., -3., 0., 0.]
g = generate_nplayer_navigation_game(Unicycle4D, NPlayerUnicycleCost, T_horizon, ΔT,
                                     xg1, xg2)
@test LinearizationStyle(dynamics(g)) isa FeedbackLinearization

dyn = dynamics(g)
nx = n_states(dyn)
nu = n_controls(dyn)
pcs = player_costs(g)
h = horizon(g)
zero_op = zero(SystemTrajectory{h, ΔT, nx, nu})
# solve the lq game
solver = iLQSolver(g)
# - setup initial_strategy
steer_init(k::Int) = cos(k/h*pi) * deg2rad(0)
acc_init(k::Int) = -cos(k/h*pi)*0.3
γ_init = SizedVector{h}([AffineStrategy((@SMatrix zeros(nu, nx)),
                                 (@SVector [steer_init(k), 0.7*acc_init(k),
                                            steer_init(k), acc_init(k)])) for k in 1:h])
@info "Benchmark *without* feedback linearization:"
@benchmark(solve!(copy(zero_op), copy(γ_init), $g, $solver, $x0)) |> display

steer_init(k::Int) = cos(k/h*pi) * deg2rad(0)
acc_init(k::Int) = -cos(k/h*pi)*0.3
γξ_init = SizedVector{h}([AffineStrategy((@SMatrix zeros(nu, nx)),
                                 (@SVector [steer_init(k), 0.7*acc_init(k),
                                            steer_init(k), acc_init(k)])) for k in 1:h])

gξ, ξ0 = transform_to_feedbacklin(g, x0)
@test LinearizationStyle(dynamics(gξ)) isa TrivialLinearization
solverξ = iLQSolver(gξ)
γξ_init = SizedVector{h}([AffineStrategy((@SMatrix zeros(nu, nx)),
                                 (@SVector [-0.05, 0.02,
                                            0.02, 0.05])) for k in 1:h])

quad_sanity_check(gξ)

@info "Benchmark *with* feedback linearization:"
@benchmark(solve!(copy(zero_op), copy(γξ_init), $gξ, $solverξ, $ξ0)) |> display

@testset "state conversion" begin
    for i in 1:100
        ξ = randn(SVector{n_states(dyn)})
        x = x_from(dyn, ξ)
        if !λ_issingular(dyn, ξ)
            @test ξ ≈ ξ_from(dyn, x)
        end
    end
end;

@testset "sprase feedback linearization" begin
    dyn_fblin_dense = feedbacklin(dyn)
    dyn_fblin_sparse = sparse_feedbacklin(dyn)
    # these should have the same behavior
    for i in 1:100
        x = randn(SVector{n_states(dyn_fblin_dense)})
        u = randn(SVector{n_controls(dyn_fblin_dense)})
        t = 0.0

        @test next_x(dyn_fblin_dense, x, u, t) ≈ next_x(dyn_fblin_sparse, x, u, t)
    end
end;

ext: jl


text:
@testset "Vulkan Specification - DataFrames" begin
    include("spec/dataframes.jl")
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((:poly1, :dom1, true, -4),),
    ((:poly1, :dom1, false, -4), :SOCExpPSD),
    ]
insts["fast"] = [
    ((:poly1, :dom2, true, -4),),
    ((:poly1, :dom1, false, -4), :SOCExpPSD),
    ((:poly1, :dom2, false, -4), :SOCExpPSD),
    ((:poly2, :dom3, true, -2),),
    ((:poly2, :dom4, true, -2),),
    ((:poly2, :dom3, false, -2), :SOCExpPSD),
    ((:poly2, :dom4, false, -2), :SOCExpPSD),
    ]
insts["various"] = insts["fast"]
return (ConvexityParameterJuMP, insts)


text:
using PyPlot
using Peacock

function epf(x,y)
    # equation of a circle with radius 0.2a
    if x^2+y^2 <= 0.2^2
        # dielectric inside the circle
        return 8.9
    else
        # air outside the circle
        return 1
    end
end

# Permeability is unity everywhere
function muf(x,y)
    return 1
end

a1 = [1, 0]  # first lattice vector
a2 = [0, 1]  # second lattice vector
d1 = 0.01  # resolution along first lattice vector
d2 = 0.01  # resolution along second lattice vector
geometry = Geometry(epf, muf, a1, a2, d1, d2)
fourier_space_cutoff = 9
solver_CPU = Solver(geometry, fourier_space_cutoff)
solver_GPU = Solver(geometry, fourier_space_cutoff, GPU=true)
G = BrillouinZoneCoordinate(  0,   0, "Γ")
X = BrillouinZoneCoordinate(1/2,   0, "X")
M = BrillouinZoneCoordinate(1/2, 1/2, "M")
ks = [G,X,M,G]

@show fourier_space_cutoff
@show typeof(solver_CPU.epc)
@show typeof(solver_GPU.muc)

function test_bands(solver, dk)
    figure(figsize=(4,3))
    plot_band_diagram(solver, ks, TE, color="red",
                bands=1:4, dk=dk, frequency_scale=1/2pi)
    plot_band_diagram(solver, ks, TM, color="blue",
                bands=1:4, dk=dk, frequency_scale=1/2pi)
    ylim(0,0.8)
end

# call once to make sure functions are compiled
test_bands(solver_CPU, 2)
test_bands(solver_GPU, 2)

close("all")
@time test_bands(solver_CPU, 0.1)
@time test_bands(solver_GPU, 0.1)
show()


text:
# ---
# title: 1317. Convert Integer to the Sum of Two No-Zero Integers
# id: problem1317
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/convert-integer-to-the-sum-of-two-no-zero-integers/description/>
# hidden: true
# ---
#
# Given an integer `n`. No-Zero integer is a positive integer which **doesn 't
# contain any 0** in its decimal representation.
#
# Return _a list of two integers_ `[A, B]` where:
#
#   * `A` and `B` are No-Zero integers.
#   * `A + B = n`
#
# It's guarateed that there is at least one valid solution. If there are many
# valid solutions you can return any of them.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 2
#     Output: [1,1]
#     Explanation: A = 1, B = 1. A + B = n and both A and B don't contain any 0 in their decimal representation.
#
#
# **Example 2:**
#
#
#
#     Input: n = 11
#     Output: [2,9]
#
#
# **Example 3:**
#
#
#
#     Input: n = 10000
#     Output: [1,9999]
#
#
# **Example 4:**
#
#
#
#     Input: n = 69
#     Output: [1,68]
#
#
# **Example 5:**
#
#
#
#     Input: n = 1010
#     Output: [11,999]
#
#
#
#
# **Constraints:**
#
#   * `2 <= n <= 10^4`
#
#
## @lc code=start
using LeetCode

function get_no_zero_integers(n::Int)::Vector{Int}
    for i in 1:(n - 1)
        if !('0' in string(i)) && !('0' in string(n - i))
            return [i, n - i]
        end
    end
end
## @lc code=end


text:
function empiricaldescriptors(x, q, t, a, b, c, pbc, rcutmax, eta, kappa, emdescriptors)

    # compute full neighbor list
    y, alist, neighlist, neighnum = fullneighborlist(x, a, b, c, pbc, rcutmax);
    if !isempty(q)
        q = q[:,alist];
    end

    dim, N = size(x);
    d = []
    dd = []
    dv = []
    for i = 1:length(emdescriptors)

        style, bondtype, bondedatomtypes, potentialfunc, mu, rcut = getpotential(emdescriptors[i]);
        potentialfunc = getfield(Main, Symbol(potentialfunc));

        # single potentials
        if (style == "single") & (bondtype == "nonbonded")
            ilist = Array(1:N);
         elseif (style == "single") & (bondtype == "bonded")
            typei = bondedatomtypes[1];
            ilist = findatomtype(Array(1:N), t, typei);
        end
        if (style == "single")
            xi, qi, ai, ti = neighsingles(x, q, t, ilist);
            ei, fi = potentialfunc(xi, qi, ti, mu, eta, kappa);
            ea, fa = tallysingle(ei, fi, ai);
            d = [d; sum(ea)];
            if length(dd)==0
                dd = fa;
            else
                dd = cat(dd, fa, dims=3)
            end
            va = tallysinglevirial(fi, xi, ai, -1.0)
            if length(dv)==0
                dv = sum(va, dims=2)
            else
                dv = [dv sum(va, dims=2)]
            end
        end

        # pair potentials
        if (style == "pair") & (bondtype == "nonbonded")
            ilist = Int64.(Array(1:N));
            pairlist, pairnum = neighpairlist(y, ilist, neighlist, neighnum, rcut*rcut);
         elseif (style == "pair") & (bondtype == "bonded")
            typei = bondedatomtypes[1];
            typej = bondedatomtypes[2];

            ilist = findatomtype(Array(1:N), t, typei);
            pairlist, pairnum = neighpairlisttypej(y, ilist, alist, t, typej, neighlist, neighnum, rcut*rcut);
        end
        if (style == "pair")
            xij, ai, aj, ti, tj, qi, qj = neighpairs(y, q, pairlist, pairnum, t, ilist, alist);
            eij, fij = potentialfunc(xij, qi, qj, ti, tj, mu, eta, kappa);
            # tm0 = 1.0*fij
            # dim, na = size(xij)
            # for i1 = 1:dim
            #     for j1 = 1:na
            #         tm1 = 1.0*xij
            #         tm1[i1,j1] = tm1[i1,j1] - 1e-6
            #         tm2, ~ = potentialfunc(tm1, qi, qj, ti, tj, mu, eta, kappa);
            #         tm0[i1,j1,:] = (eij[j1,:] - tm2[j1,:])/1e-6
            #     end
            # end
            # display(abs.(fij-tm0))
            # error("here")
            for m = 1:size(eij,2)
                ea, fa = tallypair(eij[:,m], fij[:,:,m], ai, aj);
                d = [d; sum(ea)];
                if length(dd)==0
                    dd = fa;
                else
                    dd = cat(dd, fa, dims=3)
                end
                va = tallypairvirial(fij[:,:,m], xij, ai, aj, -1.0/2)
                if length(dv)==0
                    dv = sum(va, dims=2)
                else
                    dv = [dv sum(va, dims=2)]
                end
            end
        end

        # triplet potentials
        if (style == "triplet") & (bondtype == "nonbonded")
            ilist = Array(1:N);
            pairlist, pairnum = neighpairlist(y, ilist, neighlist, neighnum, rcut*rcut);
         elseif (style == "triplet") & (bondtype == "bonded")
            typei = bondedatomtypes[1];
            typej = bondedatomtypes[2];
            typek = bondedatomtypes[3];

            ilist = findatomtype(Array(1:N), t, typei);
            pairlist, pairnum = neighpairlisttypej(y, ilist, alist, t, [typej; typek], neighlist, neighnum, rcut*rcut);
        end
        if (style == "triplet")
            tripletlist, tripletnum = neightripletlist(pairlist, pairnum, ilist);
            xij, xik, ai, aj, ak, ti, tj, tk, qi, qj, qk = neightriplets(y, q, tripletlist, tripletnum, t, ilist, alist);
            eij, fij, fik = potentialfunc(xij, xik, qi, qj, qk, ti, tj, tk, mu, eta, kappa);
            # tmj = 1.0*fik
            # tmk = 1.0*fik
            # dim, na = size(xij)
            # na = 10;
            # for i1 = 1:dim
            #     for j1 = 1:na
            #         tm1 = 1.0*xij
            #         tm1[i1,j1] = tm1[i1,j1] - 1e-6
            #         tm2, ~ = potentialfunc(tm1, xik, qi, qj, qk, ti, tj, tk, mu, eta, kappa);
            #         tmj[i1,j1,:] = (eij[j1,:] - tm2[j1,:])/1e-6
            #         tm1 = 1.0*xik
            #         tm1[i1,j1] = tm1[i1,j1] - 1e-6
            #         tm2, ~ = potentialfunc(xij, tm1, qi, qj, qk, ti, tj, tk, mu, eta, kappa);
            #         tmk[i1,j1,:] = (eij[j1,:] - tm2[j1,:])/1e-6
            #     end
            # end
            # display(abs.(fij[:,1:na,:]-tmj[:,1:na,:]))
            # display(abs.(fik[:,1:na,:]-tmk[:,1:na,:]))
            # error("here")
            for m = 1:size(eij,2)
                ea, fa = tallytriplet(eij[:,m], fij[:,:,m], fik[:,:,m], ai, aj, ak);
                d = [d; sum(ea)];
                if length(dd)==0
                    dd = fa;
                else
                    dd = cat(dd, fa, dims=3)
                end
                va = tallytripletvirial(fij[:,:,m], fik[:,:,m], xij, xik, ai, aj, ak, -1.0/3)
                if length(dv)==0
                    dv = sum(va, dims=2)
                else
                    dv = [dv sum(va, dims=2)]
                end
            end
        end

        # quadlet potentials
        if (style == "quadlet") & (bondtype == "nonbonded")
            ilist = Array(1:N);
            pairlist, pairnum = neighpairlist(y, ilist, neighlist, neighnum, rcut*rcut);
         elseif (style == "quadlet") & (bondtype == "bonded")
            typei = bondedatomtypes[1];
            typej = bondedatomtypes[2];
            typek = bondedatomtypes[3];
            typel = bondedatomtypes[4];

            ilist = findatomtype(Array(1:N), t, typei);
            pairlist, pairnum = neighpairlisttypej(y, ilist, alist, t, [typej typek typel], neighlist, neighnum, rcut*rcut);
        end
        if (style == "quadlet")
            quadletlist, quadletnum = neighquadletlist(pairlist, pairnumsum, ilist);
            xij, xik, xil, ai, aj, ak, al, ti, tj, tk, tl, qi, qj, qk, ql =
                        neighquadlets(x, q, quadletlist, quadletnum, t, ilist, alist);
            eij, fij, fik, fil = potentialfunc(xij, xik, xil, qi, qj, qk, ql, ti, tj, tk, tl, mu, eta, kappa);
            for m = 1:size(eij,2)
                ea, fa = tallyquadlet(eij[:,m], fij[:,:,m], fik[:,:,m], fil[:,:,m], ai, aj, ak, al);
                d = [d; sum(ea)];
                if length(dd)==0
                    dd = fa;
                else
                    dd = cat(dd, fa, dims=3)
                end
                va = tallyquadletvirial(fij[:,:,m], fik[:,:,m], fil[:,:,m], xij, xik, xil, ai, aj, ak, al, -1.0/4)
                if length(dv)==0
                    dv = sum(va, dims=2)
                else
                    dv = [dv sum(va, dims=2)]
                end
            end
        end

    end

    return d, dd, dv

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
struct Difference{PA, PB}
    a::PA
    b::PB
end

Base.:*(n::Number, d::Difference) = Difference(*(n, d.a), *(n, d.b))
Base.:+(d1::Difference, d2::Difference) = Difference(d1.a + d2.a, d1.b + d2.b)
any_inside(d::Difference) = Difference(any_inside(d.a), any_inside(d.b))

struct CollisionCache{GeomA, GeomB, M, D <: Difference}
    bodyA::GeomA
    bodyB::GeomB
    simplex_points::MVector{M, D}
end

function CollisionCache(geomA, geomB)
    N = dimension(typeof(geomA))
    @assert dimension(typeof(geomA)) == dimension(typeof(geomB))
    CollisionCache(N, geomA, geomB)
end

function CollisionCache(::Val{N}, geomA, geomB) where N
    interior_point = any_inside(Difference(geomA, geomB))
    simplex = MVector(ntuple(i -> interior_point, Val(N + 1)))
    CollisionCache(geomA, geomB, simplex)
end

function reset!(cache::CollisionCache)
    interior_point = any_inside(Difference(cache.bodyA, cache.bodyB))
    simplex = cache.simplex_points
    @inbounds for i in eachindex(simplex)
        simplex[i] = interior_point
    end
    cache
end

dimension(::Type{CollisionCache{G1, G2, M, D}}) where {G1,G2,M,D} = dimension(G1)

function support_vector_max(geometry::gt.GeometryPrimitive, direction, initial_guess::Tagged)
    best_pt, score = gt.support_vector_max(geometry, gt.Vec(direction))
    Tagged(SVector(best_pt))
end

function support_vector_max(pt::gt.Vec{N, T}, direction, initial_guess::Tagged) where {N,T}
    Tagged(SVector(pt))
end

function support_vector_max(pt::gt.Point{N, T}, direction, initial_guess::Tagged) where {N,T}
    Tagged(SVector(pt))
end

function support_vector_max(simplex::Union{gt.AbstractSimplex, gt.AbstractFlexibleGeometry}, direction, initial_guess::Tagged)
    best_pt, score = gt.support_vector_max(simplex, gt.Vec(direction))
    Tagged(SVector(best_pt))
end

function support_vector_max(mesh::gt.HomogenousMesh{gt.Point{N, T}}, direction, initial_guess::Tagged) where {N,T}
    best_arg, best_value = gt.argmax(x-> dot(SVector(x), direction), gt.vertices(mesh))
    best_vec = SVector(best_arg)
    Tagged(best_vec)
end

any_inside(pt::SVector) = Tagged(pt)
support_vector_max(pt::SVector, direction, initial_guess::Tagged) = Tagged(pt)

function transform_simplex(cache::CollisionCache, poseA, poseB)
    transform_simplex(dimension(typeof(cache)), cache, poseA, poseB)
end

@generated function transform_simplex(::Val{N}, cache::CollisionCache, poseA, poseB) where {N}
    transform_simplex_impl(N, cache, poseA, poseB)
end

function transform_simplex_impl(N, cache, poseA, poseB)
    Expr(:call, :(SVector),
        [:((poseA(value(cache.simplex_points[$i].a)) -
            poseB(value(cache.simplex_points[$i].b)))) for i in 1:(N + 1)]...)
end

# Note: it looks like this can be replaced with transpose(weights) * points in Julia 1.3 (before that, it's a lot slower)
# @generated function linear_combination(weights::StaticVector{N}, points::StaticVector{N}) where {N}
#     expr = :(weights[1] * points[1])
#     for i = 2 : N
#         expr = :($expr + weights[$i] * points[$i])
#     end
#     return quote
#         Base.@_inline_meta
#         $expr
#     end
# end

# Updated for Julia 1.5
function linear_combination(weights::StaticVector{N}, points::StaticVector{N}) where {N}
    return weights' * points
end

struct GJKResult{M, N, T}
    simplex::SVector{M, SVector{N, T}}
    weights::SVector{M, T}
    in_collision::Bool
    closest_point_in_body::Difference{SVector{N, T}, SVector{N, T}}
    nvrtx::Int
    wids::SVector{M, Int}
    iterations::Int
    termination::Int
end

function Base.getproperty(result::GJKResult, sym::Symbol)
    if sym === :signed_distance
        # @warn "The `signed_distance` field was removed. Please use the `separation_distance` and `simplex_penetration_distance` functions."
        if result.in_collision
            return -simplex_penetration_distance(result)
        else
            return separation_distance(result)
        end
    else
        return Core.getfield(result, sym)
    end
end

""" Calculate active simplex points of privileged object for neighbor matching in implicit surface in body frame.  """
function body_simplex(cache::CollisionCache, wids::SVector{M, Int}, nvrtx::Int) where {M}
    simplex_points = [value(cache.simplex_points[i].a) for i in wids[1:nvrtx]]
    # T = typeof(simplex_points[1][2])  # Number type
    return SVector{nvrtx, SVector{3, Float64}}(simplex_points)
end

closest_point_in_world(result::GJKResult) = linear_combination(result.simplex, result.weights)
closest_point_in_body(result::GJKResult) = result.closest_point_in_body
separation_distance(result::GJKResult) = (@assert !result.in_collision; norm(closest_point_in_world(result)))
simplex_penetration_distance(result::GJKResult) = penetration_distance(result.simplex)

"""
Algorithm referenced from Bergen "Collision Detection in Interactive 3D Environments" textbook, page 145.
"""
function gjk!(cache::CollisionCache,
              poseA::Transformation,
              poseB::Transformation,
              max_iter=10,
              eps_tol=1e-16,  # improvement point tolerance
              eps_rel=1e-16)  # collision tolerance

    rotAinv = transform_deriv(inv(poseA), 0)
    rotBinv = transform_deriv(inv(poseB), 0)
    simplex = transform_simplex(cache, poseA, poseB)
    iter = 1

    # Initialize containers
    weights = SVector{4,Float64}(1,0,0,0)
    wids = SVector{4,Int}(1,2,3,4)  # track active simplex points returned from signed_volume (top nvrtx are active)
    prev_wids = wids  # tracks repetition termination condition
    nvrtx = 0  # number of active vertices in simplex
    best_point = simplex[1]  # initialize to first vertex

    while true
        # Compute the support function to find w_k and replace in simplex
        direction = -best_point
        direction_in_A = rotAinv * direction
        direction_in_B = rotBinv * direction

        starting_vertex_index = 1
        starting_vertex = cache.simplex_points[starting_vertex_index]
        starting_vertex_score =
            dot(value(starting_vertex.a), direction_in_A) +
            dot(value(starting_vertex.b), direction_in_B)
        for j in 2:length(cache.simplex_points)
            candidate = cache.simplex_points[j]
            candidate_score =
                dot(value(candidate.a), direction_in_A) +
                dot(value(candidate.b), direction_in_B)
            if candidate_score > starting_vertex_score
                starting_vertex_score = candidate_score
                starting_vertex = candidate
            end
        end

        improved_vertex = Difference(
            support_vector_max(cache.bodyA, direction_in_A, starting_vertex.a),
            support_vector_max(cache.bodyB, -direction_in_B, starting_vertex.b))
        improved_point = poseA(value(improved_vertex.a)) - poseB(value(improved_vertex.b))

        # Primary termination conditions
        v_len = dot(best_point, best_point)
        if iter == max_iter
            # Max iterations
            closest_point_in_body = linear_combination(weights, cache.simplex_points)
            return GJKResult(simplex, weights, false, closest_point_in_body, nvrtx, wids, iter, 1)
        elseif v_len - dot(best_point, improved_point) <= eps_rel*v_len
            # w_k is sufficiently close enough to v_k to not improve v
            closest_point_in_body = linear_combination(weights, cache.simplex_points)
            return GJKResult(simplex, weights, false, closest_point_in_body, nvrtx, wids, iter, 2)
        elseif improved_point in simplex[prev_wids]
            # w_k ∈ W_(k-1) U w_(k-1)
            closest_point_in_body = linear_combination(weights, cache.simplex_points)
            return GJKResult(simplex, weights, false, closest_point_in_body, nvrtx, wids, iter, 3)
        else
            # Add improved point to the simplex set
            nvrtx += 1
            index_to_replace = wids[nvrtx]
            cache.simplex_points[index_to_replace] = improved_vertex
            simplex = setindex(simplex, improved_point, index_to_replace)
            prev_wids = wids[1:nvrtx]
        end

        # Calculate new weights and support set W
        weights, wids, nvrtx = signed_volume(simplex, wids, nvrtx)

        # Secondary termination conditions
        # best_point = linear_combination(weights, simplex)
        # w_set = simplex[wids[1:nvrtx]]
        # w_max = maximum(dot(w_set, w_set))
        #
        # if nvrtx == 4 || dot(best_point, best_point) <= eps_tol*w_max
        #     closest_point_in_body = linear_combination(weights, cache.simplex_points)
        #     return GJKResult(simplex, weights, true, closest_point_in_body, nvrtx, wids, iter, 4)
        # end

        best_point = linear_combination(weights, simplex)
        if nvrtx == 4
            closest_point_in_body = linear_combination(weights, cache.simplex_points)
            return GJKResult(simplex, weights, true, closest_point_in_body, nvrtx, wids, iter, 4)
        end

        iter += 1
    end
end

function penetration_distance(simplex)
    best_dist_squared = nothing
    for i in eachindex(simplex)
        face = simplex_face(simplex, i)
        weights = projection_weights(face)
        closest_point = linear_combination(weights, face)
        dist_squared = closest_point ⋅ closest_point
        if best_dist_squared === nothing || dist_squared < best_dist_squared
            best_dist_squared = dist_squared
        end
    end
    return sqrt(best_dist_squared)
end

function gjk(geomA, geomB,
             poseA::Transformation=IdentityTransformation(),
             poseB::Transformation=IdentityTransformation())
    cache = CollisionCache(geomA, geomB)
    gjk!(cache, poseA, poseB)
end

""" Using ForwardDiff for gradients """
function gjk_diff(q::AbstractVector{T}) where {T}  # q is generalized coordinates
    RotA = composed(RotZ(q[6]), composed(RotY(q[5]), RotX(q[4])))
    TransA = Translation(q[1:3])
    # RotB = RotZ(RotY(RotX(q[4:6])))
    # TransB = Translation(q[1:3])
    poseA = composed(RotA, TransA)
    poseB = IdentityTransformation
    result = gjk(c1, c2, poseA, poseB)
    return separation_distance(result)
end

#---
### Original functions ###
function gjk_original!(cache::CollisionCache,
              poseA::Transformation,
              poseB::Transformation,
              max_iter=100,
              atol=1e-6)
    rotAinv = transform_deriv(inv(poseA), 0)
    rotBinv = transform_deriv(inv(poseB), 0)
    simplex = transform_simplex(cache, poseA, poseB)
    iter = 1

    while true
        weights = projection_weights(simplex)
        min_weight, index_to_replace = findmin(weights)
        best_point = linear_combination(weights, simplex)
        if min_weight > 0
            # Nominally in collision; but check for numerical issues
            separation_squared = best_point ⋅ best_point
            @assert separation_squared < sqrt(1_000_000 * eps(typeof(separation_squared)))
            closest_point_in_body = linear_combination(weights, cache.simplex_points)
            # return GJKResult(simplex, weights, true, closest_point_in_body, iter, simplex, poseA, poseB)
            return GJKResult(simplex, weights, true, closest_point_in_body, 0, SVector{4,Int}(zeros(4,1)), iter, 0)
        end

        direction = -best_point
        direction_in_A = rotAinv * direction
        direction_in_B = rotBinv * direction

        starting_vertex_index = 1
        starting_vertex = cache.simplex_points[starting_vertex_index]
        starting_vertex_score =
            dot(value(starting_vertex.a), direction_in_A) +
            dot(value(starting_vertex.b), direction_in_B)
        for j in 2:length(cache.simplex_points)
            candidate = cache.simplex_points[j]
            candidate_score =
                dot(value(candidate.a), direction_in_A) +
                dot(value(candidate.b), direction_in_B)
            if candidate_score > starting_vertex_score
                starting_vertex_score = candidate_score
                starting_vertex = candidate
            end
        end

        improved_vertex = Difference(
            support_vector_max(cache.bodyA, direction_in_A, starting_vertex.a),
            support_vector_max(cache.bodyB, -direction_in_B, starting_vertex.b))
        improved_point = poseA(value(improved_vertex.a)) - poseB(value(improved_vertex.b))
        score = dot(improved_point, direction)
        if score <= dot(best_point, direction) + atol || iter >= max_iter
            closest_point_in_body = linear_combination(weights, cache.simplex_points)
            return GJKResult(simplex, weights, false, closest_point_in_body, 0, SVector{4,Int}(zeros(4,1)), iter, 0)
        else
            cache.simplex_points[index_to_replace] = improved_vertex
            simplex = setindex(simplex, improved_point, index_to_replace)
            # simplex[index_to_replace] = improved_point
        end

        iter += 1
    end
end

function gjk_original(geomA, geomB,
             poseA::Transformation=IdentityTransformation(),
             poseB::Transformation=IdentityTransformation())
    cache = CollisionCache(geomA, geomB)
    gjk_original!(cache, poseA, poseB)
end


text:
import ..Utils

Base.:*(z::PyCall.PyObject, ps::OpTerm) = typeof(ps)(ps.paulis, ps.coeff * z)
Base.:*(ps::OpTerm, z::PyCall.PyObject) = *(z, ps)

function Utils.isapprox_zero(x::PyCall.PyObject)
    return x == 0
end


text:
using MetabolicSim
using Test

#function to create 1x1 matricies
function mat(sc)
    reshape([sc],1,1)
end

#create parameters
M,N = 1,1
u,Rm = mat(1.0), [1.0]
l,ρ = mat(0.2), [0.2]
p = make_parameters(N,M,u,Rm,l,ρ)

#test vs analytical prediction
C_eq = ρ[1] / Rm[1]
R_eq = -Rm[1] / (u[1] * (l[1] - 1))

@testset begin
    sol = @test_nowarn simulate(ones(N+M),p,t_end = 100)
    @test all(isapprox.(sol[end],[C_eq,R_eq], rtol = 1e-2))
end


text:
"""
    (f::BoundaryFill)(res::AbstractArray{T,2}, verts::Vector{CartesianIndex{2}}, x::Int, y::Int, fill_value::T, boundary_value::T) where {T <: Colorant}

"""
function (f::BoundaryFill)(
    res::AbstractArray{T,2},
    verts::Vector{CartesianIndex{2}},
    x::Int,
    y::Int,
    ) where {T<:Colorant}
    if checkbounds(Bool, res, y, x)
        if (res[y, x] != f.boundary_value && res[y, x] != f.fill_value)
            res[y, x] = f.fill_value
            f(res, verts, x + 1, y)
            f(res, verts, x, y + 1)
            f(res, verts, x - 1, y)
            f(res, verts, x, y - 1)
        end
    end
    res
end

"""
    draw!(img::AbstractArray{T,2}, verts::Vector{CartesianIndex{2}}, f::AbstractPolyFillAlgorithm; closed::Bool)

Draw on `img` using algorithm `f`.

# Output

When `img` is specified, changes are made on `img` and returned.

# Example

Just simply pass an algorithm with parameters, with image and vertices of polygon

```julia
using ImageDraw

img = zeros(RGB, 7, 7)
expected = copy(img)
expected[2:6, 2:6] .= RGB{N0f8}(1)

verts = [CartesianIndex(2, 2), CartesianIndex(2, 6), CartesianIndex(6, 6), CartesianIndex(6, 2), CartesianIndex(2,2)]

draw!(img, verts, BoundaryFill(4, 4; fill_value = RGB(1), boundary_value = RGB(1)); closed = true)
```
"""
function draw!(
    img::AbstractArray{<:Colorant,N},
    verts::Vector{CartesianIndex{N}},
    f::AbstractPolyFillAlgorithm;
    closed::Bool = false,
    ) where N

    # Boundscheck for verts
    @boundscheck map(x -> (checkbounds(Bool, img, x[1], x[2]) ? nothing : throw(DomainError(x, "Vertice $x outside the image array domain"))), verts)
    # Boundscheck for seed point
    @boundscheck checkbounds(img, f.x, f.y)

    if (closed == true)
        for i = 1:length(verts)-1
            draw!(img, LineSegment(verts[i], verts[i+1]), f.fill_value)
        end
    end
    f(img, verts, f.x, f.y)
end

"""
    draw(img::AbstractArray{T,2}, verts::Vector{CartesianIndex{2}}, f::AbstractPolyFillAlgorithm; closed::Bool)

Draw on `img` using algorithm `f`.

# Output

When `img` is specified, a copy of `img` is made and changes are made on it and returned.

# Example

Just simply pass an algorithm with parameters, with image and vertices of polygon

```julia
using ImageDraw

img = zeros(RGB, 7, 7)
expected = copy(img)
expected[2:6, 2:6] .= RGB{N0f8}(1)

verts = [CartesianIndex(2, 2), CartesianIndex(2, 6), CartesianIndex(6, 6), CartesianIndex(6, 2), CartesianIndex(2,2)]

res = draw(img, verts, BoundaryFill(4, 4; fill_value = RGB(1), boundary_value = RGB(1)); closed = true)
```
"""
function draw(
    img::AbstractArray{<:Colorant,N},
    verts::Vector{CartesianIndex{N}},
    f::AbstractPolyFillAlgorithm;
    kwargs...
    ) where N
    draw!(copy(img), verts, f; kwargs...)
end


text:
struct tikzObj

    flname::String

    lw::Float64                  # default line width
    linecolor::Vector{Float64}   # default line color
    axcolor::Vector{Float64}     # axis color

    standalone::Bool
    smooth::Bool

    fID::IOStream

end

function tikzObj( flname::String;
                  lw = 0.75,
                  linecolor = [0., 0., 0.],
                  axcolor = [0.1, 0.1, 0.1],
                  smooth = false,
                  standalone = false,
                  scale = 1.0,
                  colors=Dict{String,Vector{Float64}}(),
                  styles=Vector{tikzStyle}( 0 ),
                  textheight=12,
                  additional_packages::Vector{String} = Vector{String}(undef,0),
                  additional_packages_with_options::Vector{Tuple{String,String}} = Vector{Tuple{String,String}}(undef,0) )

    fID = open( flname, "w" )

    if standalone
        @printf( fID, "\\documentclass[%dpt]{standalone}\n", textheight )
        # Signature
        @printf( fID, "%% Generated using Tikz.jl on %s, %s\n\n", Dates.today(), Dates.format(Dates.now(), "HH:MM")  )

        # Packages
        @printf( fID, "\\usepackage{tikz,pgfplots}\n" )
        @printf( fID, "\\usepackage{xcolor}\n" )

        for package in additional_packages
            @printf( fID, "\\usepackage{%s}\n", package )
        end

        for package in additional_packages_with_options
            @printf( fID, "\\usepackage[%s]{%s}\n", package[1], package[2] )
        end

        @printf( fID, "\n" )

    end

    # Colors
    @printf( fID, "\\definecolor{defaultc}{RGB}{%3.0f, %3.0f, %3.0f}\n\n", linecolor[1], linecolor[2], linecolor[3] )
    coln = collect(keys(colors))
    for jj in 1:length(colors)
        @printf( fID, "\\definecolor{%s}{RGB}{%3.0f, %3.0f, %3.0f}\n\n",
                 coln[jj], colors[coln[jj]][1], colors[coln[jj]][2], colors[coln[jj]][3] )
    end

    # Styles
    for jj in 1:length(styles)
        #\tikzstyle{my style}=[<tikz options>]
        @printf( fID, "\\tikzstyle{%s}=[", styles[jj].name )
        if styles[jj].lw > 0.0
            @printf( fID, "line width=%3.2fpt, ", styles[jj].lw )
        end
        if length(styles[jj].color) > 0
            @printf( fID, "%s, ", styles[jj].color )
        end
        if styles[jj].linestyle == "--"
            @printf( fID, "dashed, " )
        elseif length(styles[jj].linestyle) > 0
            @printf( fID, "%s, ", styles[jj].linestyle )
        end
        @printf( fID, "%s]\n ", styles[jj].extracommands )
    end
    @printf( fID, "\n" )

    if standalone
        # Start document
        @printf( fID, "\\begin{document}\n\n" )
    end

    @printf( fID, "\\begin{tikzpicture}")
    @printf( fID, "[color=defaultc,line width=%3.2fpt, scale=%3.2f]\n\n", lw, scale )

    return tikzObj( flname, lw, linecolor, axcolor, standalone, smooth, fID )

end

function closeObj( object::tikzObj )

    @printf( object.fID, "\\end{tikzpicture}")

    if object.standalone
        @printf( object.fID, "\n\n\\end{document}")
    end

    close( object.fID )

end

max_stars_repo_name: mopg/julia-tikz

max_issues_repo_name: mopg/julia-tikz

max_forks_repo_name: mopg/julia-tikz


text:
#= TicTacToe implementation in Julia
@author : soumitradev (Soumitra Shewale)
=#

# Create 3x3 board for TicTacToe
board = [[" ", " ", " "],
         [" ", " ", " "],
         [" ", " ", " "]];

# Define Player Characters for singleplayer
# Spelled wrong for (bad) comedic purposes (I know, very funny haha)
hooman = "X"
ai = "O"

# Define win codes used in checking who won, and in minimax algorithm as a score
# -1 indicates undesirable outcome in minimax, while +1 indicates desirable, and 0 is neutral
hooman_win = -1
ai_win = 1
tie = 0

# Define Player Characters for multiplayer
player1 = "X"
player2 = "O"

# Define win codes used in checking who won
# These codes are not used in minimax, but are only used to denote the player
player1_win = -1
player2_win = 1
tie = 0

# Get array of tuples of all available (unoccupied) indices in board to allow AI to choose from these options
function get_available(board)
    available_pos = []
    # Go through all elements in array and check if they are unoccupied
    for i in 1:3, j in 1:3
        # Check if element is unoccupied
        if board[i][j] == " "
            # If the element is unoccupied, append the tuple of the indices to the available array.
            push!(available_pos, (i, j))
        end
    end
    return available_pos
end

# Check number of available (unoccupied) cells similar to how we did above.
function check_available(board)
    available = 0
    for i in 1:3, j in 1:3
        if board[i][j] == " "
            available += 1
        end
    end
    return available
end

# Check if the board is in a terminal (win or tie) situation.
function check_win(board)
    # Define p1 & p2 and their properties to avoid mixing up with player1, player2, hooman, and ai variables.
    p1 = "X"
    p2 = "O"

    p1_win = -1
    p2_win = +1
    tie = 0

    # Horizontal line checker
    for i in 1:3
        # For every row, check if the elements in a row are same, and belong to some player (instead of being unoccupied)
        if ((board[i][1] == board[i][2] == board[i][3]) && (board[i][1] in [p1, p2]))
            # If row belongs to p1, return signal showing p1 wins.
            if (board[i][1] == board[i][2] == board[i][3] == p1)
                return p1_win
            # If row belongs to p2, return signal showing p2 wins.
            elseif (board[i][1] == board[i][2] == board[i][3] == p2)
                return p2_win
            end
        end
    end

    # Vertical line checker
    for i in 1:3
        # For every column, check if the elements in a column are same, and belong to some player (instead of being unoccupied)
        if ((board[1][i] == board[2][i] == board[3][i]) && (board[1][i] in [p1, p2]))
            # If column belongs to p1, return signal showing p1 wins.
            if (board[1][i] == board[2][i] == board[3][i] == p1)
                return p1_win
            # If column belongs to p2, return signal showing p2 wins.
            elseif (board[1][i] == board[2][i] == board[3][i] == p2)
                return p2_win
            end
        end
    end

    # Diagonal (from (1, 1) to (3, 3)) line checker
    if (board[1][1] == board[2][2] == board[3][3])
        # Check if diagonal belongs to p1 and if true, return signal showing p1 wins
        if (board[1][1] == board[2][2] == board[3][3] == p1)
            return p1_win
        # Check if diagonal belongs to p2 and if true, return signal showing p2 wins
        elseif (board[1][1] == board[2][2] == board[3][3] == p2)
            return p2_win
        end
    end

    # Diagonal (from (3, 1) to (1, 3)) line checker
    if (board[1][3] == board[2][2] == board[3][1])
        # Check if diagonal belongs to p1 and if true, return signal showing p1 wins
        if (board[1][3] == board[2][2] == board[3][1] == p1)
            return p1_win
        # Check if diagonal belongs to p2 and if true, return signal showing p2 wins
        elseif (board[1][3] == board[2][2] == board[3][1] == p2)
            return p2_win
        end
    end

    # If no one won, check if available cells are 0. If they are 0, it is a tie. Else, the game hasn't ended yet, so return nothing
    if check_available(board) == 0
        return tie
    else
        return nothing
    end
end

# Use minimax algorithm to determine best/worst step
# Board is the game board
# depth is the current recursion depth the algorithm is checking after
# is_maximising shows if the player that just made the move wants to maximise his/her winning chances.

# NOTE: This can be optimised further using alpha-beta pruning to decrease computational complexity
#       and using depth as a factor in deciding score to use fastest way to win, but that's for later.
function minimax(board, depth, is_maximising)
    # Check if game is over (base case)
    result = check_win(board)

    # If game is over, return the score (the win code in this case, which also is the score in that case)
    if result !== nothing
        return result
    # That finishes our base case for recursion. Phew!

    # If player that just made move wants to maximise chances of winning, find best subcase in case of the given move
    elseif is_maximising

        # Initialise best_score at -Inf so every number is greater than that, and the best score is replaced appropriately
        best_score = -Inf

        # Check every vacant spot, place player character in the spot (or cell) we are evaluating and get the best subcase.
        for i in 1:3, j in 1:3
            if board[i][j] == " "

                # Place player character
                board[i][j] = ai
                # Evaluate best subcase (this time is_maximising is false because it is the user's turn and he wants to minimise chances of AI winning)
                # Even if Hooman doesn't play his/her best move, computer still benefits since it expected a much harder case.
                score = minimax(board, depth + 1, false)

                # Reverse temporary move we made
                board[i][j] = " "

                # If score is better, use this one
                best_score = max(score, best_score)
            end
        end
        # Return the best score so the best case is evaluated.
        return best_score
    else
        # Initialise best_score at Inf so every number is greater than that, and the best score is replaced appropriately
        # This time around, we are simulating the hooman move by minimising the chances of AI winning.
        best_score = Inf

        # Check every vacant spot, place player character in the spot (or cell) we are evaluating and get the best subcase (best for hooman, worst for AI).
        for i in 1:3, j in 1:3
            if board[i][j] == " "

                # Place player character
                board[i][j] = hooman
                # Evaluate best subcase (this time is_maximising is true because it is the AI's turn and he wants to maximise chances of AI winning)
                score = minimax(board, depth + 1, true)

                # Reverse temporary move we made
                board[i][j] = " "

                # If score is better, use this one
                best_score = min(score, best_score)
            end
        end
        # Return the best score so the best case is evaluated.
        return best_score
    end
end

# Calculate best move
function best_move(board)
    # Use maximising part of minimax at start to maximise current player's winning chances
    best_score = -Inf

    # Initialise best_pos (index of best move)
    best_pos = (0, 0)

    # Check for vacant spaces, and get best case using minimax.
    for i in 1:3, j in 1:3
        if board[i][j] == " "

            # Place temporary player character
            board[i][j] = ai

            # Evaluate if this is the best move to make using minimax
            score = minimax(board, 0, false)

            # Reverse the temporary changes we made
            board[i][j] = " "

            # If score is better than what we have right now, replace the best score, and set the best_move.
            if score > best_score
                best_score = score
                best_pos = (i, j)
            end
        end
    end

    # Return best position
    return best_pos
end

# Calculate worst move for loser mode
function worst_move(board)

    # Use minimising part of minimax at start to minimise current player's winning chances
    worst_score = Inf

    # Initialise worst_pos (index of worst move)
    worst_pos = (0, 0)

    # Check for vacant spaces, and get worst case using minimax.
    for i in 1:3, j in 1:3

        if board[i][j] == " "
            board[i][j] = ai

            # Minimise AI = Maximise Hooman. So, is_maximising is true for Hooman.
            score = minimax(board, 0, true)
            board[i][j] = " "

            # Get worst score and position of move similar to how we did in best_move()
            if score < worst_score
                worst_score = score
                worst_pos = (i, j)
            end
        end
    end

    # Return worst position
    return worst_pos
end

# AI descision tree according to difficulty level selected by Hooman.
function ai_move!(board, difficulty)

    # If difficulty selected is `loser`, make the AI perform the worst move
    if difficulty == "loser"
        move_pos = worst_move(board)
        make_move!(board, ai, move_pos)

    # If difficulty selected is `easy`, make the AI perform random moves all the time.
    elseif difficulty == "easy"
        if check_available(board) != 0
            move_pos = rand(get_available(board))
            make_move!(board, ai, move_pos)
        end

    # If difficulty selected is `medium`, make the AI perform a random move half the time, and the best move the other half of the time
    elseif difficulty == "medium"
        # Use rand to decide if to make random move or smart move
        if rand() <= 0.5
            if check_available(board) != 0
                move_pos = rand(get_available(board))
                make_move!(board, ai, move_pos)
            end
        else
            move_pos = best_move(board)
            make_move!(board, ai, move_pos)
        end

    # If difficulty selected is `hard`, make the AI perform a random move a quarter the time, and the best move the rest of the time
    elseif difficulty == "hard"
        # Use rand to decide if to make random move or smart move
        if rand() <= 0.25
            if check_available(board) != 0
                move_pos = rand(get_available(board))
                make_move!(board, ai, move_pos)
            end
        else
            move_pos = best_move(board)
            make_move!(board, ai, move_pos)
        end
    else
        move_pos = best_move(board)
        make_move!(board, ai, move_pos)
    end
end

# Function to print the board with row and column labels
function print_board(board)
    println("  1   2   3")
    println("1 ", join(board[1], " | "))
    println(" ———┼———┼———")
    println("2 ", join(board[2], " | "))
    println(" ———┼———┼———")
    println("3 ", join(board[3], " | "))
end

# Function to make a move for a player at (x, y).
# Assumes that position is not taken. So, moves are checked before being made
# If cell is taken, it just ignores the move.
function make_move!(board, player, (x, y))
    if board[x][y] == " "
        board[x][y] = player
    end
end

# Create function to get user input and recursively validate the user input
function move_input(board, text_to_print, in_same_line, required_params)

    # If we want to print on the same line, use `print()`. Else, `use println()`
    in_same_line ? print(text_to_print) : println(text_to_print)

    # Get input and split it into an array
    txtin = split(lowercase(readline()));

    # Allow users to exit
    if txtin == ["exit"]
        exit()
    end

    # If the required amount of parameters is not met or if indices are out of array, or if place is already taken, ask for input again
    if length(txtin) < required_params || !(txtin[1] in ["1", "2", "3"]) || !(txtin[2] in ["1", "2", "3"]) || board[parse(Int, txtin[1])][parse(Int, txtin[2])] != " "
        println("Please enter a valid command (see README)")
        move_input(board, text_to_print, in_same_line, required_params)
    else
        # If input is valid, return the input as an array of parameters passed to it..
        return txtin
    end
end

# Create function to get user input and recursively validate the user input
function generalised_input(text_to_print, in_same_line, allowed_options, required_params)

    # If we want to print on the same line, use `print()`. Else, `use println()`
    in_same_line ? print(text_to_print) : println(text_to_print)

    # Get input and split it into an array
    txtin = split(lowercase(readline()));

    # Allow user to exit
    if txtin == ["exit"]
        exit()
    end

    # If the required amount of parameters is not met or if the command is not in the list of allowed options,
    # ask for input again
    if length(txtin) < required_params || !(txtin[1] in allowed_options)
        println("Please enter a valid command (see README)")
        generalised_input(text_to_print, in_same_line, allowed_options, required_params)
    else
        # If input is valid, return the input as a string.
        return join(txtin, " ")
    end
end

# Ask if user wants to play again
function play_again()

    # Get user input for whether he/she wants to play again
    ans = generalised_input("Play again? (y/n): ", true, ["y", "n"], 1)

    # If user wants to play again, run the game again. Else, exit()
    if ans == "y"
        # Run this same script in a new instance, and exit this instance of this script
        run(`julia $PROGRAM_FILE`)
        exit()
    else
        # Exit
        exit()
    end
end

# Start multiplayer game
function start_multiplayer_game()
    # Print starting text
    println("Player 1 (X) goes first. Player 2 (O) goes after Player 1. \nType `exit` to exit the game whenever you want.")

    while true
        # Print current state of board
        print_board(board)

        # Get player 1 input for cell to play in after input scrubbing
        txt = move_input(board, "Player 1, enter row and column seperated by space: ", true, 2)

        # Perform move player 1 requested.
        make_move!(board, player1, (parse(Int, txt[1]), parse(Int, txt[2])))

        # Check for win, and print appropriate text.
        if (check_win(board) in [tie, player1_win, player2_win])

            # Print board after le epic ending move
            print_board(board)

            # Print appropriate text
            if check_win(board) == player1_win
                println("Player 1 Wins!")
            elseif check_win(board) == player2_win
                println("Player 2 Wins!")
            else
                println("Tie.")
            end

            # Break out of while loop to end game
            break
        end

        # If no one won yet, print the board and move to Player 2's move
        print_board(board)

        # Get player 2 input for cell to play in after input scrubbing
        txt = move_input(board, "Player 2, enter row and column seperated by space: ", true, 2)

        # Perform move player 2 requested.
        make_move!(board, player2, (parse(Int, txt[1]), parse(Int, txt[2])))

        # Check for win, and print appropriate text.
        if (check_win(board) in [tie, hooman_win, ai_win])

            # Print board after le epic ending move
            print_board(board)

            # Print appropriate text
            if check_win(board) == hooman_win
                println("Hooman wins. There is hope! The hoomans are still worthy!")
            elseif check_win(board) == ai_win
                println("AI wins. AI is going to take over. We're practically dead.")
            else
                println("Tie. Have we reached the singularity?")
            end

            # Break out of while loop to end game
            break
        end
    end

    # After ending game, ask if user wants to play again.
    play_again()
end

# Start multiplayer game
function start_singleplayer_game()
    # Print difficulty select
    difficulty = """
    You will be playing as X. The computer is O.
    Type `exit` to exit the game whenever you want.

    Select difficulty:
    - `loser`: AI always makes worst move
    - `easy`: AI plays randomly
    - `medium`: AI plays smart, but still makes mistakes.
    - `hard`: AI plays smarter, makes less mistakes.
    - `impossible`: AI always makes best move. Impossible to win this.
    """

    # Get difficulty selected after user input scrubbing
    diff_selected = generalised_input(difficulty, true, ["loser", "easy", "medium", "hard", "impossible"], 1)

    # Start singleplayer game
    while true

        # Print current state of board
        print_board(board)

        # Get hooman input after scrubbing
        txt = move_input(board, "Enter row and column seperated by space: ", true, 2)

        # Perform hooman move
        make_move!(board, hooman, (parse(Int, txt[1]), parse(Int, txt[2])))

        # Check for win, and print appropriate text.
        if (check_win(board) in [tie, hooman_win, ai_win])

            # Print current state of board after le epic final move
            print_board(board)

            # Print appropriate text after end of game
            if check_win(board) == hooman_win
                println("Hooman wins. There is hope! The hoomans are still worthy!")
            elseif check_win(board) == ai_win
                println("AI wins. AI is going to take over. We're practically dead.")
            else
                println("Tie. Have we reached the singularity?")
            end

            # Break out of while loop to end game.
            break
        end

        # If no one won yet, allow AI to make move according to difficulty selected.
        ai_move!(board, diff_selected)

        # Check for win, and print appropriate text.
        if (check_win(board) in [tie, hooman_win, ai_win])

            # Print current state of board after le epic final move
            print_board(board)

            # Print appropriate text after end of game
            if check_win(board) == hooman_win
                println("Hooman wins. There is hope! The hoomans are still worthy!")
            elseif check_win(board) == ai_win
                println("AI wins. AI is going to take over. We're practically dead.")
            else
                println("Tie. Have we reached the singularity?")
            end

            # Break out of while loop to end game.
            break
        end
    end

    # Ask if hooman wants to play again.
    play_again()
end

# Start game and ask for choice of game mode
mode = generalised_input("Welcome to TicTacToe! Would you like to play singleplayer (1) or multiplayer (2)?", false, ["1", "2"], 1)

# Check mode after scrubbing and start appropriate game mode.
if mode == "1"
    start_singleplayer_game()
else
    start_multiplayer_game()
end

max_stars_repo_name: soumitradev/TicTacToe

max_issues_repo_name: soumitradev/TicTacToe

max_forks_repo_name: soumitradev/TicTacToe


text:
using Test
using DataFrames
using CSV
using ChemistryFeaturization.Utils.ElementFeatureUtils
using ..ChemistryFeaturization.Codec: build_onehot_vec

@testset "ElementFeatureUtils" begin
    df = CSV.read(abspath(@__DIR__, "..", "test_data", "lookup_table.csv"), DataFrame)

    # fea_minmax
    @test_throws AssertionError fea_minmax("heffalump")
    @test fea_minmax("Group") == [1, 18]
    @test fea_minmax("MeaningOfLife", df) == [-1, 42]

    # default_log
    @test default_log("Block") == false # not numbers
    @test default_log("MeaningOfLife", df) == false # values span 0
    @test default_log("Valence") == false # extremal value 0
    @test default_log("Atomic mass") == true
    @test default_log("Atomic mass", threshold = 3) == false

    # default_categorical
    @test default_categorical("Block") == true
    @test get_param_vec(false, 2) == [false, false]
    @test get_param_vec([false, false], 2) == [false, false]
    @test get_param_vec([false, false], 3) == [false, false, false]
    @test get_param_vec([false, false, false], 2) == [false, false]
    @test get_param_vec([3, 2], 3, pad_val = 0) == [3, 2, 0]
    @test get_param_vec(['a'], 3, pad_val = 'b') == ['a', 'b', 'b']
    @test get_param_vec('a', 3, pad_val = 'b') == ['a', 'a', 'a']

    @test get_param_vec([false, false], 3) == [false, false, false]
    @test get_param_vec([false, false, false], 2) == [false, false]
    @test get_param_vec([3, 2], 3, pad_val = 0) == [3, 2, 0]
    @test get_param_vec(['a'], 3, pad_val = 'b') == ['a', 'b', 'b']
    @test get_param_vec('a', 3, pad_val = 'b') == ['a', 'a', 'a']

    # get_bins
    @test get_bins("Block") == ["s", "p", "d", "f"]
    @test get_bins("first_letter", df) == ["A", "C", "T"]
    @test get_bins("Group", categorical = false, nbins = 2) == 1.0:8.5:18.0
    @test get_bins("MeaningOfLife", df, categorical = true) == [-1, 0, 42]
    logbins = get_bins("Atomic mass")
    @test isapprox(logbins[2] / logbins[1], logbins[end] / logbins[end-1])
    neglogbins = get_bins("neg_nums", df, logspaced = true)
    @test neglogbins[1] == -1000
    # errors and things that should be ignored...
    @test get_bins("Block", nbins = 3) == ["s", "p", "d", "f"]
    @test get_bins("Valence", logspaced = true) == 0:1:13
    @test_throws AssertionError get_bins("Valence", categorical = false, logspaced = true)
    @test get_bins("noarsenic", df) == 1.0:0.1:2.0
    @test get_bins("noarsenic", df, categorical = true) == [1, 2]

    # build_onehot_vec
    @test build_onehot_vec("s", ["s", "p", "d", "f"], true) == [1, 0, 0, 0]
    @test_throws AssertionError build_onehot_vec("s", ["s", "p", "d", "f"], false)
    bins = 0:2:6
    @test build_onehot_vec(3, bins, false) == [0, 1, 0]
    @test build_onehot_vec(0, bins, false) == [1, 0, 0]
    @test build_onehot_vec(2, bins, false) == [0, 1, 0]
    @test build_onehot_vec(6, bins, false) == [0, 0, 1]
    @test_throws AssertionError build_onehot_vec(-1, bins, false)
    @test_throws ArgumentError build_onehot_vec(1, bins, true)

    # onehot_lookup_encoder
    @test onehot_lookup_encoder("C", "Block") == [0, 1, 0, 0]
    @test onehot_lookup_encoder("C", "Row")[2] == 1
    @test onehot_lookup_encoder("C", "MeaningOfLife", df, nbins = 4) == [0, 0, 0, 1]
    @test_throws TypeError onehot_lookup_encoder("He", "X")
    @test_throws TypeError onehot_lookup_encoder("As", "noarsenic", df)
    @test_throws AssertionError onehot_lookup_encoder("C", "heffalump")
    @test_throws AssertionError onehot_lookup_encoder("Si", "MeaningOfLife", df)

    # onecold_decoder
    @test onecold_decoder([0, 1, 0, 0], "Block") == "p"
    @test onecold_decoder([0 1; 1 0; 0 0; 0 0], "Block") == ["p", "s"]
    decoded = onecold_decoder(
        onehot_lookup_encoder("As", "MeaningOfLife", df),
        "MeaningOfLife",
        df,
    )
    @test decoded[1] <= 0 < decoded[2]
    @test_throws AssertionError onecold_decoder([1, 0, 0], "heffalump")

end


text:
@registerN Laguerre AbstractCCOP1 α
export Laguerre

"""
   Laguerre{α, T <: Number}

The  Laguerre polynomials have weight function `x^α * exp(-x)` over the domain `[0, oo)`. The parameter `α` is specified through the constructor.

```jldoctest
julia> using Polynomials, SpecialPolynomials

julia> p = Laguerre{1/2}([1,2,3])
typename(Laguerre){0.5}(1⋅Lᵅ₀(x) + 2⋅Lᵅ₁(x) + 3⋅Lᵅ₂(x))

julia> convert(Polynomial, p)
Polynomials.Polynomial(9.625 - 9.5*x + 1.5*x^2)
```

The Laguerre polynomials are the case `α=0`.

```jldoctest
julia> using Polynomials, SpecialPolynomials

julia> p = Laguerre{0}([1,2,3])
typename(Laguerre){0}(1⋅L₀(x) + 2⋅L₁(x) + 3⋅L₂(x))

julia> convert(Polynomial, p)
Polynomials.Polynomial(6.0 - 8.0*x + 1.5*x^2)

julia> phi(u, i) = derivative(u) -  u # verify Rodrigues' formula for small n; n! L_n = (d/dx-1)^n x^n
phi (generic function with 1 method)

julia> x = Polynomial(:x)
Polynomials.Polynomial(1.0*x)

julia> n = 7
7

julia> factorial(n) * basis(Laguerre{0}, n) - foldl(phi, 1:n, init=x^n)
Polynomials.Polynomial(-5.4569682106375694e-12 + 1.4551915228366852e-11*x - 7.275957614183426e-12*x^2)
```

"""
Laguerre

basis_symbol(::Type{<:Laguerre{α}}) where {α} = "Lᵅ"
basis_symbol(::Type{<:Laguerre{0}}) = "L"
Polynomials.domain(::Type{<:Laguerre}) = Polynomials.Interval(0, Inf)

abcde(::Type{<:Laguerre{α}}) where {α} =
    NamedTuple{(:a, :b, :c, :d, :e)}((0, 1, 0, -1, α + 1))

k0(P::Type{<:Laguerre}) = one(eltype(P))
k1k0(P::Type{<:Laguerre{α}}, k::Int) where {α} = -one(eltype(P)) / (k + 1) # k >=0

function norm2(::Type{<:Laguerre{α}}, n) where {α}
    iszero(α) && return one(α)
    Γ(n + α + 1) / Γ(n + 1)
end

function ω₁₀(::Type{<:Laguerre{α}}, n) where {α}
    iszero(α) && return 1.0
    val = Γ(n + 1 + α + 1) / Γ(n + α + 1)
    val /= n + 1 # Γ(n+1)/Γ(n+2)
    sqrt(val)
end

weight_function(::Type{<:Laguerre{α}}) where {α} = x -> x^α * exp(-x)
generating_function(::Type{<:Laguerre{α}}) where {α} =
    (t, x) -> begin
        exp(-t * x / (1 - t)) / (1 - t)^(α + 1)
    end
function classical_hypergeometric(::Type{<:Laguerre{α}}, n, x) where {α}
    α > -1 || throw(ArgumentError("α > -1 is required"))
    as = -n
    bs = α + 1
    Pochhammer_factorial(α + 1, n) * pFq(as, bs, x)
end

gauss_nodes_weights(p::Type{P}, n) where {α,P<:Laguerre{α}} =
    FastGaussQuadrature.gausslaguerre(n, α)

## Overrides

# default  connection between Laguerre is popping out  0s
function Base.iterate(
    o::Connection{P,Q},
    state=nothing,
) where {β,P<:Laguerre{β},α,Q<:Laguerre{α}}
    k, n = o.k, o.n

    if state == nothing
        i = k
        i > n && return nothing
    elseif state + 1 > n
        return nothing
    else
        i = state + 1
    end

    # α(i,k)
    K = i - k
    N = α - β - 1 + K

    return (i, generalized_binomial(N, K)), i
end

##
## --------------------------------------------------
##
@registerN MonicLaguerre AbstractCCOP1 α
export MonicLaguerre
ϟ(::Type{<:MonicLaguerre{α}}) where {α} = Laguerre{α}
ϟ(::Type{<:MonicLaguerre{α,T}}) where {α,T} = Laguerre{α,T}
@register_monic(MonicLaguerre)

@registerN OrthonormalLaguerre AbstractCCOP1 α
export OrthonormalLaguerre
ϟ(::Type{<:OrthonormalLaguerre{α}}) where {α} = Laguerre{α}
ϟ(::Type{<:OrthonormalLaguerre{α,T}}) where {α,T} = Laguerre{α,T}
@register_orthonormal(OrthonormalLaguerre)


text:
struct Report_Folder
    partition_id::String
    partition_report::Partition_Report
    replicate_report::Replicate_Report
    chain_reports::Dict{Chain_ID,Chain_Report}
end

function generate_reports(chains::Dict{Chain_ID,Vector{EM_step}}, test_sets::Dict{String,Vector{LongSequence{DNAAlphabet{2}}}})
    #check chains dict for problems

    length(chains)==0 && throw(ArgumentError("Empty chains dict!"))
    any(chain->length(chain)<2, values(chains)) && throw(ArgumentError("Some chains are too short (<2 EM steps)! Probably not all chains have been operated on by EM workers yet. Try EM_converge first!"))
    unconv=sum(!chain[end].converged for chain in values(chains))
    unconv > 0 && @warn "Not all chains are converged to the selected step delta."
    length(test_sets)==0 && throw(ArgumentError("Empty test_sets dict!"))

    chain_reports = report_chains(chains, test_sets)
    partition_reports = report_partitions(chain_reports)

    report_folders=Dict{String, Report_Folder}()

    for part_report in partition_reports
        rep_report=report_replicates(part_report.best_repset, chains)
        chain_subset=Dict{Chain_ID,Chain_Report}()
        for id in keys(chain_reports)
            id.obs_id==part_report.partition_id && (chain_subset[id]=chain_reports[id])
        end
        report_folders[part_report.partition_id]=Report_Folder(part_report.partition_id, part_report, rep_report, chain_subset)
    end

    return report_folders
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)
#

# Dummy data
struct EmptyObject3Ddata <: Modia3D.AbstractObject3Ddata
end
const emptyObject3DData = EmptyObject3Ddata()

struct EmptyTwoObject3DObject <: Modia3D.AbstractTwoObject3DObject
end
const emptyTwoObject3DObject = EmptyTwoObject3DObject()

# FixedJoint is used if an Object3D is rigidly fixed to its parent Object3D
struct FixedJoint <: Modia3D.AbstractJoint
end
fixedJoint = FixedJoint()

# FreeMotion is used if an Object3D is moving freely with respect to its parent Object3D
mutable struct FreeMotion <: Modia3D.AbstractJoint
   canCollide::Bool                        # = true, if collision detection is enabled between the two Object3Ds connected by the FreeMotion joint (default = true)
   isDriven::Bool                          # = true, if joint is driven; = false, if generalized coordinates are used as states

   # Definitions special to FreeMotions joints
   r::ModiaMath.RealSVector3
   v::ModiaMath.RealSVector3
   a::ModiaMath.RealSVector3
   q::ModiaMath.RealSVector{4}
   derq::ModiaMath.RealSVector{4}
   w::ModiaMath.RealSVector3
   z::ModiaMath.RealSVector3

   residue_w::ModiaMath.RealSVector3
   residue_f::ModiaMath.RealSVector3
   residue_t::ModiaMath.RealSVector3
   residue_q::ModiaMath.RealScalar

   function FreeMotion(obj::Modia3D.AbstractAssemblyComponent;
                       r_start::AbstractVector = ModiaMath.ZeroVector3D,
                       v_start::AbstractVector = ModiaMath.ZeroVector3D,
                       q_start::AbstractVector = ModiaMath.NullQuaternion,
                       w_start::AbstractVector = ModiaMath.ZeroVector3D)

      r    = ModiaMath.RealSVector3(  :r   , numericType=ModiaMath.XD_EXP                ,             unit="m"    , start=r_start, fixed=true, nominal=1.0, analysis=ModiaMath.OnlyDynamicAnalysis, info="Relative position vector from origin of obj 1 to origin of obj 2, resolved in obj1")
      v    = ModiaMath.RealSVector3(  :v   , numericType=ModiaMath.XD_IMP                , integral=r, unit="m/s"  , start=v_start, fixed=true, nominal=1.0, analysis=ModiaMath.OnlyDynamicAnalysis, info="der(r): relative velocity of origin of obj 2 with respect to origin of obj 1, resolved in obj 1")
      a    = ModiaMath.RealSVector3(  :a   , numericType=ModiaMath.DER_XD_IMP            , integral=v, unit="m/s^2",                                                                     info="der(v): relative acceleration of origin of obj 2 with respect to origin of obj 1, resolved in obj 1")

      q    = ModiaMath.RealSVector{4}(:q   , numericType=ModiaMath.XD_IMP                ,               start=q_start, fixed=false, nominal=1.0,                            info="Relative quaternion to rotate obj 1 into obj 2")
      derq = ModiaMath.RealSVector{4}(:derq, numericType=ModiaMath.DER_XD_IMP, integral=q, unit="1/s",                               analysis=ModiaMath.OnlyDynamicAnalysis, info="der(q)")
      w    = ModiaMath.RealSVector3(  :w   , numericType=ModiaMath.XD_IMP                , unit="rad/s", start=w_start, fixed=true , nominal=1.0, analysis=ModiaMath.OnlyDynamicAnalysis, info="Relative angular velocity of obj 2 with respect to obj 1, resolved in obj 2")
      z    = ModiaMath.RealSVector3(  :z   , numericType=ModiaMath.DER_XD_IMP, integral=w, unit="rad/s^2",                           analysis=ModiaMath.OnlyDynamicAnalysis, info="der(w): Relative angular acceleration of obj 2 with respect to obj 1, resolved in obj 2")

      residue_w = ModiaMath.RealSVector3(:residue_w, numericType=ModiaMath.FD_IMP, analysis=ModiaMath.OnlyDynamicAnalysis, info="Angular velocity residue")
      residue_f = ModiaMath.RealSVector3(:residue_f, numericType=ModiaMath.FD_IMP, analysis=ModiaMath.OnlyDynamicAnalysis, info="Momentum equation residue")
      residue_t = ModiaMath.RealSVector3(:residue_t, numericType=ModiaMath.FD_IMP, analysis=ModiaMath.OnlyDynamicAnalysis, info="Angular momentum equation residue")
      residue_q = ModiaMath.RealScalar(  :residue_q, numericType=ModiaMath.FC    ,                                         info="Quaternion constraint residue")

      joint = new(true, false, r, v, a, q, derq, w, z, residue_w, residue_f, residue_t, residue_q)
      joint.r._internal.within = obj
      joint.v._internal.within = obj
      joint.a._internal.within = obj
      joint.q._internal.within = obj
      joint.derq._internal.within = obj
      joint.w._internal.within = obj
      joint.z._internal.within = obj
      joint.residue_w._internal.within = obj
      joint.residue_f._internal.within = obj
      joint.residue_t._internal.within = obj
      joint.residue_q._internal.within = obj

      return joint
   end
end

#= For index reduction, a vector of the following struct is needed
mutable struct Object3DAccelerationLevel
   z_rel::SVector{3,Float64}  # Relative angular acceleration of Object3D with respect to parent frame, resolved in Object3D
   a0::SVector{3,Float64}     # Absolute acceleration of Object3D origin, resolved in world frame (= der(v0))
   z::SVector{3,Float64}      # Absolute angular acceleration of Object3D, resolved in Object3D
   f::SVector{3,Float64}      # Cut-force resolved in Object3D
   t::SVector{3,Float64}      # Cut-torque resolved in Object3D
   Object3DAccelerationLevel() = new()
end

   A::Vector{Object3DAccelerationLevel}  # A[1]: Derivative of a0/z + forces/torques at Object3D origin
                                      # A[i]: = der(A[i-1])

=#

mutable struct Object3Ddynamics
   v0::SVector{3,Float64}     # Absolute velocity of Object3D origin, resolved in world frame (= der(r_abs))
   w::SVector{3,Float64}      # Absolute angular velocity of Object3D, resolved in Object3D

   a0::SVector{3,Float64}     # Absolute acceleration of Object3D origin, resolved in world frame (= der(v0))
   z::SVector{3,Float64}      # Absolute angular acceleration of Object3D, resolved in Object3D

   f::SVector{3,Float64}      # Cut-force resolved in Object3D
   t::SVector{3,Float64}      # Cut-torque resolved in Object3D

   Object3Ddynamics() = new(ModiaMath.ZeroVector3D, ModiaMath.ZeroVector3D,
                            ModiaMath.ZeroVector3D, ModiaMath.ZeroVector3D,
                            ModiaMath.ZeroVector3D, ModiaMath.ZeroVector3D)
end

"""
    obj1 = Object3D([data]; visualizeFrame=Modia3D.Inherited)
    obj2 = Object3D(parent [, data]; fixed=true, r=zeros(3), R=nothing, q=nothing,
                    v_start=zeros(3), w_start=zeros(3),
                    visualizeFrame=Modia3D.Inherited)

Generate a new Object3D object, that is a coordinate system (= frame) with associated data.
If `parent` is present, the Object3D is defined relatively
to the parent Object3D. If `parent` is not present, the Object3D is either a reference object
(such as the world-object), or the object is connected later with a joint to another
Object3D. If `fixed=true`, the object is rigidly connect to its parent; otherwise
it is moving freely relative to its parent (mathematically described by quaternions).

Note, there are many convenience functions in ModiaMath.Frames to generate a
ModiaMath.RotationMatrix `R` or a ModiaMath.Quaternion `q`.

# Arguments
- `data::Modia3D.AbstractObject3Ddata`: Optional data associated with Object3D.

- `parent::Object3D`: Parent object.

- `fixed::Bool`:
   - `fixed = true`, if the Object3D is fixed relatively to its parent Object3D at position `r,R,q`.
     It is best to provide the rotation information via `R` in this case.
   - `fixed = false`, if Object3D can move freely relatively to its parent Object3D and is
      initially placed at `r,R,q`. The movement is internally described with Quaternion vector `q`.
      Therefore, it is best to provide the rotation information via `q` in this case.

- `r::AbstractVector`: Initial relative position vector from frame of parent object to
   origin of frame object, resolved in parent frame.

- `R::Union{ModiaMath.RotationMatrix,NOTHING}`: Initial rotation matrix defining the rotation
   from frame of parent object to frame of Object3D. If both `R = nothing` and `q = nothing`,
   a null rotation is defined.

- `q::Union{ModiaMath.Quaternion,NOTHING}`: Initial quaternion defining the rotation
   from frame of parent object to frame of Object3D. If both `R = nothing` and `q = nothing`,
   a null rotation is defined.

- `v_start::AbstractVector`: If `fixed=false`, initial velocity of the origin of Object3D with respect to
  parent, resolved in parent frame.

- `w_start::AbstractVector`: If `fixed=false`, initial angular velocity of Object3D with respect to
  parent, resolved in Object3D.

- `visualizeFrame::Union{Bool,Modia3D.Ternary}`: Coordinate system of Object3D is always (= true),
  or never (= false) visualized, or it is visualized if defined in SceneOptions(...)
  (= Modia3D.Inherited).

# Examples
```julia
using Modia3D

# Define assembly
@assembly MyAssembly begin
   world = Object3D()

   # Frame fixed in world
   frame1 = Object3D(world; r=[0.1, 0.2, 0.3])

   # Frame moving relatively to frame1
   r2     = [0.2, 0.2, 0.3]
   frame2 = Object3D(frame1; fixed=false, r=r2)

   # Frame moving relatively to world
   frame3 = Object3D(world; fixed=false, r=-r2)
end
Modia3D.visualizeAssembly!(MyAssembly())
```
"""
mutable struct Object3D <: Modia3D.AbstractAssemblyComponent
   _internal::ModiaMath.ComponentInternal       # Data common to assembly component classes

   # Tree connection structure of Object3D (for efficient processing of Object3Ds)
   parent::Object3D                             # Parent Object3D (if parent===Object3D, no parent is yet defined)
   children::Vector{Object3D}                   # All Object3Ds, where Object3D is the parent

   # Relative position of Object3D with respect to parent Object3D
   joint::Modia3D.AbstractJoint                 # Object that is used to compute r_rel and R_rel
   isDrivenRot::Bool                            # frame or joint is driven from outside
   isDrivenTrans::Bool                          # frame or joint is driven from outside
   isDriven::Bool                               # frame or joint is driven from outside
   r_rel::SVector{3,Float64}                    # Relative position vector from frame of parent Object3D to origin of Object3D frame, resolved in parent frame
                                                # (if parent===Object3D, r_rel=ModiaMath.ZeroVector3D)
   R_rel::SMatrix{3,3,Float64,9}                # Rotation matrix from frame of parent Object3D to Object3D frame.
                                                # (if parent===Object3D, R_rel=ModiaMath.NullRotation)
   # Absolute position of frame
   r_abs::SVector{3,Float64}                    # Absolute position vector from origin of world frame to origin of Object3D frame, resolved in world frame
   R_abs::SMatrix{3,3,Float64,9}                # Absolute rotation matrix from world frame to Object3D frame
   massProperties::Union{Solids.InternalMassProperties,NOTHING}

   # Additional information associated with Object3D
   data::Modia3D.AbstractObject3Ddata                      # Optional data associated with Object3D
   twoObject3Dobject::Vector{Modia3D.AbstractTwoObject3DObject}  # Optional AbstractTwoObject3DObject object associated with Object3D
   hasCutJoint::Bool
   hasForceElement::Bool
   hasChildJoint::Bool
   computeAcceleration::Bool
   visualizeFrame::Modia3D.Ternary                         # = True     : Coordinate system of Object3D is always visualized
                                                           # = False    : Coordinate system of Object3D is never visualized
                                                           # = Inherited: Coordinate system of Object3D is visualized, if SceneOptions(visualizeFrames=true)
   visualizationFrame::Union{Object3D,NOTHING}                # If to be visualized, the Object3D holding the coordinate system.
   contactVisuObj1::Union{Vector{Object3D},NOTHING}
   contactVisuObj2::Union{Vector{Object3D},NOTHING}

   supportVisuObj1A::Union{Vector{Object3D},NOTHING}
   supportVisuObj1B::Union{Vector{Object3D},NOTHING}
   supportVisuObj1C::Union{Vector{Object3D},NOTHING}
   supportVisuObj2A::Union{Vector{Object3D},NOTHING}
   supportVisuObj2B::Union{Vector{Object3D},NOTHING}
   supportVisuObj2C::Union{Vector{Object3D},NOTHING}

   # Data for dynamic simulation
   dynamics::Union{Object3Ddynamics, NOTHING}   # Data for dynamic simulation

   # Object3D constructor without parent
   function Object3D(data::Modia3D.AbstractObject3Ddata = emptyObject3DData;
                     visualizeFrame::Union{Modia3D.Ternary,Bool} = Modia3D.Inherited)::Object3D
      obj                    = new(ModiaMath.ComponentInternal())
      obj.parent             = obj
      obj.children           = Vector{Object3D}[]
      obj.joint              = fixedJoint
      obj.isDrivenRot        = false
      obj.isDrivenTrans      = false
      obj.isDriven           = false
      obj.r_rel              = ModiaMath.ZeroVector3D
      obj.R_rel              = ModiaMath.NullRotation
      obj.r_abs              = ModiaMath.ZeroVector3D
      obj.R_abs              = ModiaMath.NullRotation
      obj.massProperties     = nothing
      obj.data               = data
      obj.twoObject3Dobject  = Modia3D.AbstractTwoObject3DObject[]
      obj.hasCutJoint        = false
      obj.hasForceElement    = false
      obj.hasChildJoint      = false
      obj.computeAcceleration = false
      obj.visualizeFrame     = typeof(visualizeFrame) == Modia3D.Ternary ? visualizeFrame : (visualizeFrame ? Modia3D.True : Modia3D.False)
      obj.visualizationFrame = nothing
      obj.contactVisuObj1    = nothing
      obj.contactVisuObj2    = nothing

      obj.supportVisuObj1A   = nothing
      obj.supportVisuObj1B   = nothing
      obj.supportVisuObj1C   = nothing
      obj.supportVisuObj2A   = nothing
      obj.supportVisuObj2B   = nothing
      obj.supportVisuObj2C   = nothing

      obj.dynamics           = nothing
      return obj
   end

   # Object3D constructor with parent
   function Object3D(parent::Object3D,
                     data::Modia3D.AbstractObject3Ddata = emptyObject3DData;
                     fixed::Bool = true,
                     r::AbstractVector = ModiaMath.ZeroVector3D,
                     R::Union{ModiaMath.RotationMatrix,NOTHING} = nothing,
                     q::Union{ModiaMath.Quaternion,NOTHING} = nothing,
                     v_start::AbstractVector = ModiaMath.ZeroVector3D,
                     w_start::AbstractVector = ModiaMath.ZeroVector3D,
                     visualizeFrame::Union{Modia3D.Ternary,Bool} = Modia3D.Inherited)::Object3D

      if typeof(R) != NOTHING && typeof(q) != NOTHING
         error("Modia3D.Object3D: either R or q must be nothing but both have a value.")
      end
      if typeof(R) != NOTHING
         ModiaMath.assertRotationMatrix(R)
      elseif typeof(q) != NOTHING
         ModiaMath.assertQuaternion(q)
      end

      r_rel = SVector{3,Float64}(r)
      R_rel = typeof(R) != NOTHING ? R : (typeof(q) != NOTHING ? ModiaMath.from_q(q) : ModiaMath.NullRotation)
      r_abs = parent.r_abs + r_rel
      R_abs = R_rel*parent.R_abs

      visualizeFrame2 = typeof(visualizeFrame) == Modia3D.Ternary ? visualizeFrame : (visualizeFrame ? Modia3D.True : Modia3D.False)

      obj = new(ModiaMath.ComponentInternal(), parent, Vector{Object3D}[], fixedJoint,
                false, false, false, r_rel, R_rel, r_abs, R_abs, nothing, data,
                Modia3D.AbstractTwoObject3DObject[], false, false, false, false, visualizeFrame2, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing)

      if !fixed
         obj.joint = FreeMotion(obj; r_start = r_rel,
                                q_start = typeof(R) != NOTHING ? ModiaMath.from_R(R) :
                                          typeof(q) != NOTHING ? q : ModiaMath.NullQuaternion,
                                v_start = SVector{3,Float64}(v_start),
                                w_start = SVector{3,Float64}(w_start) )
      end

      push!(parent.children, obj)
      return obj
   end

   # Constructor used only for internal purposes (not to be directly used by the user)
   Object3D(_internal::ModiaMath.ComponentInternal,
            parent::Object3D, children::Vector{Object3D},
            joint::Modia3D.AbstractJoint, r_rel::SVector{3,Float64},
            R_rel::SMatrix{3,3,Float64,9}, r_abs::SVector{3,Float64},
            R_abs::SMatrix{3,3,Float64,9}, data::Modia3D.AbstractObject3Ddata,
            twoObject3Dobject::Vector{Modia3D.AbstractTwoObject3DObject},
            visualizeFrame::Modia3D.Ternary) =
      new(_internal, parent, children, joint, false, false, false, r_rel,
          R_rel, r_abs, R_abs, nothing, data, twoObject3Dobject, false, false, false, false, visualizeFrame, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing)
end

# Object3D constructor without parent and with convex decomposition of the solid
function Object3D(data::Solids.SolidWithConvexDecomposition;
                  visualizeFrame::Union{Modia3D.Ternary,Bool} = Modia3D.Inherited)::Object3D
   obj      = Object3D(data.solid; visualizeFrame=visualizeFrame)
   elements = data.convexDecomposition
   names    = data.names
   for i in eachindex(elements)
      obj_element = Object3D(obj, elements[i]; visualizeFrame=Modia3D.False)
      obj_element._internal.within = obj
      obj_element._internal.name   = names[i]
   end
   return obj
end

# Object3D constructor with parent and with convex decomposition of the solid
function Object3D(parent::Object3D,
                  data::Solids.SolidWithConvexDecomposition;
                  fixed::Bool = true,
                  r::AbstractVector = ModiaMath.ZeroVector3D,
                  R::Union{ModiaMath.RotationMatrix,NOTHING} = nothing,
                  q::Union{ModiaMath.Quaternion,NOTHING} = nothing,
                  v_start::AbstractVector = ModiaMath.ZeroVector3D,
                  w_start::AbstractVector = ModiaMath.ZeroVector3D,
                  visualizeFrame::Union{Modia3D.Ternary,Bool} = Modia3D.Inherited)::Object3D
   obj      = Object3D(parent, data.solid; fixed=fixed, r=r, R=R, q=q, v_start=v_start, w_start=w_start, visualizeFrame=visualizeFrame)
   elements = data.convexDecomposition
   names    = data.names
   for i in eachindex(elements)
      obj_element = Object3D(obj, elements[i]; visualizeFrame=Modia3D.False)
      obj_element._internal.within = obj
      obj_element._internal.name   = Symbol(names[i])
   end
   return obj;
end

function driveJoint!(frame::Object3D)::NOTHING
   frame.isDriven = true
   return nothing
end

# Make a direct copy of the first level of frame (only used internally in Modia3D)
copyObject3D(obj::Object3D, data::Modia3D.AbstractObject3Ddata) = Object3D(obj._internal, obj.parent, obj.children, obj.joint,
                                                                           obj.r_rel, obj.R_rel, obj.r_abs, obj.R_abs, data,
                                                                           obj.twoObject3Dobject, obj.visualizeFrame)

# Inquire properties of a Object3D
hasParent(            obj::Object3D) = !(obj.parent === obj)
hasNoParent(          obj::Object3D) =   obj.parent === obj
hasChildren(          obj::Object3D) = length(obj.children) > 0
hasNoChildren(        obj::Object3D) = length(obj.children) == 0
isWorld(              obj::Object3D) = hasNoParent(obj) && ModiaMath.isInComponent(obj) && typeof(obj._internal.within._internal.scene) == Scene && obj._internal.within._internal.scene.initAnalysis
isNotWorld(           obj::Object3D) = !(isWorld(obj))
isFixed(              obj::Object3D) = typeof(obj.joint) == FixedJoint || typeof(obj.joint) == TreeJointFixed
isNotFixed(           obj::Object3D) = typeof(obj.joint) != FixedJoint && typeof(obj.joint) != TreeJointFixed
isFree(               obj::Object3D) = typeof(obj.joint) == FreeMotion
isNotFree(            obj::Object3D) = typeof(obj.joint) != FreeMotion
hasJoint(             obj::Object3D) = isNotFixed(obj) && isNotFree(obj)
hasNoJoint(           obj::Object3D) = isFixed(obj) || isFree(obj)
hasData(              obj::Object3D) = typeof(obj.data) != EmptyObject3Ddata
isCoordinateSystem(   obj::Object3D) = typeof(obj.data) == Graphics.CoordinateSystem
isNotCoordinateSystem(obj::Object3D) = typeof(obj.data) != Graphics.CoordinateSystem
hasCutJoint(          obj::Object3D) = obj.hasCutJoint
hasForceElement(      obj::Object3D) = obj.hasForceElement
hasChildJoint(        obj::Object3D) = obj.hasChildJoint
needsAcceleration(    obj::Object3D) = obj.computeAcceleration
objectHasMass(        obj::Object3D) = typeof(obj.massProperties)!= NOTHING

# The following functions need to be defined for every type derived from AbstractObject3Ddata,
# provided the function should return true

isVisible( data::Modia3D.AbstractObject3Ddata, renderer::Modia3D.AbstractRenderer) = false
canCollide(data::Modia3D.AbstractObject3Ddata) = false
dataHasMass(   data::Modia3D.AbstractObject3Ddata) = false
dataHasMass(   data::Solids.Solid)                 = typeof(data.massProperties)!= NOTHING

# Inquire properties of an Object3D that depend on the type of frame.data
"""
    isVisible(obj::Object3D, renderer::Modia3D.AbstractRenderer)

returns true, if `obj` is visualized with renderer `renderer`.
"""
isVisible( obj::Object3D, renderer::Modia3D.AbstractRenderer) = isVisible(obj.data, renderer)

function canCollide(obj::Object3D)::Bool
  if typeof(obj.data) == Modia3D.Solids.Solid
    if typeof(obj.data.geo) != Nothing && typeof(obj.data.contactMaterial) != Nothing
        if typeof(obj.data.contactMaterial) == String
            return obj.data.contactMaterial != ""
        else
            return true
        end
    else
        return false
    end
  else
    return false
  end
end
dataHasMass(obj::Object3D) = dataHasMass(obj.data)

"""    rootObject3D(frame) - returns the root frame of all parents of frame"""
function rootObject3D(frame::Object3D)::Object3D
   frame1 = frame
   while hasParent(frame1)
      frame1 = frame1.parent
   end
   return frame1
end

"""    removeChild!(frame, child) - Remove child from frame.children"""
function removeChild!(frame::Object3D, child::Object3D)::NOTHING
   children = frame.children
   for i in eachindex(children)
      if children[i] === child
         deleteat!(children,i)
         return nothing
      end
   end
   error("\nError from Modia3D.Composition.removeChild!(", frame, ", ", child, ")\n",
         child, " is not a child of frame ", frame)
end

# Print Object3D
#   dataName = Basics.trailingPartOfTypeAsString(frame.data)
function Base.show(io::IO, frame::Object3D)
   commaNeeded = true
   if frame.parent === frame || hasJoint(frame)
      commaNeeded = false
      print(io,"Object3D(")
   else
      print(io,"Object3D(", ModiaMath.fullName(frame.parent))
   end

   if typeof(frame.data) != EmptyObject3Ddata
      if commaNeeded
         print(io,", ")
      end
      typeName = string( typeof(frame.data) )
      s = findfirst("Modia3D", typeName) # search(typeName, "Modia3D")
      if s.start == 1
         typeName = "Modia3D." * Basics.trailingPartOfName( typeName )
      end
      print(io, typeName, "(...)")
   end

   if frame.r_rel != ModiaMath.ZeroVector3D
      print(io, "; r_rel = ", frame.r_rel, ")")
   else
      print(io, ")")
   end
   print(io, " # ", ModiaMath.instanceName(frame))
end