text:
# ==============================================================================
# Max-Plus Algebra toolbox for Julia >= 1.0.3
# In the way of the ScicosLab max-plus toolbox.
# ==============================================================================

module MaxPlus

using
    LinearAlgebra, SparseArrays, Printf

# ==============================================================================
# Max-Plus core

export
    MP, SpaMP, SpvMP, ArrMP, VecMP,
    mpzero, mpone, mp0, mp1, mptop,
    mpI, mpeye, mpzeros, mpones,
    mpsparse, full, dense, array,
    plustimes, minplus,
    mptrace, mpnorm, mpstar,
    mp_change_display, LaTeX

# ==============================================================================
# Max-Plus Linear system

export
    MPSysLin,
    mpsyslin, mpsimul, mpexplicit

# ==============================================================================
# Max-Plus flowhop

export
    mpgraph,
    flowshop, LaTeXG, flowshop_graph, flowshop_simu

# ==============================================================================

"""
    MP{T}

Max-Plus immutable structure. Promote a number of type T (ideally a Float64) to
a max-plus number.

# Examples
```julia-repl
julia> a=MP(3.0); typeof(a)
MP{Float64}

julia> a=MP(3); typeof(a)
MP{Int64}
```
"""
struct MP{T <: Real} <: Real
    λ::T
end

# ==============================================================================
# Type alias but with shorter number of characters

const Sparse{T,U} = SparseMatrixCSC{T,U}
const SpaMP{T,U}  = SparseMatrixCSC{MP{T},U}
const SpvMP{T,U}  = SparseVector{MP{T},U}
const ArrMP{T,N}  = Array{MP{T},N}
const VecMP{T,N}  = Vector{MP{T}}

# ==============================================================================
# Copy constructor

MP(x::MP) = MP(x.λ)

# ==============================================================================

"""
    zero(::MP{T})

Create the max-plus constant zero equal to -Inf.
This value is the neutral for the ⨁ operator.
"""
Base.zero(::Type{MP{T}}) where T = MP(typemin(T))
Base.zero(x::MP{T})      where T = zero(typeof(x))

"""
    mpzero(::Type{T})

Create the max-plus constant zero equal to -Inf.
This value is the neutral for the ⨁ operator.
"""
mpzero(::Type{T})        where T = MP(typemin(T))

# ==============================================================================

"""
    one(::MP{T})

Create the max-plus constant one equal to 0.
This value is the neutral for operators ⨁ and ⨂.
"""
Base.one(::Type{MP{T}})  where T = MP(zero(T))
Base.one(x::MP{T})       where T = one(typeof(x))

"""
    one(::Type{T})

Create the max-plus constant one equal to 0.
This value is the neutral for operators ⨁ and ⨂.
"""
mpone(::Type{T})         where T = MP(zero(T))

# ==============================================================================

"""
    mp0

Max-plus constant zero equal to -Inf.
This value is the neutral for the ⨁ operator.

Equivalent to ScicosLab code: `%0 = #(-%inf) = .`

# Examples
```julia-repl
julia> mp0 * 5
-Inf

julia> mp0 + 5
5.0
```
"""
const global mp0 = mpzero(Float64)

# ==============================================================================

"""
    mp1

Max-plus constant one equal to 0.
This value is the neutral for operators ⨁ and ⨂.

Equivalent to ScicosLab code: `%1 = #(1) = 0`

# Examples
```julia-repl
julia> mp1 * 5
5.0

julia> mp1 + 5
5.0
```
"""
const global mp1 = mpone(Float64)

# ==============================================================================

"""
    mptop

Min-plus constant "top" equal to +Inf.

Equivalent to ScicosLab code: `%top = #(%inf)`
"""
const global mptop = MP{Float64}(Inf)

# ==============================================================================
# Conversions

Base.promote_rule(::Type{MP{T}}, ::Type{U}) where {T, U} = MP{T}
Base.convert(::MP{T}, x::Number)            where T = MP{T}(x)

# ==============================================================================

"""
    MP(x::UnitRange)

Create a max-plus dense vector.

# Examples
```julia-repl
julia> MP(1:3)
3-element Array{MP{Int64},1}:
 1
 2
 3
```
"""
MP(x::UnitRange{T}) where T = MP(Vector{T}(x))

# ==============================================================================

"""
    MP(A::Array)

Convert a dense array to a max-plus dense array.

# Examples
```julia-repl
julia> MP([1.0 -Inf; 0.0 4])
2×2 Array{MP{Float64},2}:
 1.0  -Inf
 0.0  4.0
```
"""
MP(A::Array) = map(MP, A)

# ==============================================================================

"""
    MP(A::SparseArray)

Convert a sparse array to a max-plus sparse array.
By default values are max-plus zeros values are preserved else if
parameter `preserve` is set then

# Examples
```julia-repl
julia> MP(sparse([1, 2, 3], [1, 2, 3], [-Inf, 2, 0]))
  [1, 1]  =  -Inf
  [2, 2]  =  2.0
  [3, 3]  =  0.0

julia> MP(sparse([1, 2, 3], [1, 2, 3], [-Inf, 2, 0]), preserve=false)
  [2, 2]  =  2.0
```
"""
function MP(S::Sparse{T,U}; preserve=true) where {T, U}
    if preserve
        convert(SpaMP{T,U}, S)
    else
        M = spzeros(MP{T}, size(S,1), size(S,2))
        for c = 1:size(S, 2), r = nzrange(S, c)
            if (S[r,c] != zero(T)) && (MP(S[r,c]) != mpzero(T))
                M[r,c] = convert(MP{T}, S[r,c])
            end
        end
        M
    end
end

# ==============================================================================

"""
    MP(A::SparseVector)

Convert a sparse vector to a max-plus sparse vector.
By default values are max-plus zeros values are preserved else if
parameter `preserve` is set then

# Examples
```julia-repl
julia> MP(sparse([1.0, 0.0, 1.0]))
3-element SparseVector{MP{Float64},Int64} with 2 stored entries:
  [1]  =  1.0
  [3]  =  1.0
```
"""
function MP(V::SparseVector{T,U}; preserve=true) where {T, U}
    if preserve
        convert(SparseVector{MP{T},U}, V)
    else
        S = spzeros(MP{T}, size(V,1), size(V,2))
        for c = V.n
            if (S[c] != zero(T)) && (MP(S[c]) != mpzero(T))
                M[c] = convert(MP{T}, S[c])
            end
        end
        S
    end
end

# ==============================================================================

global mpstyle = 1 # Style ScicosLab

"""
    mp_change_display(style::Int)

Change the style of display for the Base.show(io::IO, x::MP) function.
See help of show(io::IO,x::MP) for examples.
"""
mp_change_display(style::Int) = (global mpstyle = min(max(style, 0), 4))

