text:
import MultivariatePolynomials
import MultivariatePolynomials: AbstractPolynomialLike

# check if a matrix is square
@inline function issquare(A::AbstractMatrix)
    return size(A, 1) == size(A, 2)
end

for (Z, AZ) in ((:ContinuousIdentitySystem, :AbstractContinuousSystem),
                (:DiscreteIdentitySystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z) <: $(AZ)
            statedim::Int
        end
        statedim(s::$Z) = s.statedim
        inputdim(s::$Z) = 0
        noisedim(s::$Z) = 0
        state_matrix(s::$Z) = I(s.statedim)
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = true
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = false
            isconstrained(::$T) = false
        end
    end
end

@doc """
    ContinuousIdentitySystem <: AbstractContinuousSystem

Trivial identity continuous-time system of the form:

```math
    x(t)' = 0 \\; \\forall t.
```

### Fields

- `statedim` -- number of state variables
"""
ContinuousIdentitySystem

@doc """
    DiscreteIdentitySystem <: AbstractDiscreteSystem

Trivial identity discrete-time system of the form:
```math
    x_{k+1} = x_k \\; \\forall k.
```

### Fields

- `statedim` -- number of state variables
"""
DiscreteIdentitySystem

for (Z, AZ) in ((:ConstrainedContinuousIdentitySystem, :AbstractContinuousSystem),
                (:ConstrainedDiscreteIdentitySystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){ST} <: $(AZ)
            statedim::Int
            X::ST
        end
        statedim(s::$Z) = s.statedim
        inputdim(::$Z) = 0
        noisedim(::$Z) = 0
        stateset(s::$Z) = s.X
        state_matrix(s::$Z) = I(s.statedim)
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = true
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = false
            isconstrained(::$T) = true
        end
    end
end

@doc """
    ConstrainedContinuousIdentitySystem <: AbstractContinuousSystem

Trivial identity continuous-time system with domain constraints of the form:

```math
    x(t)' = 0, \\; x(t) ∈ \\mathcal{X} \\; \\forall t.
```

### Fields

- `statedim` -- number of state variables
- `X`        -- state constraints
"""
ConstrainedContinuousIdentitySystem

@doc """
    ConstrainedDiscreteIdentitySystem <: AbstractDiscreteSystem

Trivial identity discrete-time system with domain constraints of the form:

```math
    x_{k+1} = x_k, \\; x_k ∈ \\mathcal{X} \\; \\forall k.
```

### Fields

- `statedim` -- number of state variables
- `X`        -- state constraints
"""
ConstrainedDiscreteIdentitySystem

for (Z, AZ) in ((:LinearContinuousSystem, :AbstractContinuousSystem),
                (:LinearDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MT <: AbstractMatrix{T}} <: $(AZ)
            A::MT
            function $(Z)(A::MT) where {T, MT <: AbstractMatrix{T}}
                @assert issquare(A)
                return new{T, MT}(A)
            end
        end
        function $(Z)(A::Number)
             return $(Z)(hcat(A))
        end

        statedim(s::$Z) = size(s.A, 1)
        inputdim(::$Z) = 0
        noisedim(::$Z) = 0
        state_matrix(s::$Z) = s.A
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = true
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = false
            isconstrained(::$T) = false
        end
    end
end

@doc """
    LinearContinuousSystem

Continuous-time linear system of the form:

```math
    x(t)' = A x(t) \\; \\forall t.
```

### Fields

- `A` -- state matrix
"""
LinearContinuousSystem

@doc """
    LinearDiscreteSystem

Discrete-time linear system of the form:

```math
    x_{k+1} = A x_k \\; \\forall k.
```

### Fields

- `A` -- state matrix
"""
LinearDiscreteSystem

for (Z, AZ) in ((:AffineContinuousSystem, :AbstractContinuousSystem),
                (:AffineDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MT <: AbstractMatrix{T}, VT <: AbstractVector{T}} <: $(AZ)
            A::MT
            c::VT
            function $(Z)(A::MT, c::VT) where {T, MT <: AbstractMatrix{T}, VT <: AbstractVector{T}}
                @assert checksquare(A) == length(c)
                return new{T, MT, VT}(A, c)
            end
        end
        function $(Z)(A::Number, c::Number)
             return $(Z)(hcat(A), vcat(c))
        end

        statedim(s::$Z) = length(s.c)
        inputdim(::$Z) = 0
        noisedim(::$Z) = 0
        state_matrix(s::$Z) = s.A
        affine_term(s::$Z) = s.c
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = false
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = false
            isconstrained(::$T) = false
        end
    end
end

@doc """
    AffineContinuousSystem

Continuous-time affine system of the form:

```math
    x(t)' = A x(t) + c \\; \\forall t.
```

### Fields

- `A` -- state matrix
- `c` -- affine term
"""
AffineContinuousSystem

@doc """
    AffineDiscreteSystem

Discrete-time affine system of the form:

```math
    x_{k+1} = A x_k + c \\; \\forall k.
```

### Fields

- `A` -- state matrix
- `c` -- affine term
"""
AffineDiscreteSystem

for (Z, AZ) in ((:LinearControlContinuousSystem, :AbstractContinuousSystem),
                (:LinearControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}} <: $(AZ)
            A::MTA
            B::MTB
            function $(Z)(A::MTA, B::MTB) where {T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}}
                @assert checksquare(A) == size(B, 1)
                return new{T, MTA, MTB}(A, B)
            end
        end
        function $(Z)(A::Number, B::Number)
             return $(Z)(hcat(A), hcat(B))
        end

        statedim(s::$Z) = size(s.A, 1)
        inputdim(s::$Z) = size(s.B, 2)
        noisedim(::$Z) = 0
        state_matrix(s::$Z) = s.A
        input_matrix(s::$Z) = s.B
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = true
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = true
            isconstrained(::$T) = false
        end
    end
end

@doc """
    LinearControlContinuousSystem

Continuous-time linear control system of the form:

```math
    x(t)' = A x(t) + B u(t) \\; \\forall t.
```

### Fields

- `A` -- state matrix
- `B` -- input matrix
"""
LinearControlContinuousSystem

@doc """
    LinearControlDiscreteSystem

Discrete-time linear control system of the form:

```math
    x_{k+1} = A x_k + B u_k \\; \\forall k.
```

### Fields

- `A` -- state matrix
- `B` -- input matrix
"""
LinearControlDiscreteSystem

for (Z, AZ) in ((:ConstrainedLinearContinuousSystem, :AbstractContinuousSystem),
                (:ConstrainedLinearDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MT <: AbstractMatrix{T}, ST} <: $(AZ)
            A::MT
            X::ST
            function $(Z)(A::MT, X::ST) where {T, MT <: AbstractMatrix{T}, ST}
                @assert issquare(A)
                return new{T, MT, ST}(A, X)
            end
        end
        function $(Z)(A::Number, X)
             return $(Z)(hcat(A), X)
        end

        statedim(s::$Z) = size(s.A, 1)
        inputdim(::$Z) = 0
        noisedim(::$Z) = 0
        stateset(s::$Z) = s.X
        state_matrix(s::$Z) = s.A
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = true
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = false
            isconstrained(::$T) = true
        end
    end
end

@doc """
    ConstrainedLinearContinuousSystem

Continuous-time linear system with domain constraints of the form:

```math
    x(t)' = A x(t), \\; x(t) ∈ \\mathcal{X} \\; \\forall t.
```

### Fields

- `A` -- state matrix
- `X` -- state constraints
"""
ConstrainedLinearContinuousSystem

@doc """
    ConstrainedLinearDiscreteSystem

Discrete-time linear system with domain constraints of the form:

```math
    x_{k+1} = A x_k, \\; x_k ∈ \\mathcal{X} \\; \\forall k.
```

### Fields

- `A` -- state matrix
- `X` -- state constraints
"""
ConstrainedLinearDiscreteSystem

for (Z, AZ) in ((:ConstrainedAffineContinuousSystem, :AbstractContinuousSystem),
                (:ConstrainedAffineDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MT <: AbstractMatrix{T}, VT <: AbstractVector{T}, ST} <: $(AZ)
            A::MT
            c::VT
            X::ST
            function $(Z)(A::MT, c::VT, X::ST) where {T, MT <: AbstractMatrix{T}, VT <: AbstractVector{T}, ST}
                @assert checksquare(A) == length(c)
                return new{T, MT, VT, ST}(A, c, X)
            end
        end
        function $(Z)(A::Number, c::Number, X)
             return $(Z)(hcat(A), vcat(c), X)
        end

        statedim(s::$Z) = length(s.c)
        inputdim(::$Z) = 0
        noisedim(::$Z) = 0
        stateset(s::$Z) = s.X
        state_matrix(s::$Z) = s.A
        affine_term(s::$Z) = s.c
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = false
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = false
            isconstrained(::$T) = true
        end
    end
end

@doc """
    ConstrainedAffineContinuousSystem

Continuous-time affine system with domain constraints of the form:

```math
    x(t)' = A x(t) + c, \\; x(t) ∈ \\mathcal{X} \\; \\forall t.
```

### Fields

- `A` -- state matrix
- `c` -- affine term
- `X` -- state constraints
"""
ConstrainedAffineContinuousSystem

@doc """
    ConstrainedAffineDiscreteSystem

Discrete-time affine system with domain constraints of the form:

```math
    x_{k+1} = A x_k + c, \\; x_k ∈ \\mathcal{X} \\; \\forall k.
```

### Fields

- `A` -- state matrix
- `c` -- affine term
- `X` -- state constraints
"""
ConstrainedAffineDiscreteSystem

 for (Z, AZ) in ((:AffineControlContinuousSystem, :AbstractContinuousSystem),
                 (:AffineControlDiscreteSystem, :AbstractDiscreteSystem))
     @eval begin
         struct $(Z){T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, VT <: AbstractVector{T}} <: $(AZ)
             A::MTA
             B::MTB
             c::VT
             function $(Z)(A::MTA, B::MTB, c::VT) where {T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, VT <: AbstractVector{T}}
                 @assert checksquare(A) == length(c) == size(B, 1)
                 return new{T, MTA, MTB, VT}(A, B, c)
             end
         end
         function $(Z)(A::Number, B::Number, c::Number)
              return $(Z)(hcat(A), hcat(B), vcat(c))
         end

         statedim(s::$Z) = length(s.c)
         inputdim(s::$Z) = size(s.B, 2)
         noisedim(::$Z) = 0
         state_matrix(s::$Z) = s.A
         input_matrix(s::$Z) = s.B
         affine_term(s::$Z) = s.c
     end
     for T in [Z, Type{<:eval(Z)}]
         @eval begin
             islinear(::$T) = false
             isaffine(::$T) = true
             ispolynomial(::$T) = false
             isblackbox(::$T) = false
             isnoisy(::$T) = false
             iscontrolled(::$T) = true
             isconstrained(::$T) = false
         end
     end
 end

 @doc """
     AffineControlContinuousSystem

 Continuous-time affine control system of the form:

 ```math
     x(t)' = A x(t) + B u(t) + c, \\; \\forall t.
 ```

 ### Fields

 - `A` -- state matrix
 - `B` -- input matrix
 - `c` -- affine term
 """
 AffineControlContinuousSystem

 @doc """
     AffineControlDiscreteSystem

 Continuous-time affine control system of the form:

 ```math
     x_{k+1} = A x_k + B u_k + c, \\; \\forall k.
 ```

 ### Fields

 - `A` -- state matrix
 - `B` -- input matrix
 - `c` -- affine term
 """
 AffineControlDiscreteSystem

for (Z, AZ) in ((:ConstrainedAffineControlContinuousSystem, :AbstractContinuousSystem),
                (:ConstrainedAffineControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, VT <: AbstractVector{T}, ST, UT} <: $(AZ)
            A::MTA
            B::MTB
            c::VT
            X::ST
            U::UT
            function $(Z)(A::MTA, B::MTB, c::VT, X::ST, U::UT) where {T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, VT <: AbstractVector{T}, ST, UT}
                @assert checksquare(A) == length(c) == size(B, 1)
                return new{T, MTA, MTB, VT, ST, UT}(A, B, c, X, U)
            end
        end
        function $(Z)(A::Number, B::Number, c::Number, X, U)
             return $(Z)(hcat(A), hcat(B), vcat(c), X, U)
        end

        statedim(s::$Z) = length(s.c)
        inputdim(s::$Z) = size(s.B, 2)
        noisedim(::$Z) = 0
        stateset(s::$Z) = s.X
        inputset(s::$Z) = s.U
        state_matrix(s::$Z) = s.A
        input_matrix(s::$Z) = s.B
        affine_term(s::$Z) = s.c
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = false
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = true
            isconstrained(::$T) = true
        end
    end
end

@doc """
    ConstrainedAffineControlContinuousSystem

Continuous-time affine control system with domain constraints of the form:

```math
    x(t)' = A x(t) + B u(t) + c, \\; x(t) ∈ \\mathcal{X}, \\; u(t) ∈ \\mathcal{U} \\; \\forall t.
```

### Fields

- `A` -- state matrix
- `B` -- input matrix
- `c` -- affine term
- `X` -- state constraints
- `U` -- input constraints
"""
ConstrainedAffineControlContinuousSystem

@doc """
    ConstrainedAffineControlDiscreteSystem

Continuous-time affine control system with domain constraints of the form:

```math
    x_{k+1} = A x_k + B u_k + c, \\; x_k ∈ \\mathcal{X}, \\; u_k ∈ \\mathcal{U} \\; \\forall k.
```

### Fields

- `A` -- state matrix
- `B` -- input matrix
- `c` -- affine term
- `X` -- state constraints
- `U` -- input constraints
"""
ConstrainedAffineControlDiscreteSystem

for (Z, AZ) in ((:ConstrainedLinearControlContinuousSystem, :AbstractContinuousSystem),
                (:ConstrainedLinearControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, ST, UT} <: $(AZ)
            A::MTA
            B::MTB
            X::ST
            U::UT
            function $(Z)(A::MTA, B::MTB, X::ST, U::UT) where {T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, ST, UT}
                @assert checksquare(A) == size(B, 1)
                return new{T, MTA, MTB, ST, UT}(A, B, X, U)
            end
        end
        function $(Z)(A::Number, B::Number, X, U)
             return $(Z)(hcat(A), hcat(B), X, U)
        end

        statedim(s::$Z) = size(s.A, 1)
        inputdim(s::$Z) = size(s.B, 2)
        noisedim(::$Z) = 0
        stateset(s::$Z) = s.X
        inputset(s::$Z) = s.U
        state_matrix(s::$Z) = s.A
        input_matrix(s::$Z) = s.B
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = true
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = true
            isconstrained(::$T) = true
        end
    end
end

@doc """
    ConstrainedLinearControlContinuousSystem

Continuous-time linear control system with domain constraints of the form:
```math
    x(t)' = A x(t) + B u(t), \\; x(t) ∈ \\mathcal{X}, \\; u(t) ∈ \\mathcal{U} \\; \\forall t.
```

### Fields

- `A` -- state matrix
- `B` -- input matrix
- `X` -- state constraints
- `U` -- input constraints
"""
ConstrainedLinearControlContinuousSystem

@doc """
    ConstrainedLinearControlDiscreteSystem

Discrete-time linear control system with domain constraints of the form:

```math
    x_{k+1} = A x_k + B u_k, \\; x_k ∈ \\mathcal{X}, \\; u_k ∈ \\mathcal{U} \\; \\forall k.
```

### Fields

- `A` -- state matrix
- `B` -- input matrix
- `X` -- state constraints
- `U` -- input constraints
"""
ConstrainedLinearControlDiscreteSystem

for (Z, AZ) in ((:LinearAlgebraicContinuousSystem, :AbstractContinuousSystem),
                (:LinearAlgebraicDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTA <: AbstractMatrix{T}, MTE <: AbstractMatrix{T}} <: $(AZ)
            A::MTA
            E::MTE
            function $(Z)(A::MTA, E::MTE) where {T, MTA <: AbstractMatrix{T}, MTE <: AbstractMatrix{T}}
                @assert size(A) == size(E)
                return new{T, MTA, MTE}(A, E)
            end
        end
        function $(Z)(A::Number, E::Number)
             return $(Z)(hcat(A), hcat(E))
        end

        statedim(s::$Z) = size(s.A, 1)
        inputdim(::$Z) = 0
        noisedim(::$Z) = 0
        state_matrix(s::$Z) = s.A
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = true
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = false
            isconstrained(::$T) = false
        end
    end
end

@doc """
    LinearAlgebraicContinuousSystem

Continuous-time linear algebraic system of the form:

```math
    E x(t)' = A x(t) \\; \\forall t.
```

### Fields

- `A` -- state matrix
- `E` -- matrix, same size as `A`
"""
LinearAlgebraicContinuousSystem

@doc """
    LinearAlgebraicDiscreteSystem

Discrete-time linear algebraic system of the form:

```math
    E x_{k+1} = A x_k \\; \\forall k.
```

### Fields

- `A` -- state matrix
- `E` -- matrix, same size as `A`
"""
LinearAlgebraicDiscreteSystem

for (Z, AZ) in ((:ConstrainedLinearAlgebraicContinuousSystem, :AbstractContinuousSystem),
                (:ConstrainedLinearAlgebraicDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTA <: AbstractMatrix{T}, MTE <: AbstractMatrix{T}, ST} <: $(AZ)
            A::MTA
            E::MTE
            X::ST
            function $(Z)(A::MTA, E::MTE, X::ST) where {T, MTA <: AbstractMatrix{T}, MTE <: AbstractMatrix{T}, ST}
                @assert size(A) == size(E)
                return new{T, MTA, MTE, ST}(A, E, X)
            end
        end
        function $(Z)(A::Number, E::Number, X)
             return $(Z)(hcat(A), hcat(E), X)
        end

        statedim(s::$Z) = size(s.A, 1)
        inputdim(::$Z) = 0
        noisedim(::$Z) = 0
        stateset(s::$Z) = s.X
        state_matrix(s::$Z) = s.A
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = true
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = false
            isconstrained(::$T) = true
        end
    end
end

@doc """
    ConstrainedLinearAlgebraicContinuousSystem

Continuous-time linear system with domain constraints of the form:

```math
    E x(t)' = A x(t), \\; x(t) ∈ \\mathcal{X} \\; \\forall t.
```

### Fields

- `A` -- state matrix
- `E` -- matrix, same size as `A`
- `X` -- state constraints
"""
ConstrainedLinearAlgebraicContinuousSystem

@doc """
    ConstrainedLinearAlgebraicDiscreteSystem

Discrete-time linear system with domain constraints of the form:

```math
    E x_{k+1} = A x_k, \\; x_k ∈ \\mathcal{X} \\; \\forall k.
```

### Fields

- `A` -- state matrix
- `E` -- matrix, same size as `A`
- `X` -- state constraints
"""
ConstrainedLinearAlgebraicDiscreteSystem

for (Z, AZ) in ((:PolynomialContinuousSystem, :AbstractContinuousSystem),
                (:PolynomialDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, PT <: AbstractPolynomialLike{T}, VPT <: AbstractVector{PT}} <: $(AZ)
            p::VPT
            statedim::Int
            function $(Z)(p::VPT, statedim::Int) where {T, PT <: AbstractPolynomialLike{T}, VPT <: AbstractVector{PT}}
                @assert statedim == MultivariatePolynomials.nvariables(p) "the state dimension $(statedim) does not match the number of state variables"
                return new{T, PT, VPT}(p, statedim)
            end
        end
        statedim(s::$Z) = s.statedim
        inputdim(::$Z) = 0
        noisedim(::$Z) = 0
        mapping(s::$Z) = s.p

        MultivariatePolynomials.variables(s::$Z) = MultivariatePolynomials.variables(s.p)
        MultivariatePolynomials.nvariables(s::$Z) = s.statedim

        $(Z)(p::AbstractVector{<:AbstractPolynomialLike}) = $(Z)(p, MultivariatePolynomials.nvariables(p))
        $(Z)(p::AbstractPolynomialLike) = $(Z)([p])
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = false
            isaffine(::$T) = false
            ispolynomial(::$T) = true
            isblackbox(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = false
            isconstrained(::$T) = false
        end
    end
end

@doc """
    PolynomialContinuousSystem

Continuous-time polynomial system of the form:

```math
    x(t)' = p(x(t)) \\; \\forall t.
```

### Fields

- `p`        -- polynomial vector field
- `statedim` -- number of state variables
"""
PolynomialContinuousSystem

@doc """
    PolynomialDiscreteSystem

Discrete-time polynomial system of the form:

```math
    x_{k+1} = p(x_k) \\; \\forall k.
```

### Fields

- `p`        -- polynomial vector field
- `statedim` -- number of state variables
"""
PolynomialDiscreteSystem

for (Z, AZ) in ((:ConstrainedPolynomialContinuousSystem, :AbstractContinuousSystem),
                (:ConstrainedPolynomialDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, PT <: AbstractPolynomialLike{T}, VPT <: AbstractVector{PT}, ST} <: $(AZ)
            p::VPT
            statedim::Int
            X::ST
            function $(Z)(p::VPT, statedim::Int, X::ST) where {T, PT <: AbstractPolynomialLike{T}, VPT <: AbstractVector{PT}, ST}
                @assert statedim == MultivariatePolynomials.nvariables(p) "the state dimension $(statedim) does not match the number of state variables"
                return new{T, PT, VPT, ST}(p, statedim, X)
            end
        end
        statedim(s::$Z) = s.statedim
        inputdim(::$Z) = 0
        noisedim(::$Z) = 0
        stateset(s::$Z) = s.X
        mapping(s::$Z) = s.p

        MultivariatePolynomials.variables(s::$Z) = MultivariatePolynomials.variables(s.p)
        MultivariatePolynomials.nvariables(s::$Z) = s.statedim

        $Z(p::AbstractVector{<:AbstractPolynomialLike}, X::ST) where {ST} = $(Z)(p, MultivariatePolynomials.nvariables(p), X)
        $Z(p::AbstractPolynomialLike, X::ST) where {ST} = $(Z)([p], X)
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = false
            isaffine(::$T) = false
            ispolynomial(::$T) = true
            isblackbox(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = false
            isconstrained(::$T) = true
        end
    end
end

@doc """
    ConstrainedPolynomialContinuousSystem

Continuous-time polynomial system with domain constraints:

```math
    x(t)' = p(x(t)), \\; x(t) ∈ \\mathcal{X} \\; \\forall t.
```

### Fields

- `p`        -- polynomial vector field
- `X`        -- constraint set
- `statedim` -- number of state variables
"""
ConstrainedPolynomialContinuousSystem

@doc """
    ConstrainedPolynomialDiscreteSystem

Discrete-time polynomial system with domain constraints:

```math
    x_{k+1} = p(x_k), \\;  x_k ∈ \\mathcal{X} \\; \\forall k.
```

### Fields

- `p`        -- polynomial
- `X`        -- constraint set
- `statedim` -- number of state variables
"""
ConstrainedPolynomialDiscreteSystem

for (Z, AZ) in ((:BlackBoxContinuousSystem, :AbstractContinuousSystem),
                (:BlackBoxDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){FT} <: $(AZ)
            f::FT
            statedim::Int
        end
        statedim(s::$Z) = s.statedim
        inputdim(s::$Z) = 0
        noisedim(::$Z) = 0
        mapping(s::$Z) = s.f
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = false
            isaffine(::$T) = false
            ispolynomial(::$T) = false
            isblackbox(::$T) = true
            isnoisy(::$T) = false
            iscontrolled(::$T) = false
            isconstrained(::$T) = false
        end
    end
end

@doc """
    BlackBoxContinuousSystem <: AbstractContinuousSystem

Continuous-time system defined by a right-hand side of the form:

```math
    x(t)' = f(x(t)) \\; \\forall t.
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
"""
BlackBoxContinuousSystem

@doc """
    BlackBoxDiscreteSystem <: AbstractDiscreteSystem

Discrete-time system defined by a right-hand side of the form:

```math
    x_{k+1} = f(x_k) \\; \\forall k.
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
"""
BlackBoxDiscreteSystem

for (Z, AZ) in ((:ConstrainedBlackBoxContinuousSystem, :AbstractContinuousSystem),
                (:ConstrainedBlackBoxDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){FT, ST} <: $(AZ)
            f::FT
            statedim::Int
            X::ST
        end
        statedim(s::$Z) = s.statedim
        inputdim(s::$Z) = 0
        noisedim(::$Z) = 0
        stateset(s::$Z) = s.X
        mapping(s::$Z) = s.f
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = false
            isaffine(::$T) = false
            ispolynomial(::$T) = false
            isblackbox(::$T) = true
            isnoisy(::$T) = false
            iscontrolled(::$T) = false
            isconstrained(::$T) = true
        end
    end
end

@doc """
    ConstrainedBlackBoxContinuousSystem <: AbstractContinuousSystem

Continuous-time system defined by a right-hand side with domain constraints of the
form:

```math
    x(t)' = f(x(t)), \\; x(t) ∈ \\mathcal{X} \\; \\forall t.
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
- `X`        -- state constraints
"""
ConstrainedBlackBoxContinuousSystem

@doc """
    ConstrainedBlackBoxDiscreteSystem <: AbstractDiscreteSystem

Discrete-time system defined by a right-hand side with domain constraints
of the form:

```math
    x_{k+1} = f(x_k), \\; x_k ∈ \\mathcal{X} \\; \\forall k.
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
- `X`        -- state constraints
"""
ConstrainedBlackBoxDiscreteSystem

for (Z, AZ) in ((:BlackBoxControlContinuousSystem, :AbstractContinuousSystem),
                (:BlackBoxControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){FT} <: $(AZ)
            f::FT
            statedim::Int
            inputdim::Int
        end
        statedim(s::$Z) = s.statedim
        inputdim(s::$Z) = s.inputdim
        noisedim(::$Z) = 0
        mapping(s::$Z) = s.f
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = false
            isaffine(::$T) = false
            ispolynomial(::$T) = false
            isblackbox(::$T) = true
            isnoisy(::$T) = false
            iscontrolled(::$T) = true
            isconstrained(::$T) = false
        end
    end
end

@doc """
    BlackBoxControlContinuousSystem <: AbstractContinuousSystem

Continuous-time control system defined by a right-hand side of the form:

```math
    x(t)' = f(x(t), u(t)) \\; \\forall t.
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
- `inputdim` -- number of input variables
"""
BlackBoxControlContinuousSystem

@doc """
    BlackBoxControlDiscreteSystem <: AbstractDiscreteSystem

Discrete-time control system defined by a right-hand side of the form:

```math
    x_{k+1} = f(x_k, u_k) \\; \\forall k.
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
- `inputdim` -- number of input variables
"""
BlackBoxControlDiscreteSystem

for (Z, AZ) in ((:ConstrainedBlackBoxControlContinuousSystem, :AbstractContinuousSystem),
                (:ConstrainedBlackBoxControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){FT, ST, UT} <: $(AZ)
            f::FT
            statedim::Int
            inputdim::Int
            X::ST
            U::UT
        end
        statedim(s::$Z) = s.statedim
        inputdim(s::$Z) = s.inputdim
        noisedim(::$Z) = 0
        stateset(s::$Z) = s.X
        inputset(s::$Z) = s.U
        mapping(s::$Z) = s.f
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = false
            isaffine(::$T) = false
            ispolynomial(::$T) = false
            isblackbox(::$T) = true
            isnoisy(::$T) = false
            iscontrolled(::$T) = true
            isconstrained(::$T) = true
        end
    end
end

@doc """
    ConstrainedBlackBoxControlContinuousSystem <: AbstractContinuousSystem

Continuous-time control system defined by a right-hand side with domain constraints
of the form:

```math
    x(t)' = f(x(t), u(t)), \\; x(t) ∈ \\mathcal{X}, \\; u(t) ∈ \\mathcal{U} \\; \\forall t.
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
- `inputdim` -- number of input variables
- `X`        -- state constraints
- `U`        -- input constraints
"""
ConstrainedBlackBoxControlContinuousSystem

@doc """
    ConstrainedBlackBoxControlDiscreteSystem <: AbstractDiscreteSystem

Discrete-time control system defined by a right-hand side with domain constraints
of the form:

```math
    x_{k+1} = f(x_k, u_k), \\; x_k ∈ \\mathcal{X}, \\;  u_k ∈ \\mathcal{U} \\; \\forall k.
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
- `inputdim` -- number of input variables
- `X`        -- state constraints
- `U`        -- input constraints
"""
ConstrainedBlackBoxControlDiscreteSystem

# ==============
# Noisy systems
# ==============

for (Z, AZ) in ((:NoisyLinearContinuousSystem, :AbstractContinuousSystem),
                (:NoisyLinearDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTA <: AbstractMatrix{T}, MTD <: AbstractMatrix{T}} <: $(AZ)
            A::MTA
            D::MTD
            function $(Z)(A::MTA, D::MTD) where {T, MTA <: AbstractMatrix{T}, MTD <: AbstractMatrix{T}}
                @assert checksquare(A) == size(D,1)
                return new{T, MTA, MTD}(A, D)
            end
        end
        function $(Z)(A::Number, D::Number)
           return $(Z)(hcat(A), hcat(D))
        end

        statedim(s::$Z) = size(s.A,1)
        inputdim(::$Z) = 0
        noisedim(s::$Z) = size(s.D, 2)
        state_matrix(s::$Z) = s.A
        noise_matrix(s::$Z) = s.D
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = true
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = true
            iscontrolled(::$T) = false
            isconstrained(::$T) = false
        end
    end
end

@doc """
    NoisyLinearContinuousSystem

Continuous-time linear system with additive disturbance of the form:

```math
    x(t)' = A x(t) + D w(t) \\; \\forall t.
```

### Fields

- `A` -- state matrix
- `D` -- noise matrix
"""
NoisyLinearContinuousSystem

@doc """
    NoisyLinearDiscreteSystem

Discrete-time linear system with additive disturbance of the form:

```math
    x_{k+1} = A x_k + D w_k \\; \\forall k.
```

### Fields

- `A` -- state matrix
- `D` -- noise matrix
"""
NoisyLinearDiscreteSystem

for (Z, AZ) in ((:NoisyConstrainedLinearContinuousSystem, :AbstractContinuousSystem),
                (:NoisyConstrainedLinearDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTA <: AbstractMatrix{T}, MTD <: AbstractMatrix{T}, ST, WT} <: $(AZ)
            A::MTA
            D::MTD
            X::ST
            W::WT
            function $(Z)(A::MTA, D::MTD, X::ST, W::WT) where {T, MTA <: AbstractMatrix{T}, MTD <: AbstractMatrix{T}, ST, WT}
                @assert checksquare(A) == size(D,1)
                return new{T, MTA, MTD, ST, WT}(A, D, X, W)
            end
        end
        function $(Z)(A::Number, D::Number, X, W)
           return $(Z)(hcat(A), hcat(D), X, W)
        end

        statedim(s::$Z) = size(s.A,1)
        inputdim(::$Z) = 0
        noisedim(s::$Z) = size(s.D, 2)
        stateset(s::$Z) = s.X
        noiseset(s::$Z) = s.W
        state_matrix(s::$Z) = s.A
        noise_matrix(s::$Z) = s.D
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = true
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = true
            iscontrolled(::$T) = false
            isconstrained(::$T) = true
        end
    end
end

@doc """
    NoisyConstrainedLinearContinuousSystem

Continuous-time linear system with  additive disturbance and domain constraints of the form:

```math
    x(t)' = A x(t) + D w(t), \\; x(t) ∈ \\mathcal{X}, \\; w(t) ∈ \\mathcal{W} \\; \\forall t.
```

### Fields

- `A` -- state matrix
- `D` -- noise matrix
- `X` -- state constraints
- `W` -- disturbance set
"""
NoisyConstrainedLinearContinuousSystem

@doc """
    NoisyConstrainedLinearDiscreteSystem

Discrete-time linear system with additive disturbance and domain constraints of the form:

```math
    x_{k+1} = A x_k + D w_k, \\; x_k ∈ \\mathcal{X}, \\; w(t) ∈ \\mathcal{W} \\; \\forall k.
```

### Fields

- `A` -- state matrix
- `D` -- noise matrix
- `X` -- state constraints
- `W` -- disturbance set
"""
NoisyConstrainedLinearDiscreteSystem

for (Z, AZ) in ((:NoisyLinearControlContinuousSystem, :AbstractContinuousSystem),
                (:NoisyLinearControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, MTD <: AbstractMatrix{T}} <: $(AZ)
            A::MTA
            B::MTB
            D::MTD
            function $(Z)(A::MTA, B::MTB, D::MTD) where {T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, MTD <: AbstractMatrix{T}}
                @assert checksquare(A) == size(B, 1) == size(D,1)
                return new{T, MTA, MTB, MTD}(A, B, D)
            end
        end
        function $(Z)(A::Number, B::Number, D::Number)
           return $(Z)(hcat(A), hcat(B), hcat(D))
        end

        statedim(s::$Z) = size(s.A, 1)
        inputdim(s::$Z) = size(s.B, 2)
        noisedim(s::$Z) = size(s.D, 2)
        state_matrix(s::$Z) = s.A
        input_matrix(s::$Z) = s.B
        noise_matrix(s::$Z) = s.D
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = true
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = true
            iscontrolled(::$T) = true
            isconstrained(::$T) = false
        end
    end
end

@doc """
    NoisyLinearControlContinuousSystem

Continuous-time linear control system with additive disturbance of the form:

```math
    x(t)' = A x(t) + B u(t) + D w(t) \\; \\forall t.
```

### Fields

- `A` -- state matrix
- `B` -- input matrix
- `D` -- noise matrix
"""
NoisyLinearControlContinuousSystem

@doc """
    NoisyLinearControlDiscreteSystem

Continuous-time linear control system with additive disturbance of the form:

```math
    x_{k+1} = A x_k + B u_k + D w_k \\; \\forall k.
```

### Fields

- `A` -- state matrix
- `B` -- input matrix
- `D` -- noise matrix
"""
NoisyLinearControlDiscreteSystem

for (Z, AZ) in ((:NoisyConstrainedLinearControlContinuousSystem, :AbstractContinuousSystem),
                (:NoisyConstrainedLinearControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, MTD <: AbstractMatrix{T}, ST, UT, WT} <: $(AZ)
            A::MTA
            B::MTB
            D::MTD
            X::ST
            U::UT
            W::WT
            function $(Z)(A::MTA, B::MTB, D::MTD, X::ST, U::UT, W::WT) where {T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, MTD <: AbstractMatrix{T}, ST, UT, WT}
                @assert checksquare(A) == size(B, 1) == size(D,1)
                return new{T, MTA, MTB, MTD, ST, UT, WT}(A, B, D, X, U, W)
            end
        end
        function $(Z)(A::Number, B::Number, D::Number, X, U, W)
           return $(Z)(hcat(A), hcat(B), hcat(D), X, U, W)
        end

        statedim(s::$Z) = size(s.A, 1)
        inputdim(s::$Z) = size(s.B, 2)
        noisedim(s::$Z) = size(s.D, 2)
        stateset(s::$Z) = s.X
        inputset(s::$Z) = s.U
        noiseset(s::$Z) = s.W
        state_matrix(s::$Z) = s.A
        input_matrix(s::$Z) = s.B
        noise_matrix(s::$Z) = s.D
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = true
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = true
            iscontrolled(::$T) = true
            isconstrained(::$T) = true
        end
    end
end

@doc """
    NoisyConstrainedLinearControlContinuousSystem

Continuous-time linear control system with additive disturbance and domain constraints
of the form:

```math
    x(t)' = A x(t) + B u(t) + D w(t), \\; x(t) ∈ \\mathcal{X}, \\; u(t) ∈ \\mathcal{U}, \\; w(t) ∈ \\mathcal{W} \\; \\forall t.
```

### Fields

- `A` -- state matrix
- `B` -- input matrix
- `D` -- noise matrix
- `X` -- state constraints
- `U` -- input constraints
- `W` -- disturbance set
"""
NoisyConstrainedLinearControlContinuousSystem

@doc """
    NoisyConstrainedLinearControlDiscreteSystem

Continuous-time linear control system with additive disturbance and domain constraints
of the form:

```math
    x_{k+1} = A x_k + B u_k + D w_k, \\; x_k ∈ \\mathcal{X}, \\; u_k ∈ \\mathcal{U}, \\; w_k ∈ \\mathcal{W} \\; \\forall k.
```

### Fields

- `A` -- state matrix
- `B` -- input matrix
- `D` -- noise matrix
- `X` -- state constraints
- `U` -- input constraints
- `W` -- disturbance set
"""
NoisyConstrainedLinearControlDiscreteSystem

for (Z, AZ) in ((:NoisyAffineControlContinuousSystem, :AbstractContinuousSystem),
                (:NoisyAffineControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, VT <: AbstractVector{T}, MTD <: AbstractMatrix{T}} <: $(AZ)
            A::MTA
            B::MTB
            c::VT
            D::MTD
            function $(Z)(A::MTA, B::MTB, c::VT, D::MTD) where {T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, VT <: AbstractVector{T}, MTD <: AbstractMatrix{T}}
                @assert checksquare(A) == length(c) == size(B, 1) == size(D,1)
                return new{T, MTA, MTB, VT, MTD}(A, B, c, D)
            end
        end
        function $(Z)(A::Number, B::Number, c::Number, D::Number)
           return $(Z)(hcat(A), hcat(B), vcat(c), hcat(D))
        end

        statedim(s::$Z) = length(s.c)
        inputdim(s::$Z) = size(s.B, 2)
        noisedim(s::$Z) = size(s.D, 2)
        state_matrix(s::$Z) = s.A
        input_matrix(s::$Z) = s.B
        noise_matrix(s::$Z) = s.D
        affine_term(s::$Z) = s.c
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = false
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = true
            iscontrolled(::$T) = true
            isconstrained(::$T) = false
        end
    end
end

@doc """
    NoisyAffineControlContinuousSystem

Continuous-time affine control system with additive disturbance of the form:

```math
    x(t)' = A x(t) + B u(t) + c + D w(t) \\; \\forall t.
```

### Fields

- `A` -- state matrix
- `B` -- input matrix
- `c` -- affine term
- `D` -- noise matrix
"""
NoisyAffineControlContinuousSystem

@doc """
    NoisyAffineControlDiscreteSystem

Continuous-time affine control system with additive disturbance of the form:

```math
    x_{k+1} = A x_k + B u_k + c + D w_k \\; \\forall k.
```

### Fields

- `A` -- state matrix
- `B` -- input matrix
- `c` -- affine term
- `D` -- noise matrix
"""
NoisyAffineControlDiscreteSystem

for (Z, AZ) in ((:NoisyConstrainedAffineControlContinuousSystem, :AbstractContinuousSystem),
                (:NoisyConstrainedAffineControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, VT <: AbstractVector{T}, MTD <: AbstractMatrix{T}, ST, UT, WT} <: $(AZ)
            A::MTA
            B::MTB
            c::VT
            D::MTD
            X::ST
            U::UT
            W::WT
            function $(Z)(A::MTA, B::MTB, c::VT, D::MTD, X::ST, U::UT, W::WT) where {T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, VT <: AbstractVector{T}, MTD <: AbstractMatrix{T}, ST, UT, WT}
                @assert checksquare(A) == length(c) == size(B, 1) == size(D,1)
                return new{T, MTA, MTB, VT, MTD, ST, UT, WT}(A, B, c, D, X, U, W)
            end
        end
        function $(Z)(A::Number, B::Number, c::Number, D::Number, X, U, W)
           return $(Z)(hcat(A), hcat(B), vcat(c), hcat(D), X, U, W)
        end

        statedim(s::$Z) = length(s.c)
        inputdim(s::$Z) = size(s.B, 2)
        noisedim(s::$Z) = size(s.D, 2)
        stateset(s::$Z) = s.X
        inputset(s::$Z) = s.U
        noiseset(s::$Z) = s.W
        state_matrix(s::$Z) = s.A
        input_matrix(s::$Z) = s.B
        noise_matrix(s::$Z) = s.D
        affine_term(s::$Z) = s.c
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = false
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isblackbox(::$T) = false
            isnoisy(::$T) = true
            iscontrolled(::$T) = true
            isconstrained(::$T) = true
        end
    end
end

@doc """
    NoisyConstrainedAffineControlContinuousSystem

Continuous-time affine control system with additive disturbance and domain constraints
of the form:

```math
    x(t)' = A x(t) + B u(t) + c + D w(t), \\; x(t) ∈ \\mathcal{X}, \\; u(t) ∈ \\mathcal{U}, \\; w(t) ∈ \\mathcal{W} \\; \\forall t.
```

### Fields

- `A` -- state matrix
- `B` -- input matrix
- `c` -- affine term
- `D` -- noise matrix
- `X` -- state constraints
- `U` -- input constraints
- `W` -- disturbance set
"""
NoisyConstrainedAffineControlContinuousSystem

@doc """
    NoisyConstrainedAffineControlDiscreteSystem

Continuous-time affine control system with additive disturbance and domain constraints
of the form:

```math
    x_{k+1} = A x_k + B u_k + c + D w_k, \\; x_k ∈ \\mathcal{X}, \\; u_k ∈ \\mathcal{U}, \\; w_k ∈ \\mathcal{W} \\; \\forall k.
```

### Fields

- `A` -- state matrix
- `B` -- input matrix
- `c` -- affine term
- `D` -- noise matrix
- `X` -- state constraints
- `U` -- input constraints
- `W` -- disturbance set
"""
NoisyConstrainedAffineControlDiscreteSystem

for (Z, AZ) in ((:NoisyBlackBoxControlContinuousSystem, :AbstractContinuousSystem),
                (:NoisyBlackBoxControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){FT} <: $(AZ)
            f::FT
            statedim::Int
            inputdim::Int
            noisedim::Int
        end
        statedim(s::$Z) = s.statedim
        inputdim(s::$Z) = s.inputdim
        noisedim(s::$Z) = s.noisedim
        mapping(s::$Z) = s.f
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = false
            isaffine(::$T) = false
            ispolynomial(::$T) = false
            isblackbox(::$T) = true
            isnoisy(::$T) = true
            iscontrolled(::$T) = true
            isconstrained(::$T) = false
        end
    end
end

@doc """
    NoisyBlackBoxControlContinuousSystem <: AbstractContinuousSystem

Continuous-time disturbance-affected control system defined by a right-hand side
of the form:

```math
    x(t)' = f(x(t), u(t), w(t)) \\; \\forall t.
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
- `inputdim` -- number of input variables
- `noisedim` -- number of noise variables
"""
NoisyBlackBoxControlContinuousSystem

@doc """
    NoisyBlackBoxControlDiscreteSystem <: AbstractDiscreteSystem

Discrete-time disturbance-affected control system defined by a right-hand side
of the form:

```math
    x_{k+1} = f(x_k, u_k) \\quad \\forall k.
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
- `inputdim` -- number of input variables
- `noisedim` -- number of noise variables
"""
NoisyBlackBoxControlDiscreteSystem

for (Z, AZ) in ((:NoisyConstrainedBlackBoxControlContinuousSystem, :AbstractContinuousSystem),
                (:NoisyConstrainedBlackBoxControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){FT, ST, UT,WT} <: $(AZ)
            f::FT
            statedim::Int
            inputdim::Int
            noisedim::Int
            X::ST
            U::UT
            W::WT
        end
        statedim(s::$Z) = s.statedim
        inputdim(s::$Z) = s.inputdim
        noisedim(s::$Z) = s.noisedim
        stateset(s::$Z) = s.X
        inputset(s::$Z) = s.U
        noiseset(s::$Z) = s.W
        mapping(s::$Z) = s.f
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = false
            isaffine(::$T) = false
            ispolynomial(::$T) = false
            isblackbox(::$T) = true
            isnoisy(::$T) = true
            iscontrolled(::$T) = true
            isconstrained(::$T) = true
        end
    end
end

@doc """
    NoisyConstrainedBlackBoxControlContinuousSystem <: AbstractContinuousSystem

Continuous-time disturbance-affected control system defined by a right-hand side
with domain constraints of the form:

```math
    x(t)' = f(x(t), u(t), w(t)), \\quad x(t) ∈ \\mathcal{X}, \\quad u(t) ∈ \\mathcal{U}, \\quad w(t) ∈ \\mathcal{W} \\quad \\forall t.
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
- `inputdim` -- number of input variables
- `noisedim` -- number of noise variables
- `X`        -- state constraints
- `U`        -- input constraints
- `W`        -- disturbance set
"""
NoisyConstrainedBlackBoxControlContinuousSystem

@doc """
    NoisyConstrainedBlackBoxControlDiscreteSystem <: AbstractDiscreteSystem

Discrete-time disturbance-affected control system defined by a right-hand side
with domain constraints of the form:

```math
    x_{k+1} = f(x_k, u_k), \\quad x_k ∈ \\mathcal{X}, \\quad u_k ∈ \\mathcal{U}, \\quad w_k ∈ \\mathcal{W} \\quad \\forall k.
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
- `inputdim` -- number of input variables
- `noisedim` -- number of noise variables
- `X`        -- state constraints
- `U`        -- input constraints
- `W`        -- disturbance set
"""
NoisyConstrainedBlackBoxControlDiscreteSystem

@doc """
    SecondOrderLinearContinuousSystem

Continuous-time second order linear system of the form:

```math
    Mx''(t) + Cx'(t) + Kx(t) = 0\\; \\forall t.
```

### Fields

- `M` -- mass matrix
- `C` -- viscosity matrix
- `K` -- stiffness matrix
"""
SecondOrderLinearContinuousSystem

@doc """
    SecondOrderLinearDiscreteSystem

Discrete-time second order linear system of the form:

```math
    Mx_{k+2} + Cx_{k+1} + Kx_k = 0\\; \\forall k.
```

### Fields

- `M` -- mass matrix
- `C` -- viscosity matrix
- `K` -- stiffness matrix
"""
SecondOrderLinearDiscreteSystem

for (Z, AZ) in ((:SecondOrderLinearContinuousSystem, :AbstractContinuousSystem),
                (:SecondOrderLinearDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTM <: AbstractMatrix{T},
                       MTC <: AbstractMatrix{T},
                       MTK <: AbstractMatrix{T}} <: $(AZ)
            M::MTM
            C::MTC
            K::MTK
            function $(Z)(M::MTM, C::MTC, K::MTK) where {T, MTM <: AbstractMatrix{T},
                                                            MTC <: AbstractMatrix{T},
                                                            MTK <: AbstractMatrix{T}}
                @assert checksquare(M) == checksquare(C) == checksquare(K)
                return new{T, MTM, MTC, MTK}(M, C, K)
            end
        end
        function $(Z)(M::Number, C::Number, K::Number)
             return $(Z)(hcat(M), hcat(C), hcat(K))
        end

        statedim(s::$Z) = size(s.C, 1)
        inputdim(::$Z) = 0
        noisedim(::$Z) = 0
        mass_matrix(s::$Z) = s.M
        viscosity_matrix(s::$Z) = s.C
        stiffness_matrix(s::$Z) = s.K
        affine_term(s::$Z) = zeros(eltype(s.C), size(s.C, 1))
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = true
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = false
            isconstrained(::$T) = false
        end
    end
end

@doc """
    SecondOrderAffineContinuousSystem

Continuous-time second order affine system of the form:

```math
    Mx''(t) + Cx'(t) + Kx(t) = b\\; \\forall t.
```

### Fields

- `M` -- mass matrix
- `C` -- viscosity matrix
- `K` -- stiffness matrix
- `b` -- affine term
"""
SecondOrderAffineContinuousSystem

@doc """
    SecondOrderAffineDiscreteSystem

Discrete-time second order affine system of the form:

```math
    Mx_{k+2} + Cx_{k+1} + Kx_k = b\\; \\forall k.
```

### Fields

- `M` -- mass matrix
- `C` -- viscosity matrix
- `K` -- stiffness matrix
- `b` -- affine term
"""
SecondOrderAffineDiscreteSystem

for (Z, AZ) in ((:SecondOrderAffineContinuousSystem, :AbstractContinuousSystem),
                (:SecondOrderAffineDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTM <: AbstractMatrix{T},
                       MTC <: AbstractMatrix{T},
                       MTK <: AbstractMatrix{T},
                       VT <: AbstractVector{T}} <: $(AZ)
            M::MTM
            C::MTC
            K::MTK
            b::VT
            function $(Z)(M::MTM, C::MTC, K::MTK, b::VT) where {T, MTM <: AbstractMatrix{T},
                                                                   MTC <: AbstractMatrix{T},
                                                                   MTK <: AbstractMatrix{T},
                                                                   VT <: AbstractVector{T}}
                @assert checksquare(M) == checksquare(C) == checksquare(K) == length(b)
                return new{T, MTM, MTC, MTK, VT}(M, C, K, b)
            end
        end
        function $(Z)(M::Number, C::Number, K::Number, b::Number)
             return $(Z)(hcat(M), hcat(C), hcat(K), [b])
        end

        statedim(s::$Z) = size(s.C, 1)
        inputdim(::$Z) = 0
        noisedim(::$Z) = 0
        mass_matrix(s::$Z) = s.M
        viscosity_matrix(s::$Z) = s.C
        stiffness_matrix(s::$Z) = s.K
        affine_term(s::$Z) = s.b
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = true
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = false
            isconstrained(::$T) = false
        end
    end
end

@doc """
    SecondOrderConstrainedLinearControlContinuousSystem

Continuous-time second order constrained linear control system of the form:

```math
    Mx''(t) + Cx'(t) + Kx(t) = Bu(t), \\; x(t) ∈ \\mathcal{X}, \\; u(t) ∈ \\mathcal{U} \\; \\forall t.
```

### Fields

- `M` -- mass matrix
- `C` -- viscosity matrix
- `K` -- stiffness matrix
- `B` -- input matrix
- `X` -- state constraints
- `U` -- input constraints
"""
SecondOrderConstrainedLinearControlContinuousSystem

@doc """
    SecondOrderConstrainedLinearControlDiscreteSystem

Discrete-time second order constrained linear control system of the form:

```math
    Mx_{k+2} + Cx_{k+1} + Kx_k = Bu_k, \\; x_k ∈ \\mathcal{X}, \\; u_k ∈ \\mathcal{U} \\; \\forall k.
```

### Fields

- `M` -- mass matrix
- `C` -- viscosity matrix
- `K` -- stiffness matrix
- `B` -- input matrix
- `X` -- state constraints
- `U` -- input constraints
"""
SecondOrderConstrainedLinearControlDiscreteSystem

for (Z, AZ) in ((:SecondOrderConstrainedLinearControlContinuousSystem, :AbstractContinuousSystem),
                (:SecondOrderConstrainedLinearControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTM <: AbstractMatrix{T},
                       MTC <: AbstractMatrix{T},
                       MTK <: AbstractMatrix{T},
                       MTB <: AbstractMatrix{T},
                       ST,
                       UT} <: $(AZ)
            M::MTM
            C::MTC
            K::MTK
            B::MTB
            X::ST
            U::UT
            function $(Z)(M::MTM, C::MTC, K::MTK, B::MTB, X::ST, U::UT) where {T, MTM <: AbstractMatrix{T},
                                                                                MTC <: AbstractMatrix{T},
                                                                                MTK <: AbstractMatrix{T},
                                                                                MTB <: AbstractMatrix{T},
                                                                                ST,
                                                                                UT}
                @assert checksquare(M) == checksquare(C) == checksquare(K) == size(B, 1)
                return new{T, MTM, MTC, MTK, MTB, ST, UT}(M, C, K, B, X, U)
            end
        end
        function $(Z)(M::Number, C::Number, K::Number, B::Number, X, U)
             return $(Z)(hcat(M), hcat(C), hcat(K), hcat(B), X, U)
        end

        statedim(s::$Z) = size(s.C, 1)
        inputdim(::$Z) = size(s.B, 2)
        noisedim(::$Z) = 0
        mass_matrix(s::$Z) = s.M
        viscosity_matrix(s::$Z) = s.C
        stiffness_matrix(s::$Z) = s.K
        affine_term(s::$Z) = zeros(eltype(s.C), size(s.C, 1))
        input_matrix(s::$Z) = s.B
        stateset(s::$Z) = s.X
        inputset(s::$Z) = s.U
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = true
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = true
            isconstrained(::$T) = true
        end
    end
end

@doc """
    SecondOrderConstrainedAffineControlContinuousSystem

Continuous-time second order constrained affine control system of the form:

```math
    Mx''(t) + Cx'(t) + Kx(t) = Bu(t) + d, \\; x(t) ∈ \\mathcal{X}, \\; u(t) ∈ \\mathcal{U} \\; \\forall t.
```

### Fields

- `M` -- mass matrix
- `C` -- viscosity matrix
- `K` -- stiffness matrix
- `B` -- input matrix
- `d` -- affine term
- `X` -- state constraints
- `U` -- input constraints
"""
SecondOrderConstrainedAffineControlContinuousSystem

@doc """
    SecondOrderConstrainedAffineControlDiscreteSystem

Discrete-time second order constrained affine control system of the form:

```math
    Mx_{k+2} + Cx_{k+1} + Kx_k = Bu_k + d, \\; x_k ∈ \\mathcal{X}, \\; u_k ∈ \\mathcal{U} \\; \\forall k.
```

### Fields

- `M` -- mass matrix
- `C` -- viscosity matrix
- `K` -- stiffness matrix
- `B` -- input matrix
- `d` -- affine term
- `X` -- state constraints
- `U` -- input constraints
"""
SecondOrderConstrainedAffineControlDiscreteSystem

for (Z, AZ) in ((:SecondOrderConstrainedAffineControlContinuousSystem, :AbstractContinuousSystem),
                (:SecondOrderConstrainedAffineControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTM <: AbstractMatrix{T},
                       MTC <: AbstractMatrix{T},
                       MTK <: AbstractMatrix{T},
                       MTB <: AbstractMatrix{T},
                       VT <: AbstractVector{T},
                       ST,
                       UT} <: $(AZ)
            M::MTM
            C::MTC
            K::MTK
            B::MTB
            d::VT
            X::ST
            U::UT
            function $(Z)(M::MTM, C::MTC, K::MTK, B::MTB, d::VT, X::ST, U::UT) where {T, MTM <: AbstractMatrix{T},
                                                                                MTC <: AbstractMatrix{T},
                                                                                MTK <: AbstractMatrix{T},
                                                                                MTB <: AbstractMatrix{T},
                                                                                VT <: AbstractVector{T},
                                                                                ST,
                                                                                UT}
                @assert checksquare(M) == checksquare(C) == checksquare(K) == size(B, 1) == length(d)
                return new{T, MTM, MTC, MTK, MTB, VT, ST, UT}(M, C, K, B, d, X, U)
            end
        end
        function $(Z)(M::Number, C::Number, K::Number, B::Number, d::Number, X, U)
             return $(Z)(hcat(M), hcat(C), hcat(K), hcat(B), [d], X, U)
        end

        statedim(s::$Z) = size(s.C, 1)
        inputdim(::$Z) = size(s.B, 2)
        noisedim(::$Z) = 0
        mass_matrix(s::$Z) = s.M
        viscosity_matrix(s::$Z) = s.C
        stiffness_matrix(s::$Z) = s.K
        affine_term(s::$Z) = s.d
        input_matrix(s::$Z) = s.B
        stateset(s::$Z) = s.X
        inputset(s::$Z) = s.U
    end
    for T in [Z, Type{<:eval(Z)}]
        @eval begin
            islinear(::$T) = true
            isaffine(::$T) = true
            ispolynomial(::$T) = false
            isnoisy(::$T) = false
            iscontrolled(::$T) = true
            isconstrained(::$T) = true
        end
    end
end


text:
struct QBSolver <: Solver
    solver
end

POMDPs.solve(sol::QBSolver, p::Union{POMDP,MDP}) = QBPlanner(solve(sol.solver, p))

struct QBPlanner{P<:DPWPlanner} <: Policy
    p::P
end

srand(p::QBPlanner, s) = srand(p.p, s)

function POMDPModelTools.action_info(qp::QBPlanner, b; tree_in_info=false)
    p = qp.p
    dpw = p
    sol = p.solver

    local a::actiontype(p.mdp)
    info = Dict{Symbol, Any}()
    try
        total_n = 0
        @assert p.solver.enable_action_pw == false
        S = statetype(p.mdp)
        A = actiontype(p.mdp)
        tree = MCTS.DPWTree{S,A}(p.solver.n_iterations)
        p.tree = tree

        banodes = Int[]

        for a in actions(dpw.mdp, b)
            n0 = init_N(sol.init_N, dpw.mdp, b, a)
            anode = insert_action_orphan!(tree, a, n0,
                                          init_Q(sol.init_Q, dpw.mdp, b, a),
                                          false)
            push!(banodes, anode)
            total_n += n0
        end

        i = 0
        start_us = CPUtime_us()
        for i = 1:p.solver.n_iterations
            s = rand(p.rng, b)
            @assert !isterminal(p.mdp, s)

            best_UCB = -Inf
            sanode = 0
            ltn = log(total_n)
            for child in banodes
                n = tree.n[child]
                q = tree.q[child]
                if ltn <= 0 && n == 0
                    UCB = q
                else
                    c = sol.exploration_constant # for clarity
                    UCB = q + c*sqrt(ltn/n)
                end
                @assert !isnan(UCB)
                @assert !isequal(UCB, -Inf)
                if UCB > best_UCB
                    best_UCB = UCB
                    sanode = child
                end
            end

            a = tree.a_labels[sanode]

            new_node = false
            if tree.n_a_children[sanode] <= sol.k_state*tree.n[sanode]^sol.alpha_state
                sp, r = @gen(:sp,:r)(dpw.mdp, s, a, dpw.rng)

                spnode = sol.check_repeat_state ? get(tree.s_lookup, sp, 0) : 0

                if spnode == 0 # there was not a state node for sp already in the tree
                    spnode = MCTS.insert_state_node!(tree, sp, sol.keep_tree || sol.check_repeat_state)
                    new_node = true
                end
                push!(tree.transitions[sanode], (spnode, r))

                if !sol.check_repeat_state
                    tree.n_a_children[sanode] += 1
                elseif !((sanode,spnode) in tree.unique_transitions)
                    push!(tree.unique_transitions, (sanode,spnode))
                    tree.n_a_children[sanode] += 1
                end
            else
                spnode, r = rand(dpw.rng, tree.transitions[sanode])
            end

            if new_node
                q = r + discount(dpw.mdp)*MCTS.estimate_value(dpw.solved_estimate, dpw.mdp, sp, dpw.solver.depth-1)
            else
                q = r + discount(dpw.mdp)*MCTS.simulate(dpw, spnode, dpw.solver.depth-1)
            end

            tree.n[sanode] += 1
            total_n += 1

            tree.q[sanode] += (q - tree.q[sanode])/tree.n[sanode]

            if CPUtime_us() - start_us >= p.solver.max_time * 1e6
                break
            end
        end

        info[:search_time_us] = CPUtime_us() - start_us
        info[:tree_queries] = i
        if p.solver.tree_in_info || tree_in_info
            info[:tree] = tree
        end

        best_Q = -Inf
        sanode = 0
        for child in banodes
            if tree.q[child] > best_Q
                best_Q = tree.q[child]
                sanode = child
            end
        end
        # XXX some publications say to choose action that has been visited the most
        a = tree.a_labels[sanode] # choose action with highest approximate value
    catch ex
        a = convert(actiontype(p.mdp), default_action(p.solver.default_action, p.mdp, b, ex))
        info[:exception] = ex
    end

    return a, info
end

POMDPs.action(p::QBPlanner, b) = first(action_info(p, b))

function insert_action_orphan!(tree::MCTS.DPWTree{S,A}, a::A, n0::Int, q0::Float64, maintain_a_lookup=true) where {S,A}
    push!(tree.n, n0)
    push!(tree.q, q0)
    push!(tree.a_labels, a)
    push!(tree.transitions, Vector{Tuple{Int,Float64}}[])
    sanode = length(tree.n)
    push!(tree.n_a_children, 0)
    @assert !maintain_a_lookup # not supported
    return sanode
end

ext: jl


text:
using Documenter
using GadgetIO

makedocs(
    sitename = "GadgetIO.jl",
    format = Documenter.HTML(),
    modules = [GadgetIO],
    pages = [
            "Table of Contents" => "index.md",
            "Install"           => "install.md",
            "File Infos"        => "file_infos.md",
            "Reading Snapshots" => "read_snapshots.md",
            "Reading Subfind"   => "read_subfind.md",
            "Writing Data"      => "write_data.md",
            "API reference"     => "api.md"
            ]
        )

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/LudwigBoess/GadgetIO.jl.git"
)


text:
using Crystalline
using Crystalline: formatirreplabel
using SymmetryBases
using PhotonicBandConnectivity
using MPBUtils

# --- setup info ---
sgnum  = 147
has_tr = false
id     = 7716
checkfragile = true
calcname = "dim3-sg147-breaktr-detfix-g1.0_symeigs_"*string(id)*"-res32"

# --- hilbert/ebr bases ---
sb, brs = compatibility_basis(sgnum, 3; timereversal=has_tr)
B = matrix(brs)
F = smith(B)

# --- load and process data ---
bandirsd, lgirsd = extract_individual_multiplicities(
                        calcname,
                        timereversal=has_tr,
                        dir = "../../mpb-ctl/output/sg147/",
                        atol=2e-2)
length(lgirsd) ≠ length(sb.klabs) && error("missing k-point data")

bands, nds = collect_separable(bandirsd, lgirsd)
μs = length.(bands)

isempty(bands) && error("   ... found no isolable band candidates ...")

# --- build symmetry vectors ---
# find the permutation between sorting in `sb` and `lgirsd`
permd = Dict(klab => Vector{Int}(undef, length(lgirsd[klab])) for klab ∈ sb.klabs)
for klab in sb.klabs
    lgirs = lgirsd[klab]
    for (i, lgir) in enumerate(lgirs)
        irlab = formatirreplabel(label(lgir))
        j = findfirst(==(irlab), sb.irlabs)
        j === nothing && error("Could not find irrep label $irlab")

        permd[klab][i] = j
    end
end

# use permutation to construct symmetry vectors, in `sb`'s sorting
ns = [Vector{Int}(undef, length(first(sb))) for _ in 1:length(bands)]
for (b, (nd, μ)) in enumerate(zip(nds, μs))
    for (klab, nᵏ) in nd
        permᵏ = permd[klab]
        ns[b][permᵏ] .= nᵏ
    end
    ns[b][end] = μ
end

# --- find which band combinations are in {BS} and then test their topology ---
band′ = 0:0
n′ = similar(first(ns))
is_bs = true # at every new iteration, `is_bs` effectively means `was_prev_iter_bs`
for (band, n) in zip(bands, ns)
    global band′, n′, is_bs

    band′ = is_bs ? band : (minimum(band′):maximum(band))
    n′    = is_bs ? n : n′ + n

    # test if `n′` is in {BS} or not
    is_bs = if first(band′) == 1                        # singular bands
        is_transverse_bandstruct(n′, sb, lgirsd["Γ"], F)
    else                                                # regular bands
        isbandstruct(n′, F)
    end
    is_bs || continue

    # calculate topology of bands
    if first(band′) == 1
        topo = calc_topology_singular(n′, sb, lgirsd["Γ"], B)
    else
        topo = checkfragile ? calc_detailed_topology(n′, B, F) : calc_topology(n′, F)
    end

    println(band′, " ⇒ ", topo)
end


text:
# parser

function read_instance(MyFileName::String)
  path = "./" * MyFileName
  # Si le fichier path existe
  if isfile(path)
    # L’ouvrir
    myFile = open(path)
    readline(myFile)

    n = parse(Int64, split(readline(myFile), " ")[3])
    P = parse(Int64, split(readline(myFile), " ")[3])

    cout = Array{Int64,2}(ones(n,n)) # on met a zero pour l'instant
    for i in 1:n
      for j in 1:n
        line =  split(readline(myFile), " ")
        cout[parse(Int64,line[1]),parse(Int64,line[2])] = parse(Int64,line[3])
      end
    end

    alpha1 = Array{Float64,1}(zeros(P))
    for i in 1:P
      alpha1[i] = 0.5 # donnees du poly
    end
    alpha2 = Array{Float64,1}(zeros(P))
    for i in 1:3
      alpha2[i] = 0.9 # donnees du poly
    end
    for i in 4:6
      alpha2[i] = 0.5 # donnees du poly
    end
    alpha3 = Array{Float64,1}(zeros(P))
    for i in 1:3
      alpha3[i] = 0.5 # donnees du poly
    end
    for i in 4:6
      alpha3[i] = 0.9 # donnees du poly
    end
    alpha4 = Array{Float64,1}(zeros(P))
    for i in 1:3
      alpha4[i] = 0.8 # donnees du poly
    end
    for i in 4:6
      alpha4[i] = 0.6 # donnees du poly
    end

    indice_danger = Array{Float64,1}(zeros(P))
    for i in 1:3
      indice_danger[i] = 1
    end

    Survie = Array{Float64,3}(zeros(P,n,n))
    data = readlines(myFile)
		for datum in data
      line = split(datum, " ")
      Survie[parse(Int64,line[1]),parse(Int64,line[2]),parse(Int64,line[3])] = parse(Float64,line[4])
    end
    return n, P, alpha4, Survie, indice_danger, cout
  end
end

max_stars_repo_name: boumpteryx/RODD

max_issues_repo_name: boumpteryx/RODD

max_forks_repo_name: boumpteryx/RODD


text:
struct HestonModel{T}
    v0::T
    κ::T
    θ::T
    ρ::T
    σ::T
    r::T
    q::T
end


text:
using Caesar, RoME
using TransformUtils, Rotations, CoordinateTransformations
using PyCall, PyLCM
using LibBSON
# using Colors, FixedPointNumbers #needed if pushing to DrakeVisualizer

# generate python bindings for lcmtypes and add them to PYTHONPATH
@show lcmtpath = joinpath(dirname(@__FILE__),"lcmtypes")
run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/hauv_submap_points_t.lcm`)
println("Adding lcmtypes dir to Python path: $(lcmtpath)")
unshift!(PyVector(pyimport("sys")["path"]),lcmtpath)
@pyimport hauv

function lcmpointcloudmsg!(vc, slaml::SLAMWrapper,
                          msgdata  )
  msg = hauv.submap_points_t[:decode](msgdata)

  vert = getVert(slaml.fg, slaml.lastposesym, api=IncrementalInference.dlapi) # fetch from database

  # 2d arrays of points and colors (from LCM data into arrays{arrays})
  ptarr = [[pt[1], pt[2], pt[3]] for pt in msg[:points_local]]
  carr = [[UInt8(c.data[1]),UInt8(c.data[2]),UInt8(c.data[3])] for c in msg[:colors]]

  # push to mongo (using BSON as a quick fix)
  # (for deserialization, see src/DirectorVisService.jl:cachepointclouds!)
  serialized_point_cloud = BSONObject(Dict("pointcloud" => ptarr))
  appendvertbigdata!(slaml.fg, vert, "BSONpointcloud", string(serialized_point_cloud).data)
  serialized_colors = BSONObject(Dict("colors" => carr))
  appendvertbigdata!(slaml.fg, vert, "BSONcolors", string(serialized_colors).data)

  # push to DrakeVisualizer
  #=
  position = Translation(msg[:pos]...)
  orientation = Rotations.Quat(msg[:orientation]...)
  Tf = position ∘ LinearMap(orientation)
  vsym = Symbol(vert.label)
  setgeometry!(vc[:submaps][vsym], Triad())
  settransform!(vc[:submaps][vsym], Tf)

  pointcloud = PointCloud(ptarr)
  pointcloud.channels[:rgb] = map(x -> RGB(reinterpret(N0f8, x)...), carr)
  setgeometry!(vc[:submaps][vsym][:pc], pointcloud)
  =#

  nothing
end

# this function is on notice!
function setupSLAMinDB(;cloudGraph=nothing, addrdict=nothing)
  if cloudGraph != nothing
    return SLAMWrapper(initfg(sessionname=addrdict["session"], cloudgraph=cloudGraph), nothing, 0)
  else
    return SLAMWrapper(RoME.initfg(), nothing, 0)
  end
end

# prepare the factor graph with just one node
# (will prompt on stdin for db credentials)
cloudGraph, addrdict = standardcloudgraphsetup()

slam = setupSLAMinDB(cloudGraph=cloudGraph, addrdict=addrdict)
slam.fg = identitypose6fg(fg=slam.fg)

function listenone(slaml)
  lc = LCM()

  # lcm callback closure
  pchdl = (channel, msg_data) -> lcmpointcloudmsg!( slaml, msg_data )

  subscribe(lc, "SUBMAPS", pchdl)

  # run this thing for one message only....
  handle(lc)
end

listenone(slam)

# Juno.breakpoint(@__FILE__(), 52)


text:
using Statistics, DataFrames, CSV
data = CSV.read("data/purchaseData.csv", DataFrame)

println(mean(data.Price),"\n")
println(mean(skipmissing(data.Price)),"\n")
println(coalesce.(data.Grade, "QQ")[1:4],"\n")
println(first(dropmissing(data,:Price), 4),"\n")
println(sum(ismissing.(data.Name)),"\n")
println(findall(completecases(data))[1:4])

max_stars_repo_name: newptcai/StatsWithJuliaBook

max_issues_repo_name: newptcai/StatsWithJuliaBook

max_forks_repo_name: newptcai/StatsWithJuliaBook


text:
Base.Enums.basetype(::Type{CE}) where {CE<:Cenum} = Base.Enums.basetype(bitstype(CE))

Base.Enums.namemap(ce::Cenum) = Base.Enums.namemap(typeof(ce))
function Base.Enums.namemap(::Type{CE}) where {CE<:Cenum}
	result = ()
	for (n, v) in values(bitstype(CE))
		result = (result..., (n, v))
	end
	return result
end

(::Type{CE})(val::Integer = 0) where {CE<:Cenum} = Core.Intrinsics.bitcast(bitstype(CE), convert(Base.Enums.basetype(CE), val))::CE
(::Type{T})(ce::CE) where {T<:Integer, CE<:Cenum} = convert(T, Core.Intrinsics.bitcast(Base.Enums.basetype(CE), ce))::T

Base.sizeof(::Type{CE}) where {CE<:Cenum} = invoke(sizeof, Tuple{Any}, bitstype(CE))
Base.zero(::Type{CE}) where {CE<:Cenum} = CE()

Base.convert(::Type{CE}, val::Integer) where {CE<:Cenum} = CE(val)
Base.convert(::Type{T}, ce::CE) where {CE<:Cenum, T<:Integer} = T(ce)

Base.promote_rule(::Type{T}, ::Type{CE}) where {T<:Integer, CE<:Cenum} = promote_type(T, Base.Enums.basetype(CE))
Base.:(~)(ce::CE) where {CE<:Cenum} = ~convert(Base.Enums.basetype(CE), ce)
for op in (:+, :-, :&, :|, :xor, :(==))
	@eval begin
		Base.$(op)(x::CE, y::CE) where {CE<:Cenum} = ($op)(Core.Intrinsics.bitcast(Base.Enums.basetype(CE), x), Core.Intrinsics.bitcast(Base.Enums.basetype(CE), y))
		Base.$(op)(x::CE, y::T) where {CE<:Cenum, T<:Integer} = ($op)(promote_type(T, CE)(x), promote_type(T, CE)(y))
		Base.$(op)(x::T, y::CE) where {T<:Integer, CE<:Cenum} = ($op)(promote_type(T, CE)(x), promote_type(T, CE)(y))
	end
end

function Base.typemin(::Type{CE}) where {CE<:Cenum}
	x = nothing
	for (n, v) in Base.Enums.namemap(CE)
		x = x === nothing || v < x ? v : x
	end
	return x === nothing ? CE() : CE(x)
end

function Base.typemax(::Type{CE}) where {CE<:Cenum}
	x = nothing
	for (n, v) in Base.Enums.namemap(CE)
		x = x === nothing || v >= x ? v : x
	end
	return x === nothing ? CE() : CE(x)
end

function Base.instances(::Type{CE}) where {CE<:Cenum}
	result = ()
	for (n, v) in Base.Enums.namemap(CE)
		result = (result..., CE(v))
	end
	return result
end

function Base.Symbol(ce::Cenum)
	for (n, v) in Base.Enums.namemap(ce)
		typeof(v)(ce) == v && return n::Symbol
	end
	error("Value of $(typeof(ce)) is not one of $(join(map(first, Base.Enums.namemap(ce)), ", "))")
end

function Base.show(io::IO, ce::Cenum)
	for (ind, val) in enumerate(sort(collect(instances(typeof(ce))), by = val -> convert(Base.Enums.basetype(typeof(ce)), val)))
		if ce == val
			print(io, Symbol(ce))
			return
		end
	end

	show(io, typeof(ce))
	print(io, " (")
	show(io, convert(Base.Enums.basetype(typeof(ce)), ce))
	println(io, ')')
end

function Base.show(io::IO, m::MIME"text/plain", ce::Cenum)
	show(io, typeof(ce))
	print(io, " (")
	show(io, convert(Base.Enums.basetype(typeof(ce)), ce))
	println(io, ')')

	numvals = length(instances(typeof(ce)))
	maxwidth = mapreduce(val -> length(String(Symbol(val))), max, instances(typeof(ce)), init = 0)

	indent = get(io, :indent, 1)
	prefix = repeat("  ", indent-1)
	for (ind, val) in enumerate(sort(collect(instances(typeof(ce))), by = val -> convert(Base.Enums.basetype(typeof(ce)), val)))
		select = ce == val ? "->" : "  "
		print(io, prefix, select, rpad(String(Symbol(val)), maxwidth), " (")
		show(io, convert(Base.Enums.basetype(typeof(ce)), val))
		print(io, ')')
		ind < numvals && println(io)
	end
end

function Base.show(io::IO, m::MIME"text/plain", ce::Type{<:Cenum})
	show(io, ce)

	numvals = length(instances(ce))
	println(io, " (", numvals, " option", numvals > 1 ? "s)" : ")")

	maxwidth = mapreduce(val -> length(String(Symbol(val))), max, instances(ce), init = 0)

	indent = get(io, :indent, 1)
	prefix = repeat("  ", indent)
	for (ind, val) in enumerate(sort(collect(instances(ce)), by = val -> convert(Base.Enums.basetype(ce), val)))
		print(io, prefix, rpad(String(Symbol(val)), maxwidth), " (")
		show(io, convert(Base.Enums.basetype(ce), val))
		print(io, ')')
		ind < numvals && println(io)
	end
end


text:
## COMSOL POSITION TYPES ==================================================== ##
abstract type CPosition end

struct CPosition2Dsym <: CPosition
    r::Vector{Float64}
    z::Vector{Float64}
end

## COMSOL EXPORTED DATA TYPES =============================================== ##
abstract type CData end

struct CDataRaw <: CData
    path::String                # file path
    meta::Dict{String, Any}     # meta data (from header from file)
    field::Vector{String}       # names of the exported fields
    v::Dict{String, Any}        # values of fields
end

struct CDataTable <: CData
    path::String                # file path
    meta::Dict{String, Any}     # meta data (from header from file)
    n::Vector{String}           # names
    u::Vector{String}           # units
    v::Dict{String, Any}        # values
end

struct CDataDSet <: CData
    path::String                # file path
    meta::Dict{String, Any}     # meta data (from header from file)
    n::Vector{String}           # names (without time and position)
    u::Vector{String}           # units (without time and position)
    t::Vector{Float64}          # time [s]
    p::CPosition                # position [m]
    v::Dict{String, Any}        # values
end


text:
# ## Heatmap

#md # ![FILE_NAME.png](images/FILE_NAME.png)

using Random
Random.seed!(2018)

a = rand(5,5)
xlabel = string.(collect('A':'E'))
ylabel = string.(collect('a':'e'))
heatmap(a, xticks=(1:5, xlabel), yticks=(1:5, ylabel),
    fill_z=a, aspect_ratio=:equal)

fontsize = 15
nrow, ncol = size(a)
ann = [(i,j, text(round(a[i,j], digits=2), fontsize, :white, :center))
            for i in 1:nrow for j in 1:ncol]
annotate!(ann, linecolor=:white)


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

@doc raw"""DeleteOptions may be provided when deleting an API object.

    IoK8sApimachineryPkgApisMetaV1DeleteOptions(;
        apiVersion=nothing,
        dryRun=nothing,
        gracePeriodSeconds=nothing,
        kind=nothing,
        orphanDependents=nothing,
        preconditions=nothing,
        propagationPolicy=nothing,
    )

    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - dryRun::Vector{String} : When present, indicates that modifications should not be persisted. An invalid or unrecognized dryRun directive will result in an error response and no further processing of the request. Valid values are: - All: all dry run stages will be processed
    - gracePeriodSeconds::Int64 : The duration in seconds before the object should be deleted. Value must be non-negative integer. The value zero indicates delete immediately. If this value is nil, the default grace period for the specified type will be used. Defaults to a per object value if not specified. zero means delete immediately.
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - orphanDependents::Bool : Deprecated: please use the PropagationPolicy, this field will be deprecated in 1.7. Should the dependent objects be orphaned. If true/false, the \"orphan\" finalizer will be added to/removed from the object's finalizers list. Either this field or PropagationPolicy may be set, but not both.
    - preconditions::IoK8sApimachineryPkgApisMetaV1Preconditions : Must be fulfilled before a deletion is carried out. If not possible, a 409 Conflict status will be returned.
    - propagationPolicy::String : Whether and how garbage collection will be performed. Either this field or OrphanDependents may be set, but not both. The default policy is decided by the existing finalizer set in the metadata.finalizers and the resource-specific default policy. Acceptable values are: 'Orphan' - orphan the dependents; 'Background' - allow the garbage collector to delete the dependents in the background; 'Foreground' - a cascading policy that deletes all dependents in the foreground.
"""
mutable struct IoK8sApimachineryPkgApisMetaV1DeleteOptions <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    dryRun::Any # spec type: Union{ Nothing, Vector{String} } # spec name: dryRun
    gracePeriodSeconds::Any # spec type: Union{ Nothing, Int64 } # spec name: gracePeriodSeconds
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    orphanDependents::Any # spec type: Union{ Nothing, Bool } # spec name: orphanDependents
    preconditions::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Preconditions } # spec name: preconditions
    propagationPolicy::Any # spec type: Union{ Nothing, String } # spec name: propagationPolicy

    function IoK8sApimachineryPkgApisMetaV1DeleteOptions(;apiVersion=nothing, dryRun=nothing, gracePeriodSeconds=nothing, kind=nothing, orphanDependents=nothing, preconditions=nothing, propagationPolicy=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1DeleteOptions, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1DeleteOptions, Symbol("dryRun"), dryRun)
        setfield!(o, Symbol("dryRun"), dryRun)
        validate_property(IoK8sApimachineryPkgApisMetaV1DeleteOptions, Symbol("gracePeriodSeconds"), gracePeriodSeconds)
        setfield!(o, Symbol("gracePeriodSeconds"), gracePeriodSeconds)
        validate_property(IoK8sApimachineryPkgApisMetaV1DeleteOptions, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApimachineryPkgApisMetaV1DeleteOptions, Symbol("orphanDependents"), orphanDependents)
        setfield!(o, Symbol("orphanDependents"), orphanDependents)
        validate_property(IoK8sApimachineryPkgApisMetaV1DeleteOptions, Symbol("preconditions"), preconditions)
        setfield!(o, Symbol("preconditions"), preconditions)
        validate_property(IoK8sApimachineryPkgApisMetaV1DeleteOptions, Symbol("propagationPolicy"), propagationPolicy)
        setfield!(o, Symbol("propagationPolicy"), propagationPolicy)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1DeleteOptions

const _property_map_IoK8sApimachineryPkgApisMetaV1DeleteOptions = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("dryRun")=>Symbol("dryRun"), Symbol("gracePeriodSeconds")=>Symbol("gracePeriodSeconds"), Symbol("kind")=>Symbol("kind"), Symbol("orphanDependents")=>Symbol("orphanDependents"), Symbol("preconditions")=>Symbol("preconditions"), Symbol("propagationPolicy")=>Symbol("propagationPolicy"))
const _property_types_IoK8sApimachineryPkgApisMetaV1DeleteOptions = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("dryRun")=>"Vector{String}", Symbol("gracePeriodSeconds")=>"Int64", Symbol("kind")=>"String", Symbol("orphanDependents")=>"Bool", Symbol("preconditions")=>"IoK8sApimachineryPkgApisMetaV1Preconditions", Symbol("propagationPolicy")=>"String")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1DeleteOptions }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1DeleteOptions))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1DeleteOptions }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1DeleteOptions[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1DeleteOptions }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1DeleteOptions[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1DeleteOptions)
    true
end

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


text:
# # Solve BTPDE
#
# We start by loading SpinDoctor and a Makie plotting backend.

# LSP indexing solution                                                          #src
# https://github.com/julia-vscode/julia-vscode/issues/800#issuecomment-650085983 #src
if isdefined(@__MODULE__, :LanguageServer)                                       #src
    include("../src/SpinDoctor.jl")                                              #src
    using .SpinDoctor                                                            #src
end                                                                              #src

using SpinDoctor
using LinearAlgebra

if haskey(ENV, "GITHUB_ACTIONS")
    using CairoMakie
else
    using GLMakie
end

# The built in geometry recipes allow for making various cell configuration. We here consider
# the case of three twisted axons immersed in an extracellular space (ECS).

setup = CylinderSetup(;
    name = "some-very-real-axons",
    ncell = 3,
    rmin = 2.0,
    rmax = 6.0,
    dmin = 0.2,
    dmax = 0.3,
    height = 40.0,
    bend = 0.0,
    twist = π / 4,
    include_in = false,
    in_ratio = 0.6,
    ecs_shape = :convex_hull,
    ecs_ratio = 0.5,
)

# We also define coefficients for the different cell compartments `:in` (axon), `:out`
# (myelin), and `:ecs` (ECS).

coeffs = coefficients(
    setup;
    D = (; in = 0.002 * I(3), out = 0.002 * I(3), ecs = 0.002 * I(3)),
    T₂ = (; in = Inf, out = Inf, ecs = Inf),
    ρ = (; in = 1.0, out = 1.0, ecs = 1.0),
    κ = (; in_out = 1e-4, out_ecs = 1e-4, in = 0.0, out = 0.0, ecs = 0.0),
    γ = 2.67513e-4,
)

# The following line creates a random cell configuration for our cylinder setup, generates a
# surface triangulation and calls TetGen to create a tetrahedral finite element mesh. The
# compartments and boundaries will be ordered in the same way as `coeffs`.

mesh, = create_geometry(setup)

# The resulting mesh can be plotted in 3D provided a Makie backend is loaded.

plot_mesh(mesh)

# The mesh looks good, so we can proceed with the assembly our biological model and the
# associated finite element matrices.

model = Model(; mesh, coeffs...)
matrices = assemble_matrices(model)

# The Bloch-Torrey PDE takes a magnetic field gradient pulse sequence as an input. Here
# we consider a [`ScalarGradient`](@ref) with a [`PGSE`](@ref) time profile.

dir = [1.0, 0.0, 0.0]
profile = PGSE(2000.0, 6000.0)
b = 1000
g = √(b / int_F²(profile)) / coeffs.γ
gradient = ScalarGradient(dir, profile, g)

# SpinDoctor provides a [`solve`](@ref) function, which has the same base signature for all
# diffusion MRI problems. The BTPDE is one such problem. They generally take a gradient
# sequence as an input.

btpde = BTPDE(; model, matrices)
ξ = solve(btpde, gradient)

# Here, `ξ` is a vector containing the complex-valued magnetization at all degrees of freedom
# at the echo time `TE`. We may compute the resulting signal as follows:

compute_signal(matrices.M, ξ)

# The global mass matrix `M` is used to compute the integral. We may however be interested in
# the compartment-wise signals. This requires splitting the magnetization field into the
# respective compartments. The compartment mass matrices are also available.

ξ_cmpts = split_field(mesh, ξ)
compute_signal.(matrices.M_cmpts, ξ_cmpts)

# The final magnetization can be visualized using the [`plot_field`](@ref) function.

plot_field(mesh, ξ)

# In this example, we have computed the complex transverse water proton magnetization field
# using the finite element method. The measured diffusion MRI signal is the integral of this
# field, and other quantities of interest, such as the apparent diffusion coefficient (ADC),
# or the effective diffusion tensor, may easily be obtained from this reference field.
# Directly solving the BTPDE is thus considered to be the "gold standard" for computing these
# quantities, as arbitrary precision may be obtained.

# However, this is also often the most computationally expensive approach. In the following
# examples, we will consider some other specialized methods provided by SpinDoctor, each
# having their own domains of validity, use cases, and computational footprints.


text:
# Bitpair operations
# ==================
#
# Copyright Ben J. Ward and Luis Yanes.
# License is MIT.

@inline function count_nonzero_bitpairs(x::UInt64)
    out = x & 0x5555555555555555
    out |= (x & 0xAAAAAAAAAAAAAAAA) >> 1
    return count_ones(out)
end

@inline function count_zero_bitpairs(x::UInt64)
    return 32 - count_nonzero_bitpairs(x)
end

@inline function count_one_bitpairs(x::UInt64)
    out = x & 0x5555555555555555
    out &= (x & 0xAAAAAAAAAAAAAAAA) >> 1
    return count_ones(out)
end


text:
module NeuralNets

for (dir, filename) in [
	# CORE
	("","constants.jl"),
	("","types.jl"),
	# NEURON
	("neurons","simple_spiking.jl"),
	("neurons","hodgkin_huxley.jl"),
	# Synapses
	("synapses","spiking_synapse.jl"),
	("synapses","spiking_conn.jl"),
	("synapses","stdp.jl"),
	# Layers
	("layers","spiking_layer.jl"),
	# Network
	("nets","spiking_layers_net.jl"),
	("nets","small_net.jl"),
	# simulation
	("","run_results.jl"),
	# APMA2821V-specific codes
	("apma2821v","connectivities.jl"),
	("apma2821v","buildnet.jl"),
	("apma2821v","stimuli.jl")
	#("apma2821v","plotting.jl")
	]
	include(joinpath(dir, filename))
end

# NEURON
export HHNeuron, SimpleSpikingNeuron
# SYNAPSE
export SpikingSynapse, SynapseType, EXC, INH
export SpikingConnection, stdp
# LAYERS
export SpikingLayer, update!
# I/O
export RunSettings, RunResults
# NETWORKS
export SpikingLayersNet, add_layer!, connect!, runSpikingNet!, save, loadNet, resetLayers!, normalizeWeights!
# APMA2821V
export all_to_all, p_lgn_connectivity, forward_connectivity, lateral_connectivity
export buildnet
export image_stim
#export plot_all_rasters, plot_raster, plot_stim

end # module NeuralNets


text:
arg_dict = Dict{Any, Any}(
"physics" => "Euler",
"run_type" => 5,
"jac_type" => 2,
"jac_method" => 2,
"newton_globalize_euler" => true,
"order" => 1,
"real_time" => false,
"variable_type" => :conservative,
"smb_name" => "SRCMESHES/2D_airfoil.smb",
"dmg_name" => "SRCMESHES/2D_airfoil.dmg", # "src/mesh_files/2D_airfoil.dmg",
"IC_name" => "ICFreeStream",
"numBC" => 2,
"BC1" => [8],
"BC1_name" => "FreeStreamBC",
"BC2" => [5],
"BC2_name" => "noPenetrationBC",
"aoa" => 2.0,
"Ma" => 0.5,# default working = 0.3
"use_DG" => true,
"Flux_name" => "RoeFlux",
"t_max" => 500000.0,
"delta_t" => 0.0001,
"res_tol" => 1e-10,
"res_abstol" => 1e-11,
"res_reltol" => 1e-11,
"step_tol" => 1e-10,
"itermax" => 5000000000, # For RK4
"writeq" => false,
"write_edge_vertnums" => false,
"write_face_vertnums" => false,
"write_qic" => false,
"writeboundary" => false,
"write_res" => false,
"print_cond" => false,
"write_counts" => false,
"write_vis" => false,
"output_freq" => 100000,
"solve" => true,
"do_postproc" => true,
"force_solution_complex" => true,
"need_adjoint" => true,

# functional options
"functional_name1" => "lift",
"functional_bcs1" => [2],

"exact_soln_func" => "ICFreeStream",
)


text:
@testset "Data" begin
    @test artificialIn_SAR_image(2) == 2 * π * ones(2, 2)

    @test artificial_S1_slope_signal(20, 0.0) == repeat([-π / 2], 20)

    @test ismissing(artificial_S1_signal(-1.0))
    @test ismissing(artificial_S1_signal(2.0))
    @test artificial_S1_signal(2) == [-3 * π / 4, -3 * π / 4]

    # for the remainder check data types only
    @test length(artificial_S1_signal(20)) == 20

    @test size(artificial_S2_whirl_image(64)) == (64, 64)
    @test length(artificial_S2_whirl_image(64)[1, 1]) == 3

    @test size(artificial_S2_rotation_image(64)) == (64, 64)
    @test length(artificial_S2_rotation_image(64)[1, 1]) == 3

    @test size(artificial_S2_whirl_patch(8)) == (8, 8)
    @test length(artificial_S2_whirl_patch(8)[1, 1]) == 3

    @test size(artificial_SPD_image(8)) == (8, 8)
    @test size(artificial_SPD_image(8)[1, 1]) == (3, 3)

    @test size(artificial_SPD_image2(8)) == (8, 8)
    @test size(artificial_SPD_image2(8)[1, 1]) == (3, 3)
    @test eltype(artificial_SPD_image2(8)) == Array{Float64,2}

    @test length(artificial_S2_lemniscate([0.0, 0.0, 1.0], 20)) == 20
    @test length(artificial_S2_lemniscate([0.0, 0.0, 1.0], 20)[1]) == 3
end


text:
"""
    abstract type Rewriter

Pluggable extension to a [`Project`] with hooks to transform
individual documents, create new documents, register file update
handlers and perform additional build steps.

See the following methods:
- [`rewritedoc`](#)
- [`createdocs`](#)
- [`reset!`](#)
- [`postbuild`](#)
- [`getfilehandlers`](#)
"""
abstract type Rewriter end

"""
    rewritedoc(rewriter, p, doc) -> doc'

Rewrite `doc` at path `p`. Return rewritten `doc`.
"""
function rewritedoc(rewriter::Rewriter, p, doc)
    return doc
end

"""
    reset!(rewriter)

Clears internal state of `rewriter`. Does nothing if not overwritten.
"""
function reset!(rewriter::Rewriter) end

"""
    postbuild(rewriter, project, dst, format)

Post-build callback for [`Rewriter`](#)s.
"""
function postbuild(rewriter, project, builder) end

createsources!(::Rewriter) = Dict{AbstractPath, XTree}()

struct Replacer <: Rewriter
    fn
    selector::Selector
end

Base.show(io::IO, replacer::Replacer) = print(io, "Replacer($(replacer.selector))")

function rewritedoc(replace::Replacer, p::AbstractPath, doc::XTree)
    return cata(replace.fn, doc, replace.selector)
end


text:
module EMLiae
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...types.EMLtypes
	import EMLtypes.length
	include("iae/IAE.jl")
end


text:
using NQCDynamics.Calculators: DiabaticFrictionCalculator

abstract type AbstractMDEF <: DynamicsMethods.Method end

"""
```math
dr = v dt\\\\
dv = -\\Delta U/M dt - \\Gamma v dt + \\sigma \\sqrt{2\\Gamma} dW
```
``\\Gamma`` is the friction tensor with units of inverse time.
For thermal dynamics we set ``\\sigma = \\sqrt{kT / M}``,
where ``T`` is the electronic temperature.

This is integrated using the BAOAB algorithm where the friction "O" step is performed
in the tensor's eigenbasis. See `src/dynamics/mdef_baoab.jl` for details.
"""
struct MDEF{M} <: AbstractMDEF
    mass_scaling::M
end

MDEF(masses::AbstractVector, DoFs::Integer) = MDEF(get_mass_scale_matrix(masses, DoFs))

function get_mass_scale_matrix(masses::AbstractVector, DoFs::Integer)
    vect = repeat(sqrt.(masses), inner=DoFs)
    vect * vect'
end

function NQCDynamics.Simulation{MDEF}(atoms::Atoms, model::Model; kwargs...)
    NQCDynamics.Simulation(atoms, model, MDEF(atoms.masses, ndofs(model));  kwargs...)
end
function NQCDynamics.RingPolymerSimulation{MDEF}(atoms::Atoms, model::Model, n_beads::Integer; kwargs...)
    NQCDynamics.RingPolymerSimulation(atoms, model, MDEF(atoms.masses, ndofs(model)), n_beads; kwargs...)
end

"""
    acceleration!(dv, v, r, sim::Simulation{MDEF,<:DiabaticFrictionCalculator}, t)

Sets acceleration due to ground state force when using a `DiabaticFrictionModel`.
"""
function acceleration!(dv, v, r, sim::Simulation{MDEF,<:DiabaticFrictionCalculator}, t)
    Calculators.update_electronics!(sim.calculator, r)
    for I in eachindex(dv)
        dv[I] = -sim.calculator.adiabatic_derivative[I][1,1]
    end
    DynamicsUtils.divide_by_mass!(dv, sim.atoms.masse)
end

"""
    friction!(g, r, sim, t)

Evaluates friction tensor
"""
function friction!(g, r, sim::AbstractSimulation{<:AbstractMDEF}, t)
    Calculators.evaluate_friction!(sim.calculator, r)
    g .= sim.calculator.friction ./ sim.method.mass_scaling
end

function DynamicsMethods.create_problem(u0, tspan::Tuple, sim::AbstractSimulation{<:AbstractMDEF})
    StochasticDiffEq.DynamicalSDEProblem(acceleration!, DynamicsUtils.velocity!, friction!,
        DynamicsUtils.get_velocities(u0), DynamicsUtils.get_positions(u0), tspan, sim)
end


text:
using Cairo
using LSystems
using DataStructures

c = CairoRGBSurface(1024, 896);
cr = CairoContext(c);

save(cr);

set_source_rgb(cr, 1.0, 1.0, 1.0);
rectangle(cr, 0.0, 0.0, 1024.0, 896.0);
fill(cr);

restore(cr);
save(cr);

x, y, a, d = 64., 896., -pi / 6, 0.8
set_line_width(cr, 0.6);
move_to(cr, x, y);

plant_start = "X"
plant_trans = Dict([('X', "F-[[X]+X]+F[+FX]-X"), ('F', "FF")])

function determine_new_position(x, y, a, d)
    x -= sin(a) * d
    y -= cos(a) * d
    return x, y
end

s = Stack(Vector{Float64}) # [x, y, a]

level = 9
drawing_instructions = @task lindenmayer(plant_start, plant_trans, level)

for v in drawing_instructions
    if v == 'F'
        x, y = determine_new_position(x, y, a, d)
        line_to(cr, x, y)
    elseif v == '+'
        a -= pi / 7
    elseif v == '-'
        a += pi / 7
    elseif v == '['
        push!(s, [x, y, a])
    elseif v == ']'
        x, y, a = pop!(s)
        move_to(cr, x, y)
    end
end

stroke(cr);
write_to_png(c, "plant_$level.png");