"""
    show(io::IO,x::MP)

Display a max-plus number depending on the currently set style:
- -Inf are displayed either with 'ε' (style 2 or 3) or '.' symbols (style 1).
- 0 are displayed either with 'e' (style 3) or '0' symbols (style 1 or 2).
- else -Inf and 0 are displayed in Julia default sytle (style 0).
By default the ScicosLab style is used (style 1).

# Examples
```julia-repl
julia> mp_change_display(0); mpeye(Float64, 2,2)
2×2 Array{MP{Float64},2}:
  0.0  -Inf
 -Inf   0.0

julia> mp_change_display(0); mpeye(Int64, 2,2)
2×2 Array{MP{Int64},2}:
                    0  -9223372036854775808
 -9223372036854775808                     0

julia> mp_change_display(1); mpeye(Float64, 2,2) # or mpeye(Int64, 2,2)
2×2 Array{MP{Float64},2}:
 0  .
 .  0

julia> mp_change_display(2); mpeye(Float64, 2,2) # or mpeye(Int64, 2,2)
2×2 Array{MP{Float64},2}:
 e  .
 .  e

julia> mp_change_display(3); mpeye(Float64, 2,2) # or mpeye(Int64, 2,2)
2×2 Array{MP{Float64},2}:
 0  ε
 ε  0

julia> mp_change_display(4); mpeye(Float64, 2,2) # or mpeye(Int64, 2,2)
2×2 Array{MP{Float64},2}:
 e  ε
 ε  e
```
"""
function Base.show(io::IO, x::MP{T}) where T
    if mpstyle == 0
        show(io, x.λ)
    elseif x == mpzero(T)
        (mpstyle == 1 || mpstyle == 2) ? (@printf io ".") : (@printf io "ε")
    elseif x == mpone(T)
        (mpstyle == 1 || mpstyle == 3) ? (@printf io "0") : (@printf io "e")
    elseif x.λ == trunc(x.λ)
        (@printf io "%d" x.λ)
    else
        show(io, x.λ)
    end
end

# ==============================================================================
# FIXME: not a column-major traversal
"""
    LaTeX(io::IO, A::Array{MP})

Convert a max-plus dense matrix to a LaTeX formula. The display depends on
mp_change_display(style).

# Examples
```julia-repl
julia> LaTeX(stdout, MP([4 3; 7 -Inf]))
\\left[
\\begin{array}{*{20}c}
4 & 3 \\\\
7 & . \\\\
\\end{array}
\\right]
```
"""
function LaTeX(io::IO, A::ArrMP{T}) where T
    (@printf io "\\left[\n\\begin{array}{*{20}c}\n")
    for i in 1:size(A,1)
        for j in 1:size(A,2)
            if A[i,j] == mpzero(T)
                if mpstyle == 0
                    (@printf io "-\\infty")
                elseif mpstyle == 3 || mpstyle == 4
                    (@printf io "\\varepsilon")
                else
                    (@printf io ".")
                end
            elseif A[i,j] == mpone(T)
                if mpstyle == 2 || mpstyle == 4
                    (@printf io "e")
                else
                    (@printf io "%d" A[i,j].λ)
                end
            elseif A[i,j].λ == trunc(A[i,j].λ)
                (@printf io "%d" A[i,j].λ)
            else
                show(io, A[i,j].λ)
            end
            if j < size(A, 2)
                (@printf io " & ")
            end
        end
        (@printf io " \\\\\n")
    end
    (@printf io "\\end{array}\n\\right]\n")
end

LaTeX(io::IO, S::SpaMP{T}) where {T,U} = LaTeX(io, full(S))

# ==============================================================================

"""
    +(x::MP, y::MP)

Max operator ⨁. Return the maximum of `x` and `y`.

# Examples
```julia-repl
julia> MP(1.0) + MP(3.0)
MP{Float64}(3.0)
```
"""
Base.:(+)(x::MP,   y::MP)   = MP(max(x.λ, y.λ))
Base.:(+)(x::MP,   y::Real) = MP(max(x.λ, y))
Base.:(+)(x::Real, y::MP)   = MP(max(x,   y.λ))

# ==============================================================================

"""
    *(x::MP, y::MP)

Addition operator ⨂. Return the summation of `x` and `y`.

# Examples
```julia-repl
julia> MP(1.0) * MP(3.0)
MP{Float64}(4.0)
```
"""
Base.:(*)(x::MP,   y::MP)   = MP(x.λ + y.λ)
Base.:(*)(x::MP,   y::Real) = MP(x.λ + y)
Base.:(*)(x::Real, y::MP)   = MP(x   + y.λ)

# ==============================================================================
# These functions fix Julia bugs

@inline Base.literal_pow(::typeof(^), x::MP, ::Val{0}) = one(x)
@inline Base.literal_pow(::typeof(^), x::MP, ::Val{p}) where {p} = MP(x.λ * p)
@inline Base.:(^)(x::MP, y::Number) = MP(x.λ * y)
@inline Base.literal_pow(::typeof(^), A::ArrMP{T}, ::Val{p}) where {T, p} = A^p
@inline Base.abs2(x::MP) = x.λ + x.λ

# ==============================================================================

"""
    /(x::MP, y::MP)

Divisor operator. Return the difference between `x` and `y`.

# Examples
```julia-repl
julia> MP(1.0) / MP(2.0)
MP{Float64}(-1.0)
```
"""
Base.:(/)(x::MP,   y::MP)   = MP(x.λ - y.λ)
Base.:(/)(x::MP,   y::Real) = MP(x.λ - b)
Base.:(/)(x::Real, y::MP)   = MP(a   - y.λ)

# ==============================================================================

"""
    -(x::MP, y::MP)

Minus operator (not used for max-plus but for min-plus).
Return the difference between `x` and `y`.

# Examples
```julia-repl
julia> MP(1.0) - MP(2.0)
MP{Float64}(-1.0)
```
"""
Base.:(-)(x::MP,   y::MP)   = MP(x.λ - y.λ)
Base.:(-)(x::MP,   y::Real) = MP(x.λ - y)
Base.:(-)(x::Real, y::MP)   = MP(x   - y.λ)

# ==============================================================================

"""
    min(x::MP, y::MP)

Return the minimun of `x` and `y`.

# Examples
```julia-repl
julia> min(MP(1), -3)
MP{Int64}(-3)

julia> min(MP([10 1; 10 1]), MP([4 5; 6 5]))
2×2 Array{MP{Int64},2}:
 4  1
 6  1
```
"""
Base.:min(x::MP,   y::MP)   = MP(min(x.λ, y.λ))
Base.:min(x::MP,   y::Real) = MP(min(x.λ, y))
Base.:min(x::Real, y::MP)   = MP(min(x,   y.λ))
Base.:min(A::ArrMP, B::ArrMP) = map(Base.:min, A, B)
Base.:min(A::SpaMP, B::SpaMP) = map(Base.:min, A, B)

# ==============================================================================

Base.:(==)(x::MP,   y::MP)   = (x.λ == y.λ)
Base.:(==)(x::MP,   y::Real) = (x.λ == y)
Base.:(==)(x::Real, y::MP)   = (x   == y.λ)

Base.:(<=)(x::MP,   y::MP)   = (x.λ <= y.λ)
Base.:(<=)(x::MP,   y::Real) = (x.λ <= y)
Base.:(<=)(x::Real, y::MP)   = (x   <= y.λ)

Base.:(<)(x::MP,   y::MP)   = (x.λ < y.λ)
Base.:(<)(x::MP,   y::Real) = (x.λ < y)
Base.:(<)(x::Real, y::MP)   = (x   < y.λ)

Base.isless(x::MP,   y::MP)   = x.λ < y.λ
Base.isless(x::MP,   y::Real) = x.λ < y
Base.isless(x::Real, y::MP)   = x < y.λ

# ==============================================================================

"""
    mpeye(::Type{T}, n::Int64)

Construct a max-plus identity n-by-n matrix.

# Examples
```julia-repl
julia> mpeye(Float64, 2)
2×2 Array{MP{Float64},2}:
  0.0  -Inf
 -Inf   0.0
```
"""
mpeye(::Type{T}, n::Int64) where T = Matrix{MP{T}}(mpI, n, n)

# ==============================================================================

"""
    mpeye(::Type{T}, m::Int64, n::Int64)

Construct a max-plus identity m-by-n matrix.

# Examples
```julia-repl
julia> mpeye(Float64, 2, 2)
2×2 Array{MP{Float64},2}:
  0.0  -Inf
 -Inf   0.0
```
"""
mpeye(::Type{T}, m::Int64, n::Int64) where T = Matrix{MP{T}}(mpI, m, n)

# ==============================================================================

"""
    mpzeros(::Type{T}, n::Int64)

Construct a max-plus zero n-by-m sparse matrix.

# Examples
```julia-repl
julia> mpzeros(Float64, 2)
2-element SparseVector{MP{Float64},Int64} with 0 stored entries
```
"""
mpzeros(::Type{T}, n::Int64) where T = spzeros(MP{T}, n)

# ==============================================================================

"""
    mpzeros(::Type{T}, n::Int64)

Construct a sparse max-plus zero m-by-n matrix.

# Examples
```julia-repl
julia> mpzeros(Float64, 2,5)
2×5 SparseMatrixCSC{MP{Float64},Int64} with 0 stored entries
```
"""
mpzeros(::Type{T}, m::Int64, n::Int64) where T = spzeros(MP{T}, m, n)

"""
    mpfzeros(::Type{T}, n::Int64)

Construct a dense max-plus zero m-by-n matrix.

# Examples
```julia-repl
julia> mpfzeros(Float64, 2,2)
2×2 Array{MP{Float64},2}:
 -Inf  -Inf
 -Inf  -Inf
```
"""
mpfzeros(::Type{T}, m::Int64, n::Int64) where T = full(mpzeros(T, m, n))

# ==============================================================================

"""
    mpzeros(x::UnitRange{T})

Construct a max-plus zero vector.
# Examples
```julia-repl
julia> mpzeros(1:3)
3-element Array{MP{Int64},1}:
 -Inf
 -Inf
 -Inf
```
"""
mpzeros(x::UnitRange{T}) where {T} = fill!(MP(similar(x)), mpzero(T))

# ==============================================================================

"""
    mpones(::Type{T}, n::Int64)

Construct a max-plus one n-by-n matrix.

# Examples
```julia-repl
julia> mpones(Float64, 2)
2-element Array{MP{Float64},1}:
 0.0
 0.0
```
"""
mpones(::Type{T}, n::Int64) where T = ones(MP{T}, n)

# ==============================================================================

"""
    mpones(::Type{T}, m::Int64, n::Int64)

Construct a max-plus one m-by-n matrix.

# Examples
```julia-repl
julia> mpones(Float64, 2,2)
2×2 Array{MP{Float64},2}:
 0.0  0.0
 0.0  0.0
```
"""
mpones(::Type{T}, m::Int64, n::Int64) where T = ones(MP{T}, m, n)

# ==============================================================================

"""
    mpones(x::UnitRange{T})

Construct a max-plus one vector.
# Examples
```julia-repl
julia> mpones(1:3)
3-element Array{MP{Int64},1}:
 0
 0
 0
```
"""
mpones(x::UnitRange{T}) where {T} = fill!(MP(similar(x)), mpone(T))

# ==============================================================================

"""
    mpsparse(A::Array{T}; preserve=false)

Transform a dense matrix to a max-plus sparse matrix.

# Examples
```julia-repl
julia> mpsparse([4 0; 7 -Inf])
2×2 SparseMatrixCSC{MP{Float64},Int64} with 3 stored entries:
  [1, 1]  =  4.0
  [2, 1]  =  7.0
  [1, 2]  =  0.0
```
"""
mpsparse(A::Array{T}) where T = mpsparse(MP(A))
# ; preserve=false) where T = mpsparse(MP(A), preserve)

"""
    mpsparse(A::ArrMP{T}; preserve=false)

Transform a dense matrix to a max-plus sparse matrix. By default
Remove -Inf values. To keep them set param preserve to true.

# Examples
```julia-repl
julia> mpsparse([4 0; 7 -Inf])
2×2 SparseMatrixCSC{MP{Float64},Int64} with 3 stored entries:
  [1, 1]  =  4.0
  [2, 1]  =  7.0
  [1, 2]  =  0.0

mpsparse([4 0; 7 -Inf], preserve=true)
2×2 SparseMatrixCSC{MP{Float64},Int64} with 3 stored entries:
  [1, 1]  =  4.0
  [2, 1]  =  7.0
  [1, 2]  =  0.0
  [2, 2]  = -Inf
```
"""
function mpsparse(M::ArrMP{T}) where T #; preserve=false) where T
    #if preserve
    #    sparse(M)
    #else
        R = Array{Int64, 1}[];
        C = Array{Int64, 1}[];
        V = Array{MP{T}, 1}[];

        for j = 1:size(M, 2), i=1:size(M, 1)
            M[i,j] != zero(MP{T}) && (R = [R; i]; C = [C; j]; V = [V; M[i,j]])
        end

        sparse(R, C, V)
    #end
end

# ==============================================================================

"""
    array(::ArrMP{T})

Convert a max-plus array to a standard type array.

# Examples
```julia-repl
julia> array([MP(1.0) 2.0; 3.0 4.0])
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0
```
"""
array(A::ArrMP{T}) where T = map(x -> x.λ, A)

# ==============================================================================

"""
    array(::SpaMP{T}})

Convert a sparse max-plus array to a dense non max-plus array.

# Examples
```julia-repl
julia> array(mpzeros(Float64, 2,2))
2×2 SparseMatrixCSC{Float64,Int64} with 4 stored entries:
  [1, 1]  =  -Inf
  [2, 1]  =  -Inf
  [1, 2]  =  -Inf
  [2, 2]  =  -Inf
```
"""
array(S::SpaMP{T,U}) where {T, U} = map(x -> x.λ, S)

# ==============================================================================

SparseArrays.sparse(S::SpaMP{T,U}) where {T, U} = map(x -> x.λ, S)

# ==============================================================================

"""
    full(::SpaMP{T}})

Convert a sparse max-plus array to a dense max-plus array.
Alternative function name: dense.

# Examples
```julia-repl
julia> full(mpzeros(Float64, 2,5))
2×5 Array{MP{Float64},2}:
 -Inf  -Inf  -Inf  -Inf  -Inf
 -Inf  -Inf  -Inf  -Inf  -Inf
```
"""
full(S::SpaMP{T,U}) where {T, U} = Matrix(S)

# ==============================================================================

"""
    dense(::SpaMP{T}})

Convert a sparse max-plus array to a dense max-plus array.
Alternative function name: full.

# Examples
```julia-repl
julia> dense(mpzeros(Float64, 2,5))
2×5 Array{MP{Float64},2}:
 -Inf  -Inf  -Inf  -Inf  -Inf
 -Inf  -Inf  -Inf  -Inf  -Inf
```
"""
dense(S::SpaMP{T,U}) where {T, U} = Matrix(S)

# ==============================================================================

"""
    plustimes(x::MP{T})

Convert max-plus number to the standard type.
Note: the function name comes from ScicosLab.

# Examples
```julia-repl
julia> typeof(plustimes(MP(1)))
Int64

julia> plustimes([MP(1.0) 2.0; 3.0 4.0])
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0
```
"""
plustimes(n::MP{T})      where T      = n.λ
plustimes(A::ArrMP{T})   where T      = array(A)
plustimes(S::SpaMP{T,U}) where {T, U} = array(S)

# ==============================================================================
# Conversion to min-plus algebra. Convert +Inf and -Inf

"""
    minplus(A::Array)

Conversion to min-plus algebra. Convert +Inf and -Inf

# Examples
```julia-repl
julia> A = MP([0 3 Inf 1; 1 2 2 -Inf; -Inf Inf 1 0])
minplus(A);
```
"""
minplus(n::MP{T})      where T      = map(x -> (x.λ == mp0) ? mptop : ((x.λ == mptop) ? mp0 : x), n)
minplus(A::ArrMP{T})   where T      = map(x -> minplus(x), A)
minplus(S::SpaMP{T,U}) where {T, U} = map(x -> minplus(x), S)

# ==============================================================================

"""
    mptrace(A::Array)

Compute the trace of the matrix (summation of diagonal elements).

# Examples
```julia-repl
julia> mptrace([MP(1) 2; 3 4])
4
```
"""
mptrace(A::ArrMP{T})    where T = isempty(A) ? mp0 : sum(diag(A))
mptrace(A::Array{T})    where T = isempty(A) ? mp0 : sum(MP(diag(A)))
mptrace(S::SpaMP{T,U})  where {T, U} = isempty(S) ? mp0 : sum(diag(S))
mptrace(S::Sparse{T,U}) where {T, U} = sum(MP(diag(S)))

# ==============================================================================

"""
    mpnorm(A)

Compute the norm of the full or sparce matrix A.
Return the largest entry minus smallest entry of A.

# Examples
```
julia-repl
julia> A = MP([1 20 2;30 400 4;4 50 10]);
S = mpsparse(A);
[mpnorm(A) mpnorm(S)]
1×2 Array{MP{Int64},2}:
 399  399
```
"""
mpnorm(A::Array{T}) where T = A[argmax(A)].λ - A[argmin(A)].λ
mpnorm(S::SpaMP{T}) where T = S[argmax(S)].λ - S[argmin(S)].λ

# ==============================================================================
function hstar(B::Array{T}) where T
    n = size(B, 1)
    if n != size(B, 2) error("Matrix shall be squared") end
    C = B
    for k in 1:n
        t = (C[k,k] <= zero(T)) ? zero(T) : typemax(T);
        for j in 1:n, i in 1:n
            C[i,j] = max(C[i,j], C[i,k] + C[k,j] + t)
        end
    end

    for k in 1:n
        C[k,k] = max(C[k,k], zero(T));
    end
    C
end

"""
    mpstar(A::Array{T})

TODO
"""
mpstar(x::MP{T})    where T = mpstar(plustimes(x))
mpstar(x::T)        where T = MP(hstar([x]))[1,1]
mpstar(A::Array{T}) where T = MP(hstar(A))
mpstar(A::ArrMP{T}) where T = MP(hstar(plustimes(A)))
#mpstar(S::Sparse{T,U}) where {T, U} = MP(hstar(S))
#mpstar(S::SpaMP{T,U}) where {T, U} = MP(hstar(plustimes(S)))

# ==============================================================================
# TODO insertion out of bounds => Sparse:  d=MP([1.0 2; 3 4]); d[5,5] = MP(6.0)

# ==============================================================================
include("juliabugs.jl")
include("howard.jl")
include("syslin.jl")
include("flowshop.jl")

end # MaxPlus module


text:
global defaults = Dict(:maxlen_widget_array => 9,
                       :minlen_plottable_array => 10,
                       :missing_label => "unknown",
                       :plot_range => 1:10000000
                  )

"""
    set_defaults(args...)

Set global defaults to control behaviour.

# Examples

```julia
set_defaults(:maxlen_slider_array => 10,
             :minlen_plot_array => 10,
             :missing_label => "List"
            )
```
"""
function set_defaults(args...)
    global defaults
    for (key, val) in args
        defaults[key] = val
    end
end

"""
    get_box(::T)

Returns the vbox or hbox function depending on input type.
Add methods to control the layout. Associtave structures are displayed
vertically by default, otherwise horizontally.

# Examples

Use a vertical layout for tuples:

```julia
using AutoInteract.get_box
AutoInteract.get_box(::Tuple) = vbox
```
"""
get_box(::Associative) = vbox
get_box(::Any) = hbox

"""
    get_range(x::Number, label)

Set the range of a slider. Add methods to customise.
The widget label provided if you need to modify specific widgets.

# Examples

Make all sliders range from 0.0 to 100.0, with a step of 2.0 for Float64:

```julia
using AutoInteract.get_range
AutoInteract.get_range(x::Float64, label) = 0.0:2.0:100.0
```
"""
function get_range(x::AbstractFloat, label)
    low = 0.0
    if x == 0.0
        high = 1.0
    else
        high = x * 2
    end
    step = (high - low) / 100
    return low:step:high
end
function get_range(x::Integer, label)
    step = 1
    if x == 0
        range = 0:step:100
    elseif x < 0
        range = 2x:step:0
    else
        range = 0:step:2x
    end
    return range
end

@require Unitful begin
function get_range(x::Unitful.Quantity, label)
    u = Unitful.unit(x)
    low = 0.0u
    if x == 0.0u
        high = 1.0u
    else
        high = x * 2
    end
    step = (high - low) / 100
    return low:step:high
end
end

"""
    plotit!(p, parent, data, label)

Common plotting function. Can extend specific methods for a parent type,
or customise for a particular label.

p is a Plots object to plot to.
"""
plotit!(p, parent, data, label) =
    plot!(p, data[inbounds_range(data)], label=label, layout=deflayout)

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct StorageAccountMicrosoftEndpoints <: SwaggerModel
    blob::Any # spec type: Union{ Nothing, String } # spec name: blob
    queue::Any # spec type: Union{ Nothing, String } # spec name: queue
    table::Any # spec type: Union{ Nothing, String } # spec name: table
    file::Any # spec type: Union{ Nothing, String } # spec name: file
    web::Any # spec type: Union{ Nothing, String } # spec name: web
    dfs::Any # spec type: Union{ Nothing, String } # spec name: dfs

    function StorageAccountMicrosoftEndpoints(;blob=nothing, queue=nothing, table=nothing, file=nothing, web=nothing, dfs=nothing)
        o = new()
        validate_property(StorageAccountMicrosoftEndpoints, Symbol("blob"), blob)
        setfield!(o, Symbol("blob"), blob)
        validate_property(StorageAccountMicrosoftEndpoints, Symbol("queue"), queue)
        setfield!(o, Symbol("queue"), queue)
        validate_property(StorageAccountMicrosoftEndpoints, Symbol("table"), table)
        setfield!(o, Symbol("table"), table)
        validate_property(StorageAccountMicrosoftEndpoints, Symbol("file"), file)
        setfield!(o, Symbol("file"), file)
        validate_property(StorageAccountMicrosoftEndpoints, Symbol("web"), web)
        setfield!(o, Symbol("web"), web)
        validate_property(StorageAccountMicrosoftEndpoints, Symbol("dfs"), dfs)
        setfield!(o, Symbol("dfs"), dfs)
        o
    end
end # type StorageAccountMicrosoftEndpoints

const _property_map_StorageAccountMicrosoftEndpoints = Dict{Symbol,Symbol}(Symbol("blob")=>Symbol("blob"), Symbol("queue")=>Symbol("queue"), Symbol("table")=>Symbol("table"), Symbol("file")=>Symbol("file"), Symbol("web")=>Symbol("web"), Symbol("dfs")=>Symbol("dfs"))
const _property_types_StorageAccountMicrosoftEndpoints = Dict{Symbol,String}(Symbol("blob")=>"String", Symbol("queue")=>"String", Symbol("table")=>"String", Symbol("file")=>"String", Symbol("web")=>"String", Symbol("dfs")=>"String")
Base.propertynames(::Type{ StorageAccountMicrosoftEndpoints }) = collect(keys(_property_map_StorageAccountMicrosoftEndpoints))
Swagger.property_type(::Type{ StorageAccountMicrosoftEndpoints }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_StorageAccountMicrosoftEndpoints[name]))}
Swagger.field_name(::Type{ StorageAccountMicrosoftEndpoints }, property_name::Symbol) =  _property_map_StorageAccountMicrosoftEndpoints[property_name]

function check_required(o::StorageAccountMicrosoftEndpoints)
    true
end

function validate_property(::Type{ StorageAccountMicrosoftEndpoints }, name::Symbol, val)
end


text:
struct Sub
    # sourcenode::Union{Nothing, AbstractENode}
    ids::Vector{EClassId}
    nodes::Vector{Union{Nothing, ENodeLiteral}}
end

haseclassid(sub::Sub, p::PatVar) = sub.ids[p.idx] >= 0
geteclassid(sub::Sub, p::PatVar) = sub.ids[p.idx]

hasliteral(sub::Sub, p::PatVar) = sub.nodes[p.idx] !== nothing
getliteral(sub::Sub, p::PatVar) = sub.nodes[p.idx]

## ====================== Instantiation =======================

function instantiate(g::EGraph, pat::PatVar, sub::Sub, rule::AbstractRule)
    if haseclassid(sub, pat)
        ec = geteclass(g, geteclassid(sub, pat))
        if hasliteral(sub, pat)
            node = getliteral(sub, pat)
            return node.value
        end
        return ec
    else
        error("unbound pattern variable $pat in rule $rule")
    end
end

function instantiate(g::EGraph, pat::PatLiteral{T}, sub::Sub, rule::AbstractRule) where T
    pat.val
end

function instantiate(g::EGraph, pat::PatTypeAssertion, sub::Sub, rule::AbstractRule)
    instantiate(g, pat.name, sub, rule)
end

function instantiate(g::EGraph, pat::PatTerm, sub::Sub, rule::AbstractRule)
    eh = exprhead(pat)
    op = operation(pat)
    ar = arity(pat)

    T = gettermtype(g, op, ar)
    children = map(x -> instantiate(g, x, sub, rule), arguments(pat))
    similarterm(T, op, children; exprhead=eh)
end


text:
##################
# Trie #
##################

struct Trie{N}
    # Can be any nested (key => value) mapping, e.g. nested dict, static dict, named tuples, tuples, arrays or a combination
    nodes::N
end
Trie() = Trie(Dict{Any,Any}())
function Base.keys(t::Trie)
    return TrieMultiKey(StaticDict(k => t isa Trie ? keys(t) : nothing for (k, t) in t.nodes))
end

function Base.getindex(trie::Trie, key::TrieMultiKey)
    return Trie(
        StaticDict(
            map(keys(key.keys)) do k1
                k2 = key.keys[k1]
                k1 => (k2 === nothing ? trie.nodes[k1] : trie.nodes[k1][k2])
            end
        )
    )
end

Base.getindex(trie::Trie, key) = trie[TrieKey(key)]
function Base.getindex(trie::Trie, key::AbstractTrieKey)
    temp = get(trie.nodes, first(key), missing)
    (temp === missing || length(key) == 1) && return temp
    return temp[Base.tail(key)]
end

Base.setindex!(trie::Trie, value, key) = trie[TrieKey(key)] = value
function Base.setindex!(trie::Trie, value, key::AbstractTrieKey)
    if length(key) == 1
        setindex!(trie.nodes, value, first(key))
    else
        _trie = get!(trie.nodes, first(key), Trie{Any}())
        setindex!(_trie, value, Base.tail(key))
    end
    return value
end

function dummy(::Type{Trie})
    return Trie(
        Dict{Any,Any}(
            :x => rand(),
            :y => Trie(
                Dict{Any, Any}(
                    1 => rand(),
                    2 => rand(),
                )
            ),
        )
    )
end

import JSON
Base.show(io::IO, trie::Trie) = print(io, JSON.json(trie, 2))
Base.println(trie::Trie) = println(JSON.json(trie, 2))
Base.print(trie::Trie) = print(JSON.json(trie, 2))

function StaticDict(trie::Trie)
    keys, values = _dict(trie)
    return StaticDict(StaticTrieKey.(_static(keys)), _static(values))
end
function Base.Dict(trie::Trie)
    keys, values = _dict(trie)
    return Dict(TrieKey(k) => v for (k, v) in zip(keys, values))
end
@inline function _dict(trie, old_keys=(), old_values=(), header=())
    new_keys_values = map(collect(keys(trie.nodes))) do k
        if trie.nodes[k] isa Trie
            return _dict(trie.nodes[k], old_keys, old_values, (header..., k))
        else
            return (header..., k), trie.nodes[k]
        end
    end
    new_keys = map(kv -> kv[1], new_keys_values)
    new_values = map(kv -> kv[2], new_keys_values)
    return vcat(old_keys..., new_keys...), vcat(old_values..., new_values...)
end

export Trie, dummy

# Missing functionality
# - Dict to Trie
#      sort keys by length
#      find deps using subsume - single parent per node
# - Reshaping
#      shape function and reshape with a shape arg - shape defaults to size for arrays but can be anything
# - Remove an element - delete
# - Linearise
# - De-linearise
# - Merge
# - getindex with a bunch of keys
# - type specialise
# - explore https://github.com/andyferris/Dictionaries.jl
# - isempty
# - merge


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

import JSON

"""
    to_graphviz(network)

* `network::SymbolicNode`: the network to visualize.
* `title::AbstractString:` keyword argument, default "Network Visualization",
          the title of the GraphViz graph.
* `input_shapes`: keyword argument, default `nothing`. If provided,
          will run shape inference and plot with the shape information. Should
          be either a dictionary of name-shape mapping or an array of shapes.

Returns the graph description in GraphViz `dot` language.
"""
function to_graphviz(network :: SymbolicNode; title="Network Visualization", input_shapes=nothing)
  if !isa(input_shapes, Void)
    internals = get_internals(network)
    if isa(input_shapes, Dict)
      _, out_shapes, _ = infer_shape(internals; input_shapes...)
    else
      _, out_shapes, _ = infer_shape(internals, input_shapes...)
    end
    @assert(!isa(out_shapes, Void), "Failed to do shape inference, input shapes are incomplete")
    shape_dict = Dict(zip(list_outputs(internals), out_shapes))
    draw_shape = true
  else
    draw_shape = false
  end

  conf = JSON.parse(to_json(network))
  nodes = conf["nodes"]
  heads = unique([x[1]+1 for x in conf["heads"]])
  node_attr = Dict(:shape => :box, :fixedsize => true, :width => 1.3,
                   :height => 0.8034, :style => (:rounded, :filled), :penwidth => 2)
  io = IOBuffer()
  println(io, "digraph $(_simple_escape(title)) {")
  println(io, "node [fontsize=10];")
  println(io, "edge [fontsize=10];")

  # color map
  fillcolors = ("#8dd3c7", "#fb8072", "#ffffb3", "#bebada", "#80b1d3",
                "#fdb462", "#b3de69", "#fccde5")
  edgecolors = ("#245b51", "#941305", "#999900", "#3b3564", "#275372",
                "#975102", "#597d1c", "#90094e")

  # make nodes
  for i = 1:length(nodes)
    node  = nodes[i]
    op    = node["op"]
    name  = node["name"]
    attr  = deepcopy(node_attr)
    label = op

    # Up to 0.11.0 version of mxnet additional info was stored in
    # node["attr"]. Staring from 0.12 `attr` was changed to `attrs`.
    # See: https://github.com/dmlc/nnvm/pull/152
    if haskey(node, "attrs")
      node_info = node["attrs"]
    elseif haskey(node, "attr")
      node_info = node["attr"]
    end

    if op == "null"
      if i ∈ heads
        # heads are output nodes
        label = node["name"]
        colorkey = 1
      else
        # otherwise, input nodes, might be data, label or parameters
        continue
      end
    elseif op == "Convolution"
      if haskey(node_info,"stride")
        stride_info=_extract_shape(node_info["stride"])
      else
        stride_info="1"
      end

      label = format("Convolution\nkernel={1}\nstride={2}\nn-filter={3}",
                     _extract_shape(node_info["kernel"]),
                     stride_info,
                     node_info["num_filter"])
      colorkey = 2
    elseif op == "FullyConnected"
      label = format("FullyConnected\nnum-hidden={1}", node_info["num_hidden"])
      colorkey = 2
    elseif op == "Activation"
      label = format("Activation\nact-type={1}", node_info["act_type"])
      colorkey = 3
    elseif op == "BatchNorm"
      colorkey = 4
    elseif op == "Pooling"
      if haskey(node_info,"stride")
        stride_info=_extract_shape(node_info["stride"])
      else
        stride_info="1"
      end
      label = format("Pooling\ntype={1}\nkernel={2}\nstride={3}",
                     node_info["pool_type"],
                     _extract_shape(node_info["kernel"]),
                     stride_info)
      colorkey = 5
    elseif op ∈ ("Concat", "Flatten", "Reshape")
      colorkey = 6
    elseif endswith(op, "Output") || op == "BlockGrad"
      colorkey = 7
    else
      colorkey = 8
    end

    if op != "null"
      label = "$name\n$label"
    end
    attr[:fillcolor] = fillcolors[colorkey]
    attr[:color]     = edgecolors[colorkey]
    attr[:label]     = label
    _format_graphviz_node(io, name, attr)
  end

  # add edges
  for i = 1:length(nodes)
    node  = nodes[i]
    op    = node["op"]
    name  = node["name"]
    if op == "null"
      continue
    end
    inputs = node["inputs"]
    for item in inputs
      input_node = nodes[item[1]+1]
      input_name = input_node["name"]
      if input_node["op"] != "null" || (item[1]+1) ∈ heads
        attr = Dict(:dir => :back, :arrowtail => :open, :color => "#737373")
        if draw_shape
          if input_node["op"] != "null"
            key   = Symbol(input_name, "_output")
            shape = shape_dict[key][1:end-1]
          else
            key   = Symbol(input_name)
            shape = shape_dict[key][1:end-1]
          end
          label = "(" * join([string(x) for x in shape], ",") * ")"
          attr[:label] = label
        end
        _format_graphviz_edge(io, name, input_name, attr)
      end
    end
  end
  println(io, "}")

  return String(take!(io))
end

function _format_graphviz_attr(io::IOBuffer, attrs)
  label = get(attrs, :label, nothing)
  if isa(label, Void)
    print(io, " [")
  else
    print(io, " [label=$(_simple_escape(label)),")
  end
  first_attr = true
  for (k,v) in attrs
    if k != :label
      if !first_attr
        print(io, ",")
      end
      first_attr = false

      if isa(v, AbstractString) && v[1] == '#'
        # color
        v = _simple_escape(v)
      elseif isa(v, Tuple)
        v = _simple_escape(join([string(x) for x in v], ","))
      end
      print(io, "$k=$v")
    end
  end
  println(io, "];")
end
function _simple_escape(str)
  str = replace(string(str), r"\n", "\\n")
  return "\"$str\""
end
function _format_graphviz_node(io::IOBuffer, name::AbstractString, attrs)
  print(io, "$(_simple_escape(name)) ")
  _format_graphviz_attr(io, attrs)
end
function _format_graphviz_edge(io::IOBuffer, head, tail, attrs)
  print(io, """$(_simple_escape(head)) -> $(_simple_escape(tail)) """)
  _format_graphviz_attr(io, attrs)
end
function _extract_shape(str :: AbstractString)
  shape = matchall(r"\d+", str)
  shape = reverse(shape) # JSON in libmxnet has reversed shape (column vs row majoring)
  return "(" * join(shape, ",") * ")"
end

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet


text:
# License is MIT: https://github.com/genkuroki/MyUtils.jl/blob/main/LICENSE

"""
    RemoveOneFrom(A::AbstractVector, k::Integer) <: AbstractVector

Equivalent to A[[1:k-1; k+1:end]] without unnecessary memory allocation.

Example

```julia
A = [1, 2, 3, 4, 5]
R = RemoveOneFrom(A, 3)
@show collect(R)
# -> collect(R) = [1, 2, 4, 5]
R[3] = 99
@show R
# -> R = [1, 2, 99, 5]
@show A
# -> A = [1, 2, 3, 99, 5]
```
"""
struct RemoveOneFrom{T, A<:AbstractVector{T}, K<:Integer} <: AbstractVector{T}
    A::A
    k::K
end

Base.size(x::RemoveOneFrom) = (size(x.A, 1) - 1,)
Base.getindex(x::RemoveOneFrom, i) = i < x.k ? x.A[i] : x.A[i+1]
Base.getindex(x::RemoveOneFrom, ::Colon) = collect(x)
Base.setindex!(x::RemoveOneFrom, v, i) = i < x.k ? (x.A[i] = v) : (x.A[i+1] = v)


text:
# lambda: uniform
function prior_pdf_lambda(lambda, lambdamin, lambdamax)
    return 1/(lambdamax - lambdamin)
end

# theta: uniform
function prior_pdf_theta(theta, thetamin, thetamax)
    l1 = thetamax[1] - thetamin[1]
    l2 = thetamax[2] - thetamin[2]
    return 1/(l1*l2)
end

thetamin = [0., 0.]
thetamax = [1., 2.]
lambdamin = -3.
lambdamax = 3.

pdf_lambda = lambda -> prior_pdf_lambda(lambda, lambdamin, lambdamax)
pdf_theta = theta -> prior_pdf_theta(theta, thetamin, thetamax)

param_prior_lambda = (pdf = pdf_lambda, min = lambdamin, max = lambdamax)
param_prior_theta = (pdf = pdf_theta, min = thetamin, max = thetamax)

param_priorInfo = (lambda = param_prior_lambda, theta = param_prior_theta)

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
struct ProdKernel{K1<:Kernel, K2<:Kernel} <: PairKernel{K1,K2}
    kleft::K1
    kright::K2
end
leftkern(prodkern::ProdKernel) = prodkern.kleft
rightkern(prodkern::ProdKernel) = prodkern.kright

get_param_names(prodkern::ProdKernel) = composite_param_names(components(prodkern), :pk)

function cov(sk::ProdKernel, x::AbstractVector, y::AbstractVector)
    cov(sk.kleft, x, y) * cov(sk.kright, x, y)
end

@inline cov_ij(k::ProdKernel, X1::AbstractMatrix, X2::AbstractMatrix, data::EmptyData, i::Int, j::Int, dim::Int) = cov_ij(k.kleft, X1, X2, data, i, j, dim) * cov_ij(k.kright, X1, X2, data, i, j, dim)
@inline cov_ij(k::ProdKernel, X1::AbstractMatrix, X2::AbstractMatrix, data::PairData, i::Int, j::Int, dim::Int) = cov_ij(k.kleft, X1, X2, data.data1, i, j, dim) * cov_ij(k.kright, X1, X2, data.data2, i, j, dim)

@inline function dKij_dθp(prodkern::ProdKernel, X1::AbstractMatrix, X2::AbstractMatrix, data::EmptyData, i::Int, j::Int, p::Int, dim::Int)
    np = num_params(prodkern.kleft)
    if p<=np
        cK_other = cov_ij(prodkern.kright, X1, X2, data, i, j, dim)
        return dKij_dθp(prodkern.kleft,    X1, X2, data, i,j,p,dim)*cK_other
    else
        cK_other = cov_ij(prodkern.kleft, X1, X2, data, i, j, dim)
        return dKij_dθp(prodkern.kright,  X1, X2, data, i,j,p-np,dim)*cK_other
    end
end
@inline function dKij_dθp(prodkern::ProdKernel, X1::AbstractMatrix, X2::AbstractMatrix, data::PairData, i::Int, j::Int, p::Int, dim::Int)
    np = num_params(prodkern.kleft)
    if p<=np
        cK_other = cov_ij(prodkern.kright,  X1, X2, data.data2, i, j, dim)
        dKij_sub = dKij_dθp(prodkern.kleft, X1, X2, data.data1, i, j, p, dim)
        return dKij_sub * cK_other
    else
        cK_other = cov_ij(prodkern.kleft,    X1, X2, data.data1, i, j, dim)
        dKij_sub = dKij_dθp(prodkern.kright, X1, X2, data.data2, i, j, p-np, dim)
        return dKij_sub * cK_other
    end
end
@inline @inbounds function dKij_dθ!(dK::AbstractVector, prodkern::ProdKernel, X1::AbstractMatrix, X2::AbstractMatrix,
                                    data::EmptyData, i::Int, j::Int, dim::Int, npars::Int)
    cov_left  = cov_ij(prodkern.kleft,  X1, X2, data, i, j, dim)
    cov_right = cov_ij(prodkern.kright, X1, X2, data, i, j, dim)
    npright = num_params(prodkern.kright)
    npleft = num_params(prodkern.kleft)
    dKij_dθ!(dK, prodkern.kright, X1, X2, data, i, j, dim, npars-npleft)
    for ipar in npright:-1:1
        dK[npleft+ipar] = dK[ipar]*cov_left
    end
    dKij_dθ!(dK,  prodkern.kleft,  X1, X2, data, i, j, dim, npleft)
    for ipar in 1:npleft
        dK[ipar] *= cov_right
    end
end
@inline @inbounds function dKij_dθ!(dK::AbstractVector, prodkern::ProdKernel, X1::AbstractMatrix, X2::AbstractMatrix,
                                    data::PairData, i::Int, j::Int, dim::Int, npars::Int)
    cov_left  = cov_ij(prodkern.kleft,  X1, X2, data.data1, i, j, dim)
    cov_right = cov_ij(prodkern.kright, X1, X2, data.data2, i, j, dim)
    npright = num_params(prodkern.kright)
    npleft = num_params(prodkern.kleft)
    dKij_dθ!(dK, prodkern.kright, X1, X2, data.data2, i, j, dim, npars-npleft)
    for ipar in npright:-1:1
        dK[npleft+ipar] = dK[ipar]*cov_left
    end
    dKij_dθ!(dK,  prodkern.kleft,  X1, X2, data.data1, i, j, dim, npleft)
    for ipar in 1:npleft
        dK[ipar] *= cov_right
    end
end

# Multiplication operators
Base.:*(kleft::Kernel, kright::Kernel) = ProdKernel(kleft,kright)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to the API function to compose rotations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/compose_rotations.jl
# ================================

# Functions: compose_rotations
# ----------------------------

@testset "Compose rotations" begin
    for T in (Float32, Float64)
        # DCMs
        # ======================================================================

        # Sample 4 DCMs.
        D1 = rand(DCM{T})
        D2 = rand(DCM{T})
        D3 = rand(DCM{T})
        D4 = rand(DCM{T})

        # Test the function `compose_rotation`.
        @test                D1 === compose_rotation(D1)
        @test           D2 * D1  ≈  compose_rotation(D1,D2)
        @test      D3 * D2 * D1  ≈  compose_rotation(D1,D2,D3)
        @test D4 * D3 * D2 * D1  ≈  compose_rotation(D1,D2,D3,D4)

        # Euler angle and axis
        # ==========================================================================

        # Sample 4 Euler angle and axis.
        ea1 = rand(EulerAngleAxis{T})
        ea2 = rand(EulerAngleAxis{T})
        ea3 = rand(EulerAngleAxis{T})
        ea4 = rand(EulerAngleAxis{T})

        ear1 = ea2 * ea1
        eac1 = compose_rotation(ea1, ea2)

        ear2 = ea3 * ea2 * ea1
        eac2 = compose_rotation(ea1, ea2, ea3)

        ear3 = ea4 * ea3 * ea2 * ea1
        eac3 = compose_rotation(ea1, ea2, ea3, ea4)

        # Test the function `compose_rotation`.
        @test  ea1 === compose_rotation(ea1)
        @test ear1  ≈  eac1
        @test ear2  ≈  eac2
        @test ear3  ≈  eac3

        # Euler angles
        # ==========================================================================

        # Sample 4 Euler angles.
        Θ1 = rand(EulerAngles{T})
        Θ2 = rand(EulerAngles{T})
        Θ3 = rand(EulerAngles{T})
        Θ4 = rand(EulerAngles{T})

        Θr1 = Θ2 * Θ1
        Θc1 = compose_rotation(Θ1, Θ2)

        Θr2 = Θ3 * Θ2 * Θ1
        Θc2 = compose_rotation(Θ1, Θ2, Θ3)

        Θr3 = Θ4 * Θ3 * Θ2 * Θ1
        Θc3 = compose_rotation(Θ1, Θ2, Θ3, Θ4)

        # Test the function `compose_rotation`.
        @test  Θ1 === compose_rotation(Θ1)
        @test Θr1  ≈  Θc1
        @test Θr2  ≈  Θc2
        @test Θr3  ≈  Θc3

        # Quaternions
        # ==========================================================================

        # Sample 4 quaternions.
        q1 = rand(Quaternion{T})
        q2 = rand(Quaternion{T})
        q3 = rand(Quaternion{T})
        q4 = rand(Quaternion{T})

        # Test the function `compose_rotation`.
        @test               q1 === compose_rotation(q1)
        @test       (q1*q2)[:]  ≈  compose_rotation(q1,q2)[:]
        @test    (q1*q2*q3)[:]  ≈  compose_rotation(q1,q2,q3)[:]
        @test (q1*q2*q3*q4)[:]  ≈  compose_rotation(q1,q2,q3,q4)[:]
    end
end

# Operator ∘
# ----------

@testset "Operator ∘" begin
    for T in (Float32, Float64)
        D  = rand(DCM{T})
        ea = rand(EulerAngles{T})
        av = rand(EulerAngleAxis{T})
        q  = rand(Quaternion{T})

        D_ea = convert(DCM, ea)
        D_av = convert(DCM, av)
        D_q  = convert(DCM, q)

        # DCM
        # ======================================================================

        R = D ∘ ea ∘ av ∘ q
        R_exp = D * D_ea * D_av * D_q
        @test R ≈ R_exp atol = √(eps(T))

        # Euler angle and axis
        # ======================================================================

        R = av ∘ D ∘ ea ∘ q
        R_exp = convert(EulerAngleAxis, D_av * D * D_ea * D_q)
        @test R ≈ R_exp atol = 100 * √(eps(T))

        # Euler angles
        # ======================================================================

        R = ea ∘ D ∘ av ∘ q
        R_exp = convert(EulerAngles(R.rot_seq), D_ea * D * D_av * D_q)

        # When converting to Euler angles, we normalize the angles. We need to
        # make the same with the composition for the sake of testing.
        R = convert(EulerAngles, R)

        @test R ≈ R_exp atol = 100 * √(eps(T))

        # Quaternion
        # ======================================================================

        R = q ∘ ea ∘ av ∘ D
        R_exp = convert(Quaternion, D_q * D_ea * D_av * D)

        # When converting to quaternion, we make sure that the real part is
        # always positive. Hence, we need to do the same with the composition
        # for the sake of testing.
        if R.q0 < 0
            R = -R
        end

        @test R ≈ R_exp atol = √(eps(T))
    end
end