text:
# This file contains ODESystem prototypes

import DifferentialEquations: Tsit5, ODEProblem
import UUIDs: uuid4

"""
    @def_ode_system ex

where `ex` is the expression to define to define a new AbstractODESystem component type. The usage is as follows:
```julia
@def_ode_system mutable struct MyODESystem{T1,T2,T3,...,TN,OP,RH,RO,ST,IP,OP} <: AbstractODESystem
    param1::T1 = param1_default                     # optional field
    param2::T2 = param2_default                     # optional field
    param3::T3 = param3_default                     # optional field
        ⋮
    paramN::TN = paramN_default                     # optional field
    righthandside::RH = righthandeside_function     # mandatory field
    readout::RO = readout_function                  # mandatory field
    state::ST = state_default                       # mandatory field
    input::IP = input_default                       # mandatory field
    output::OP = output_default                     # mandatory field
end
```
Here, `MyODESystem` has `N` parameters. `MyODESystem` is represented by the `righthandside` and `readout` function. `state`, `input` and `output` is the state, input port and output port of `MyODESystem`.

!!! warning
    `righthandside` must have the signature
    ```julia
    function righthandside(dx, x, u, t, args...; kwargs...)
        dx .= .... # update dx
    end
    ```
    and `readout` must have the signature
    ```julia
    function readout(x, u, t)
        y = ...
        return y
    end
    ```

!!! warning
    New ODE system must be a subtype of `AbstractODESystem` to function properly.

!!! warning
    New ODE system must be mutable type.

# Example
```jldoctest
julia> @def_ode_system mutable struct MyODESystem{RH, RO, IP, OP} <: AbstractODESystem
       α::Float64 = 1.
       β::Float64 = 2.
       righthandside::RH = (dx, x, u, t, α=α) -> (dx[1] = α * x[1] + u[1](t))
       readout::RO = (x, u, t) -> x
       state::Vector{Float64} = [1.]
       input::IP = Inport(1)
       output::OP = Outport(1)
       end

julia> ds = MyODESystem();

julia> ds.input
1-element Inport{Inpin{Float64}}:
 Inpin(eltype:Float64, isbound:false)
```
"""
macro def_ode_system(ex)
    checksyntax(ex, :AbstractODESystem)
    fields = quote
        trigger::$(TRIGGER_TYPE_SYMBOL) = Inpin()
        handshake::$(HANDSHAKE_TYPE_SYMBOL) = Outpin{Bool}()
        callbacks::$(CALLBACKS_TYPE_SYMBOL) = nothing
        name::Symbol = Symbol()
        id::$(ID_TYPE_SYMBOL) = Causal.uuid4()
        t::Float64 = 0.
        modelargs::$(MODEL_ARGS_TYPE_SYMBOL) = ()
        modelkwargs::$(MODEL_KWARGS_TYPE_SYMBOL) = NamedTuple()
        solverargs::$(SOLVER_ARGS_TYPE_SYMBOL) = ()
        solverkwargs::$(SOLVER_KWARGS_TYPE_SYMBOL) = NamedTuple()
        alg::$(ALG_TYPE_SYMBOL) = Causal.Tsit5()
        integrator::$(INTEGRATOR_TYPE_SYMBOL) = Causal.construct_integrator(Causal.ODEProblem, input, righthandside, state, t, modelargs,
            solverargs; alg=alg, modelkwargs=modelkwargs, solverkwargs=solverkwargs, numtaps=3)
    end, [TRIGGER_TYPE_SYMBOL, HANDSHAKE_TYPE_SYMBOL, CALLBACKS_TYPE_SYMBOL, ID_TYPE_SYMBOL, MODEL_ARGS_TYPE_SYMBOL, MODEL_KWARGS_TYPE_SYMBOL, SOLVER_ARGS_TYPE_SYMBOL, SOLVER_KWARGS_TYPE_SYMBOL, ALG_TYPE_SYMBOL, INTEGRATOR_TYPE_SYMBOL]
    _append_common_fields!(ex, fields...)
    deftype(ex)
end

#### Define ODE system library.

"""
    ODESystem(;righthandside, readout, state, input, output)

Constructs a generic ODE system.

# Example
```jldoctest
julia> ds = ODESystem(righthandside=(dx,x,u,t)->(dx.=-x), readout=(x,u,t)->x, state=[1.],input=nothing, output=Outport(1));

julia> ds.state
1-element Array{Float64,1}:
 1.0
```
"""
@def_ode_system mutable struct ODESystem{RH, RO, ST, IP, OP} <: AbstractODESystem
    righthandside::RH
    readout::RO
    state::ST
    input::IP
    output::OP
end

@doc raw"""
    ContinuousLinearSystem(input, output, modelargs=(), solverargs=();
        A=fill(-1, 1, 1), B=fill(0, 1, 1), C=fill(1, 1, 1), D=fill(0, 1, 1), state=rand(size(A,1)), t=0.,
        alg=ODEAlg, modelkwargs=NamedTuple(), solverkwargs=NamedTuple())

Constructs a `ContinuousLinearSystem` with `input` and `output`. `state` is the initial state and `t` is the time. `modelargs` and `modelkwargs` are passed into `ODEProblem` and `solverargs` and `solverkwargs` are passed into `solve` method of `DifferentialEquations`. `alg` is the algorithm to solve the differential equation of the system.

The `ContinuousLinearSystem` is represented by the following state and output equations.
```math
\begin{array}{l}
    \dot{x} = A x + B u \\[0.25cm]
    y = C x + D u
\end{array}
```
where ``x`` is `state`. `solver` is used to solve the above differential equation.
"""
@def_ode_system mutable struct ContinuousLinearSystem{IP, OP, RH, RO} <: AbstractODESystem
    A::Matrix{Float64} = fill(-1., 1, 1)
    B::Matrix{Float64} = fill(1., 1, 1)
    C::Matrix{Float64} = fill(1., 1, 1)
    D::Matrix{Float64} = fill(0., 1, 1)
    input::IP = Inport(1)
    output::OP = Outport(1)
    state::Vector{Float64} = rand(size(A, 1))
    righthandside::RH = input === nothing ? (dx, x, u, t) -> (dx .= A * x) :
        (dx, x, u, t) -> (dx .= A * x + B * map(ui -> ui(t), u.itp))
    readout::RO = input === nothing ? (x, u, t) -> (C * x) :
           ( (C === nothing || D === nothing) ? nothing : (x, u, t) -> (C * x + D * map(ui -> ui(t), u)) )
end

@doc raw"""
    LorenzSystem(input, output, modelargs=(), solverargs=();
        sigma=10, beta=8/3, rho=28, gamma=1, outputfunc=allstates, state=rand(3), t=0.,
        alg=ODEAlg, cplmat=diagm([1., 1., 1.]), modelkwargs=NamedTuple(), solverkwargs=NamedTuple())

Constructs a `LorenzSystem` with `input` and `output`. `sigma`, `beta`, `rho` and `gamma` is the system parameters. `state` is the initial state and `t` is the time. `modelargs` and `modelkwargs` are passed into `ODEProblem` and `solverargs` and `solverkwargs` are passed into `solve` method of `DifferentialEquations`. `alg` is the algorithm to solve the differential equation of the system.

If `input` is `nothing`, the state equation of `LorenzSystem` is
```math
\begin{array}{l}
    \dot{x}_1 = \gamma (\sigma (x_2 - x_1)) \\[0.25cm]
    \dot{x}_2 = \gamma (x_1 (\rho - x_3) - x_2) \\[0.25cm]
    \dot{x}_3 = \gamma (x_1 x_2 - \beta x_3)
\end{array}
```
where ``x`` is `state`. `solver` is used to solve the above differential equation. If `input` is not `nothing`, then the state eqaution is
```math
\begin{array}{l}
    \dot{x}_1 = \gamma (\sigma (x_2 - x_1)) + \sum_{j = 1}^3 \alpha_{1j} u_j \\[0.25cm]
    \dot{x}_2 = \gamma (x_1 (\rho - x_3) - x_2) + \sum_{j = 1}^3 \alpha_{2j} u_j \\[0.25cm]
    \dot{x}_3 = \gamma (x_1 x_2 - \beta x_3) + \sum_{j = 1}^3 \alpha_{3j} u_j
\end{array}
```
where ``A = [\alpha_{ij}]`` is `cplmat` and ``u = [u_{j}]`` is the value of the `input`. The output function is
```math
    y = g(x, u, t)
```
where ``t`` is time `t`, ``y`` is the value of the `output` and ``g`` is `outputfunc`.
"""
@def_ode_system mutable struct LorenzSystem{RH, RO, IP, OP} <: AbstractODESystem
    σ::Float64 = 10.
    β::Float64 = 8 / 3
    ρ::Float64 = 28.
    γ::Float64 = 1.
    righthandside::RH = function lorenzrhs(dx, x, u, t, σ=σ, β=β, ρ=ρ, γ=γ)
        dx[1] = σ * (x[2] - x[1])
        dx[2] = x[1] * (ρ - x[3]) - x[2]
        dx[3] = x[1] * x[2] - β * x[3]
        dx .*= γ
    end
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(3)
    input::IP = nothing
    output::OP = Outport(3)
end

"""
    ForcedLorenzSystem()

Constructs a LorenzSystem that is driven by its inputs.
"""
@def_ode_system mutable struct ForcedLorenzSystem{CM, RH, RO, IP, OP} <: AbstractODESystem
    σ::Float64 = 10.
    β::Float64 = 8 / 3
    ρ::Float64 = 28.
    γ::Float64 = 1.
    cplmat::CM = I(3)
    righthandside::RH = function lorenzrhs(dx, x, u, t, σ=σ, β=β, ρ=ρ, γ=γ, cplmat=cplmat)
        dx[1] = σ * (x[2] - x[1])
        dx[2] = x[1] * (ρ - x[3]) - x[2]
        dx[3] = x[1] * x[2] - β * x[3]
        dx .*= γ
        dx .+= cplmat * map(ui -> ui(t), u.itp)   # Couple inputs
    end
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(3)
    input::IP = Inport(3)
    output::OP = Outport(3)
end

@doc raw"""
    ChenSystem(input, output, modelargs=(), solverargs=();
        a=35, b=3, c=28, gamma=1, outputfunc=allstates, state=rand(3), t=0.,
        alg=ODEAlg, cplmat=diagm([1., 1., 1.]), modelkwargs=NamedTuple(), solverkwargs=NamedTuple())

Constructs a `ChenSystem` with `input` and `output`. `a`, `b`, `c` and `gamma` is the system parameters. `state` is the initial state and `t` is the time. `modelargs` and `modelkwargs` are passed into `ODEProblem` and `solverargs` and `solverkwargs` are passed into `solve` method of `DifferentialEquations`. `alg` is the algorithm to solve the differential equation of the system.

If `input` is `nothing`, the state equation of `ChenSystem` is
```math
\begin{array}{l}
    \dot{x}_1 = \gamma (a (x_2 - x_1)) \\[0.25cm]
    \dot{x}_2 = \gamma ((c - a) x_1 + c x_2 + x_1 x_3) \\[0.25cm]
    \dot{x}_3 = \gamma (x_1 x_2 - b x_3)
\end{array}
```
where ``x`` is `state`. `solver` is used to solve the above differential equation. If `input` is not `nothing`, then the state eqaution is
```math
\begin{array}{l}
    \dot{x}_1 = \gamma (a (x_2 - x_1)) + \sum_{j = 1}^3 \alpha_{1j} u_j \\[0.25cm]
    \dot{x}_2 = \gamma ((c - a) x_1 + c x_2 + x_1 x_3) + \sum_{j = 1}^3 \alpha_{2j} u_j \\[0.25cm]
    \dot{x}_3 = \gamma (x_1 x_2 - b x_3) + \sum_{j = 1}^3 \alpha_{3j} u_j
\end{array}
```
where ``A = [\alpha_{ij}]`` is `cplmat` and ``u = [u_{j}]`` is the value of the `input`. The output function is
```math
    y = g(x, u, t)
```
where ``t`` is time `t`, ``y`` is the value of the `output` and ``g`` is `outputfunc`.
"""
@def_ode_system mutable struct ChenSystem{RH, RO, IP, OP} <: AbstractODESystem
    a::Float64 = 35.
    b::Float64 = 3.
    c::Float64 = 28.
    γ::Float64 = 1.
    righthandside::RH = function chenrhs(dx, x, u, t, a=a, b=b, c=c, γ=γ)
        dx[1] = a * (x[2] - x[1])
        dx[2] = (c - a) * x[1] + c * x[2] - x[1] * x[3]
        dx[3] = x[1] * x[2] - b * x[3]
        dx .*= γ
    end
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(3)
    input::IP = nothing
    output::OP = Outport(3)
end

"""
    ForcedChenSystem()

Constructs Chen system driven by its inputs.
"""
@def_ode_system mutable struct ForcedChenSystem{CM, RH, RO, IP, OP} <: AbstractODESystem
    a::Float64 = 35.
    b::Float64 = 3.
    c::Float64 = 28.
    γ::Float64 = 1.
    cplmat::CM = I(3)
    righthandside::RH = function chenrhs(dx, x, u, t, a=a, b=b, c=c, γ=γ)
        dx[1] = a * (x[2] - x[1])
        dx[2] = (c - a) * x[1] + c * x[2] - x[1] * x[3]
        dx[3] = x[1] * x[2] - b * x[3]
        dx .*= γ
        dx .+= cplmat * map(ui -> ui(t), u.itp)   # Couple inputs
    end
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(3)
    input::IP = Inport(3)
    output::OP = Outport(3)
end

Base.@kwdef struct PiecewiseLinearDiode
    m0::Float64 = -1.143
    m1::Float64 = -0.714
    m2::Float64 = 5.
    bp1::Float64 = 1.
    bp2::Float64 = 5.
end
@inline function (d::PiecewiseLinearDiode)(x)
    m0, m1, m2, bp1, bp2 = d.m0, d.m1, d.m2, d.bp1, d.bp2
    if x < -bp2
        return m2 * x + (m2 - m1) * bp2 + (m1 - m0) * bp1
    elseif -bp2 <= x < -bp1
        return m1 * x + (m1 - m0) * bp1
    elseif -bp1 <= x < bp1
        return m0 * x
    elseif bp1 <= x < bp2
        return m1 * x +  (m0 - m1) * bp1
    else
        return m2 * x + (m1 - m2) * bp2 + (m0 - m1) * bp1
    end
end

Base.@kwdef struct PolynomialDiode
    a::Float64 = 1 / 16
    b::Float64 = - 1 / 6
end
(d::PolynomialDiode)(x) = d.a * x^3 + d.b * x

@doc raw"""
    ChuaSystem(input, output, modelargs=(), solverargs=();
        diode=PiecewiseLinearDiode(), alpha=15.6, beta=28., gamma=1., outputfunc=allstates, state=rand(3), t=0.,
        alg=ODEAlg, cplmat=diagm([1., 1., 1.]), modelkwargs=NamedTuple(), solverkwargs=NamedTuple())

Constructs a `ChuaSystem` with `input` and `output`. `diode`, `alpha`, `beta` and `gamma` is the system parameters. `state` is the initial state and `t` is the time. `modelargs` and `modelkwargs` are passed into `ODEProblem` and `solverargs` and `solverkwargs` are passed into `solve` method of `DifferentialEquations`. `alg` is the algorithm to solve the differential equation of the system.

If `input` is `nothing`, the state equation of `ChuaSystem` is
```math
\begin{array}{l}
    \dot{x}_1 = \gamma (\alpha (x_2 - x_1 - h(x_1))) \\[0.25cm]
    \dot{x}_2 = \gamma (x_1 - x_2 + x_3 ) \\[0.25cm]
    \dot{x}_3 = \gamma (-\beta x_2)
\end{array}
```
where ``x`` is `state`. `solver` is used to solve the above differential equation. If `input` is not `nothing`, then the state eqaution is
```math
\begin{array}{l}
    \dot{x}_1 = \gamma (\alpha (x_2 - x_1 - h(x_1))) + \sum_{j = 1}^3 \theta_{1j} u_j \\[0.25cm]
    \dot{x}_2 = \gamma (x_1 - x_2 + x_3 ) + \sum_{j = 1}^3 \theta_{2j} u_j \\[0.25cm]
    \dot{x}_3 = \gamma (-\beta x_2) + \sum_{j = 1}^3 \theta_{3j} u_j
\end{array}
```
where ``\Theta = [\theta_{ij}]`` is `cplmat` and ``u = [u_{j}]`` is the value of the `input`. The output function is
```math
    y = g(x, u, t)
```
where ``t`` is time `t`, ``y`` is the value of the `output` and ``g`` is `outputfunc`.
"""
@def_ode_system mutable struct ChuaSystem{DT,RH, RO, IP, OP} <: AbstractODESystem
    diode::DT = PiecewiseLinearDiode()
    α::Float64 = 15.6
    β::Float64 = 28.
    γ::Float64 = 1.
    righthandside::RH = function chuarhs(dx, x, u, t, diode=diode, α=α, β=β, γ=γ)
        dx[1] = α * (x[2] - x[1] - diode(x[1]))
        dx[2] = x[1] - x[2] + x[3]
        dx[3] = -β * x[2]
        dx .*= γ
    end
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(3)
    input::IP = nothing
    output::OP = Outport(3)
end

"""
    ForcedChuaSystem()

Constructs a Chua system with inputs.
"""
@def_ode_system mutable struct ForcedChuaSystem{DT, CM, RH, RO, IP, OP} <: AbstractODESystem
    diode::DT = PiecewiseLinearDiode()
    α::Float64 = 15.6
    β::Float64 = 28.
    γ::Float64 = 1.
    cplmat::CM = I(3)
    righthandside::RH = function chuarhs(dx, x, u, t, diode=diode, α=α, β=β, γ=γ)
        dx[1] = α * (x[2] - x[1] - diode(x[1]))
        dx[2] = x[1] - x[2] + x[3]
        dx[3] = -β * x[2]
        dx .*= γ
        dx .+= cplmat * map(ui -> ui(t), u.itp)   # Couple inputs
    end
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(3)
    input::IP = Inport(3)
    output::OP = Outport(3)
end

@doc raw"""
    RosslerSystem(input, output, modelargs=(), solverargs=();
        a=0.38, b=0.3, c=4.82, gamma=1., outputfunc=allstates, state=rand(3), t=0.,
        alg=ODEAlg, cplmat=diagm([1., 1., 1.]), modelkwargs=NamedTuple(), solverkwargs=NamedTuple())

Constructs a `RosllerSystem` with `input` and `output`. `a`, `b`, `c` and `gamma` is the system parameters. `state` is the initial state and `t` is the time. `modelargs` and `modelkwargs` are passed into `ODEProblem` and `solverargs` and `solverkwargs` are passed into `solve` method of `DifferentialEquations`. `alg` is the algorithm to solve the differential equation of the system.

If `input` is `nothing`, the state equation of `RosslerSystem` is
```math
\begin{array}{l}
    \dot{x}_1 = \gamma (-x_2 - x_3) \\[0.25cm]
    \dot{x}_2 = \gamma (x_1 + a x_2) \\[0.25cm]
    \dot{x}_3 = \gamma (b + x_3 (x_1 - c))
\end{array}
```
where ``x`` is `state`. `solver` is used to solve the above differential equation. If `input` is not `nothing`, then the state eqaution is
```math
\begin{array}{l}
    \dot{x}_1 = \gamma (-x_2 - x_3) + \sum_{j = 1}^3 \theta_{1j} u_j \\[0.25cm]
    \dot{x}_2 = \gamma (x_1 + a x_2 ) + \sum_{j = 1}^3 \theta_{2j} u_j \\[0.25cm]
    \dot{x}_3 = \gamma (b + x_3 (x_1 - c)) + \sum_{j = 1}^3 \theta_{3j} u_j
\end{array}
```
where ``\Theta = [\theta_{ij}]`` is `cplmat` and ``u = [u_{j}]`` is the value of the `input`. The output function is
```math
    y = g(x, u, t)
```
where ``t`` is time `t`, ``y`` is the value of the `output` and ``g`` is `outputfunc`.
"""
@def_ode_system mutable struct RosslerSystem{RH, RO, IP, OP} <: AbstractODESystem
    a::Float64 = 0.38
    b::Float64 = 0.3
    c::Float64 = 4.82
    γ::Float64 = 1.
    righthandside::RH = function rosslerrhs(dx, x, u, t, a=a, b=b, c=c, γ=γ)
        dx[1] = -x[2] - x[3]
        dx[2] = x[1] + a * x[2]
        dx[3] = b + x[3] * (x[1] - c)
        dx .*= γ
    end
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(3)
    input::IP = nothing
    output::OP = Outport(3)
end

"""
    ForcedRosslerSystem()

Constructs a Rossler system driven by its input.
"""
@def_ode_system mutable struct ForcedRosslerSystem{CM, RH, RO, IP, OP} <: AbstractODESystem
    a::Float64 = 0.38
    b::Float64 = 0.3
    c::Float64 = 4.82
    γ::Float64 = 1.
    cplmat::CM = I(3)
    righthandside::RH = function rosslerrhs(dx, x, u, t, a=a, b=b, c=c, γ=γ)
        dx[1] = -x[2] - x[3]
        dx[2] = x[1] + a * x[2]
        dx[3] = b + x[3] * (x[1] - c)
        dx .*= γ
        dx .+= cplmat * map(ui -> ui(t), u.itp)   # Couple inputs
    end
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(3)
    input::IP = Inport(3)
    output::OP = Outport(3)
end

@doc raw"""
    VanderpolSystem(input, output, modelargs=(), solverargs=();
        mu=5., gamma=1., outputfunc=allstates, state=rand(2), t=0.,
        alg=ODEAlg, cplmat=diagm([1., 1]), modelkwargs=NamedTuple(), solverkwargs=NamedTuple())

Constructs a `VanderpolSystem` with `input` and `output`. `mu` and `gamma` is the system parameters. `state` is the initial state and `t` is the time. `modelargs` and `modelkwargs` are passed into `ODEProblem` and `solverargs` and `solverkwargs` are passed into `solve` method of `DifferentialEquations`. `alg` is the algorithm to solve the differential equation of the system.

If `input` is `nothing`, the state equation of `VanderpolSystem` is
```math
\begin{array}{l}
    \dot{x}_1 = \gamma (x_2) \\[0.25cm]
    \dot{x}_2 = \gamma (\mu (x_1^2 - 1) x_2 - x_1 )
\end{array}
```
where ``x`` is `state`. `solver` is used to solve the above differential equation. If `input` is not `nothing`, then the state eqaution is
```math
\begin{array}{l}
    \dot{x}_1 = \gamma (x_2) + \sum_{j = 1}^3 \theta_{1j} u_j \\[0.25cm]
    \dot{x}_2 = \gamma (\mu (x_1^2 - 1) x_2 - x_1) + \sum_{j = 1}^3 \theta_{2j} u_j
\end{array}
```
where ``\Theta = [\theta_{ij}]`` is `cplmat` and ``u = [u_{j}]`` is the value of the `input`. The output function is
```math
    y = g(x, u, t)
```
where ``t`` is time `t`, ``y`` is the value of the `output` and ``g`` is `outputfunc`.
"""
@def_ode_system mutable struct VanderpolSystem{RH, RO, IP, OP} <: AbstractODESystem
    μ::Float64 = 5.
    γ::Float64 = 1.
    righthandside::RH = function vanderpolrhs(dx, x, u, t, μ=μ, γ=γ)
        dx[1] = x[2]
        dx[2] = -μ * (x[1]^2 - 1) * x[2] - x[1]
        dx .*= γ
    end
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(2)
    input::IP = nothing
    output::OP = Outport(3)
end

"""
    ForcedVanderpolSystem()

Constructs a Vanderpol system driven by its input.
"""
@def_ode_system mutable struct ForcedVanderpolSystem{CM, RH, RO, IP, OP} <: AbstractODESystem
    μ::Float64 = 5.
    γ::Float64 = 1.
    cplmat::CM = I(2)
    righthandside::RH = function vanderpolrhs(dx, x, u, t, μ=μ, γ=γ)
        dx[1] = x[2]
        dx[2] = -μ * (x[1]^2 - 1) * x[2] - x[1]
        dx .*= γ
        dx .+= cplmat * map(ui -> ui(t), u.itp)   # Couple inputs
    end
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(2)
    input::IP = Inport(2)
    output::OP = Outport(3)
end

@doc raw"""
    Integrator(state=zeros(0), t=0., modelargs=(), solverargs=();
        alg=ODEAlg, modelkwargs=NamedTuple(), solverkwargs=NamedTuple(), numtaps=numtaps, callbacks=nothing,
        name=Symbol())

Constructs an integrator whose input output relation is given by
```math
u(t) = ki * \int_{0}^{t} u(\tau) d\tau
```
where ``u(t)`` is the input, ``y(t)`` is the output and ``ki`` is the integration constant.
"""
@def_ode_system mutable struct Integrator{RH, RO, IP, OP} <: AbstractODESystem
    ki::Float64 = 1.
    righthandside::RH = (dx, x, u, t) -> (dx[1] = ki * u[1](t))
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(1)
    input::IP = Inport()
    output::OP = Outport()
end

##### Pretty-printing
show(io::IO, ds::ODESystem) = print(io,
    "ODESystem(righthandside:$(ds.righthandside), readout:$(ds.readout), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::ContinuousLinearSystem) = print(io,
    "ContinuousLinearSystem(A:$(ds.A), B:$(ds.B), C:$(ds.C), D:$(ds.D), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::LorenzSystem) = print(io,
    "LorenzSystem(σ:$(ds.σ), β:$(ds.β), ρ:$(ds.ρ), γ:$(ds.γ), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::ForcedLorenzSystem) = print(io,
    "LorenzSystem(σ:$(ds.σ), β:$(ds.β), ρ:$(ds.ρ), γ:$(ds.γ), cplmat:$(ds.cplmat), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::ChenSystem) = print(io,
    "ChenSystem(a:$(ds.a), b:$(ds.b), c:$(ds.c), γ:$(ds.γ), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::ForcedChenSystem) = print(io,
    "ChenSystem(a:$(ds.a), b:$(ds.b), c:$(ds.c), γ:$(ds.γ), cplmat:$(ds.cplmat), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::ChuaSystem) = print(io,
    "ChuaSystem(diode:$(ds.diode), α:$(ds.α), β:$(ds.β), γ:$(ds.γ), state:$(ds.state), ",
    "t:$(ds.t), input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::ForcedChuaSystem) = print(io,
    "ChuaSystem(diode:$(ds.diode), α:$(ds.α), β:$(ds.β), γ:$(ds.γ), cplmat:$(ds.cplmat), state:$(ds.state), ",
    "t:$(ds.t), input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::RosslerSystem) = print(io,
    "RosslerSystem(a:$(ds.a), b:$(ds.b), c:$(ds.c), γ:$(ds.γ), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::ForcedRosslerSystem) = print(io,
    "RosslerSystem(a:$(ds.a), b:$(ds.b), c:$(ds.c), γ:$(ds.γ), cplmat:$(ds.cplmat), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::VanderpolSystem) = print(io,
    "VanderpolSystem(μ:$(ds.μ), γ:$(ds.γ), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::ForcedVanderpolSystem) = print(io,
    "VanderpolSystem(μ:$(ds.μ), γ:$(ds.γ), cplmat:$(ds.cplmat), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::Integrator) = print(io,
    "Integrator(ki:$(ds.ki), state:$(ds.state), t:$(ds.t), input:$(ds.input), output:$(ds.output))")


text:
#' ---
#' title: Problems on Graphs
#' ---

#' **Originally Contributed by**: Arpit Bhatia

#' In the mathematical discipline of graph theory, a number of problems can be solved by modelling them as optimization problems.
#' We will see some examples of such problems in the tutorial.
#' These problems are also sometimes referred to as combinatorial optimization problems as
#' they consist of finding an optimal object from a finite set of objects.

#' Let's first import all the packages we will be using for this tutorial.
#' We'll be plotting all the graphs we work with in this tutorial for which we will need some additional packages.

using JuMP
using GLPK
using GraphPlot
using LightGraphs
using Colors

#' ## Representing Graphs
#' For the purpose of this tutorial, we will represent graphs using adjacency matrices.
#' An adjacency matrix, sometimes also called the connection matrix, is a square matrix used to represent a finite graph.
#' Its rows and columns are labeled by the graph vertices,
#' with a 1 or 0 in position ($v_{i}$,$v_{j}$) according to whether $v_{i}$ and $v_{j}$ are adjacent or not.

#' ## Minimum Vertex Cover
#' Given a graph $G = (V, E)$, a vertex-cover $V' \subset V$ of $G$ is a collection of vertices such that
#' each edge in $E$ is incident to at least one of the vertices in $V'$.
#' The size of a vertex-cover $|V'|$ is the number of vertices present in the cover.
#' We wish to find the minimum vertex cover of $G$ i.e. a minimum size vertex cover.
#' We model this problem as an ILP by defining a decision variable $y_{v}$ for each vertex $v \in V$ and
#' a constraint for each edge $e \in E$ as follows:

#' $$
#' \begin{align*}
#' \min && \sum_{v \in V} y_{v} \\
#' s.t. && y_{u} + y_{v} \geq 1 && \forall \{u,v\} \in E \\
#' && y_{v} \in \{0,1\} && \forall v \in V
#' \end{align*}
#' $$

G = [
0 1 0 0 0 0;
1 0 1 1 0 0;
0 1 0 0 1 1;
0 1 0 0 1 0;
0 0 1 1 0 0;
0 0 1 0 0 0
]

g = SimpleGraph(G)

gplot(g)

#+

vertex_cover = Model(with_optimizer(GLPK.Optimizer))

@variable(vertex_cover, y[1:nv(g)], Bin)
@constraint(vertex_cover, [i = 1:nv(g), j = 1:nv(g); G[i,j] == 1], y[i] + y[j] >= 1)
@objective(vertex_cover, Min, sum(y))

optimize!(vertex_cover)
@show value.(y);

#+

membership = convert(Array{Int},value.(y)) # Change to Int
membership = membership + ones(Int, nv(g)) # Make the color groups one indexed
nodecolor = [colorant"red", colorant"blue"] # Blue to represent vertices in the cover
nodefillc = nodecolor[membership]
gplot(g, nodefillc = nodefillc)

#' ## Dominating Set
#' A dominating set in a graph $G = (V, E)$ is a set $S \subset V$ such that
#' for each vertex $v \in V$ either $v$ or one of its neighbour should be in $S$.
#' Note that for some vertex $u$, $u$ and its neighbour both can be present in $S$.
#' We wish to find the smallest dominating set for a graph.
#' We model this problem as an ILP by defining a decision variable $x_{v}$ for each vertex $v \in V$ along with
#' a constraint for its neighbourhood.

#' $$
#' \begin{align*}
#' \min && \sum_{v \in V} x_{v} \\
#' s.t. && \sum_{u \in N(v)}x_{u} \geq 1 && \forall v \in V \\
#' && x_{v} \in \{0,1\} && \forall v \in V
#' \end{align*}
#' $$

G = [
0 1 0 0 0 0 0 0 0 1 0 ;
1 0 1 0 0 0 0 0 0 0 1;
0 1 0 1 0 1 0 0 0 0 0;
0 0 1 0 1 0 0 0 0 0 0;
0 0 0 1 0 1 0 0 0 0 0;
0 0 1 0 1 0 1 0 0 0 0;
0 0 0 0 0 1 0 1 0 0 0;
0 0 0 0 0 0 1 0 1 0 1;
0 0 0 0 0 0 0 1 0 1 1;
1 0 0 0 0 0 0 0 1 0 1;
0 1 0 0 0 0 0 1 1 1 0
]

g = SimpleGraph(G)

gplot(g)

#+

dominating_set = Model(with_optimizer(GLPK.Optimizer))

@variable(dominating_set, x[1:nv(g)], Bin)
@constraint(dominating_set, [i = 1:nv(g)], sum(G[i,:] .* x) >= 1)
@objective(dominating_set, Min, sum(x))

optimize!(dominating_set)
@show value.(x);

#+

membership = convert(Array{Int},value.(x)) # Change to Int
membership = membership + ones(Int, nv(g)) # Make the color groups one indexed
nodecolor = [colorant"red", colorant"blue"] # Blue to represent vertices in the set
nodefillc = nodecolor[membership]
gplot(g, nodefillc = nodefillc)

#' ## Maximum Matching Problem
#' Given a graph $G = (V, E)$, a matching $M \subset E$ of $G$ is a collection of vertex disjoint edges.
#' The size of the matching $M$ is the number of edges present in $M$ i.e. $|M|$.
#' We wish to find the Maximum matching of $G$ i.e. a matching of maximum size.
#' We can solve this problem by modelling it as an integer linear program (ILP).
#' We define a decision variable $m_{e}$ for each edge $e \in E$ and a constraint for each vertex $u \in V$ as follows:

#' $$
#' \begin{align*}
#' \max && \sum_{e \in E} m_{e} \\
#' s.t. && \sum_{e \sim u} m_{e} \leq 1 && \forall u \in V \\
#' && m_{e} \in \{0,1\} && \forall e \in E
#' \end{align*}
#' $$

#' Let's now use JuMP to solve this problem for a sample graph.

G = [
0 0 0 0 1 0 0 0;
0 0 0 0 0 1 0 0;
0 0 0 0 0 0 1 0;
0 0 0 0 0 0 0 1;
1 0 0 0 0 1 0 1;
0 1 0 0 1 0 1 0;
0 0 1 0 0 1 0 1;
0 0 0 1 1 0 1 0;
]

g = SimpleGraph(G)

gplot(g)

#+

matching = Model(with_optimizer(GLPK.Optimizer))

@variable(matching, m[i = 1:nv(g), j = 1:nv(g)], Bin)
@constraint(matching, [i = 1:nv(g)], sum(m[i,:]) <= 1)
@constraint(matching, [i = 1:nv(g), j = 1:nv(g); G[i,j] == 0], m[i,j] == 0)
@constraint(matching, [i = 1:nv(g), j = 1:nv(g)], m[i,j] == m[j,i])
@objective(matching, Max, sum(m))

optimize!(matching)
@show value.(m);

#' The edges corresponding to the Matching are marked as one in the above matrix.

#' ## k-Coloring Problem
#' A k-coloring of a graph $G=(V,E)$ is a function $c: V \rightarrow \{1,2...k\}$ such that
#' $c(u) \neq c(v)$ for every edge $(u,v) \in E$. In other words, the numbers 1,2...k represent k colors,
#' and adjacent vertices must have different colours.
#' The goal of a graph coloring problem is to find a minimum number of colours needed to colour a graph.

#' We model this problem as an ILP by defining a variable decision variable $z_{i}$ for each colour we have available.
#' Given an upper bound $k$ on the number of colors needed,
#' we use $|V| \times k$ decision variables $c_{v,k}$ denoting if vertex $v$ is assigned color $k$.
#' Our model will become:

#' $$
#' \begin{align*}
#' \min && \sum_{i=1}^{k} z_{i} \\
#' s.t. && \sum_{i=1}^{k} c_{v,i} = 1 && \forall v \in V \\
#' && c_{u,i} + c_{v,i} \leq 1 && \forall (u,v) \in V, i \in \{1,2...k\} \\
#' && c_{v,i} \in \{0,1\} && \forall v \in V, i \in \{1,2...k\} \\
#' && z_{i} \in \{0,1\} && \forall i \in \{1,2...k\}
#' \end{align*}
#' $$

G = [
0 1 0 0 1 1 0 0 0 0;
1 0 1 0 0 0 1 0 0 0;
0 1 0 1 0 0 0 1 0 0;
0 0 1 0 1 0 0 0 1 0;
1 0 0 1 0 0 0 0 0 1;
1 0 0 0 0 0 1 0 0 1;
0 1 0 0 0 1 0 1 0 0;
0 0 1 0 0 0 1 0 1 0;
0 0 0 1 0 0 0 1 0 1;
0 0 0 0 1 1 0 0 1 0;
]

g = SimpleGraph(G)

gplot(g)

#+

k = nv(g)

k_colouring = Model(with_optimizer(GLPK.Optimizer))

@variable(k_colouring, z[1:k], Bin)
@variable(k_colouring, c[1:nv(g),1:k], Bin)
@constraint(k_colouring, [i = 1:nv(g)], sum(c[i,:]) == 1)
@constraint(k_colouring, [i = 1:nv(g), j = 1:nv(g), l = 1:k; G[i,j] == 1], c[i,l] + c[j,l] <= 1)
@constraint(k_colouring, [i = 1:nv(g), l = 1:k], c[i,l] <= z[l])

@objective(k_colouring, Min, sum(z))

optimize!(k_colouring)
@show value.(z);
@show value.(c);

#+

c = value.(c)
membership = zeros(nv(g))
for i in 1:nv(g)
    for j in 1:k
        if c[i,j] == 1
            membership[i] = j
            break
        end
    end
end
membership = convert(Array{Int},membership)

nodecolor = distinguishable_colors(nv(g), colorant"green")
nodefillc = nodecolor[membership]
gplot(g, nodefillc = nodefillc)

ext: jl


text:
using BenchmarkTools
using JLD

# Define a parent BenchmarkGroup to contain our suite
const suite = BenchmarkGroup()

# Add some child groups to our benchmark suite.
suite["utf8"] = BenchmarkGroup(["string", "unicode"])
suite["trigonometry"] = BenchmarkGroup(["math", "triangles"])

# This string will be the same every time because we're seeding the RNG
teststr = UTF8String(join(rand(MersenneTwister(1), 'a':'d', 10^4)))

# Add some benchmarks to the "utf8" group
suite["utf8"]["replace"] = @benchmarkable replace($teststr, "a", "b")
suite["utf8"]["join"] = @benchmarkable join($teststr, $teststr)

# Add some benchmarks to the "trigonometry" group
for f in (sin, cos, tan)
    for x in (0.0, pi)
        suite["trigonometry"][string(f), x] = @benchmarkable $(f)($x)
    end
end

# Load the suite's cached parameters as part of including the file. This is much
# faster and more reliable than re-tuning `suite` every time the file is included
paramspath = joinpath(Pkg.dir("BenchmarkTools"), "benchmark", "params.jld")
#tune!(suite); JLD.save(paramspath, "suite", params(suite));
loadparams!(suite, JLD.load(paramspath, "suite"), :evals, :samples);


text:
using OptiMode, BenchmarkTools
include("mpb_example.jl") # for now, just to load ε⁻¹

H,kz = solve_k(ω,ε⁻¹,Δx,Δy,Δz)
@benchmark solve_k($ω,$ε⁻¹,$Δx,$Δy,$Δz)
# BenchmarkTools.Trial:
#   memory estimate:  250.60 MiB
#   allocs estimate:  147758
#   --------------
#   minimum time:     1.178 s (0.97% GC)
#   median time:      1.385 s (0.80% GC)
#   mean time:        1.419 s (0.79% GC)
#   maximum time:     1.728 s (0.82% GC)
#   --------------
#   samples:          4
#   evals/sample:     1

ωs = collect(0.5:0.01:0.6)
Hs,ks = solve_k(ωs,ε⁻¹,Δx,Δy,Δz)
@benchmark solve_k($ωs,$ε⁻¹,$Δx,$Δy,$Δz)
# BenchmarkTools.Trial:
#   memory estimate:  1.23 GiB
#   allocs estimate:  1120050
#   --------------
#   minimum time:     5.371 s (1.19% GC)
#   median time:      5.371 s (1.19% GC)
#   mean time:        5.371 s (1.19% GC)
#   maximum time:     5.371 s (1.19% GC)
#   --------------
#   samples:          1
#   evals/sample:     1

H,ω = solve_ω(kz,ε⁻¹,Δx,Δy,Δz)
@benchmark solve_ω($kz,$ε⁻¹,$Δx,$Δy,$Δz)
# BenchmarkTools.Trial:
#   memory estimate:  114.77 MiB
#   allocs estimate:  58115
#   --------------
#   minimum time:     541.258 ms (1.01% GC)
#   median time:      748.920 ms (0.73% GC)
#   mean time:        787.952 ms (0.78% GC)
#   maximum time:     1.127 s (0.91% GC)
#   --------------
#   samples:          7
#   evals/sample:     1

ks = collect(1.5:0.01:1.6)
Hs,ωs = solve_ω(ks,ε⁻¹,Δx,Δy,Δz)
@benchmark solve_ω($ks,$ε⁻¹,$Δx,$Δy,$Δz)
# BenchmarkTools.Trial:
#   memory estimate:  632.98 MiB
#   allocs estimate:  309718
#   --------------
#   minimum time:     2.887 s (0.89% GC)
#   median time:      2.931 s (0.80% GC)
#   mean time:        2.931 s (0.80% GC)
#   maximum time:     2.976 s (0.71% GC)
#   --------------
#   samples:          2
#   evals/sample:     1


text:
module Brillouin

using Reexport
using Requires
using StaticArrays
# ---------------------------------------------------------------------------------------- #
const AVec = AbstractVector
const BasisLike{D} = AVec{<:SVector{D,<:Real}}
const SHOWDIGITS = 6
# ---------------------------------------------------------------------------------------- #
@enum BasisEnum CARTESIAN LATTICE

"""
    setting(x::Union{KPath, KPathInterpolant, Cell})

Return the basis setting of coordinates in `x`. The returned value is a member of the
[`BasisEnum`](@ref) enum with member values `LATTICE` (i.e. coordinates in the basis of the
lattice vectors) or `CARTESIAN` (i.e. coordinates in the Cartesian basis).
By default, methods in Brillouin will return coordinates in the `LATTICE` setting.
"""
setting(x) = x.setting[]

set_setting!(x, new_setting) = (x.setting[] = new_setting)

# ---------------------------------------------------------------------------------------- #
@doc """
    latticize!

In-place transform object with coordinates in a Cartesian basis to an object with
coordinates in a (explicitly or implicitly specified) lattice basis.
"""
function latticize! end
@doc """
    latticize

Transform an object with coordinates in a Cartesian basis to an object with coordinates in
a (explicitly or implicitly specified) lattice basis.
"""
function latticize end
latticize(v::AVec{<:Real}, basismatrix::AbstractMatrix{<:Real}) = basismatrix\v
latticize(v::AVec{<:Real}, basis::AVec{<:AVec{<:Real}}) = latticize(v, reduce(hcat, basis))
latticize(x) = (setting(x) === CARTESIAN ? latticize!(deepcopy(x)) : deepcopy(x))

# TODO: Move "root" definition of `cartesianize(!)` and `latticize(!)` to Bravais, so we
#       can share the definition with Crystalline
@doc """
    cartesianize!

In-place transform an object with coordinates in a (explicitly or implicitly specified)
lattice basis to an object with coordinates in a Cartesian basis.
"""
function cartesianize! end
@doc """
    cartesianize

Transform an object with coordinates in a (explicitly or implicitly specified) lattice basis
to an object with coordinates in a Cartesian basis.
"""
function cartesianize end
cartesianize(v::AVec{<:Real}, basis::AVec{<:AVec{<:Real}}) = v'basis
cartesianize(x) = (setting(x) === LATTICE ? cartesianize!(deepcopy(x)) : deepcopy(x))

"""
    basis(x::Union{KPath, KPathInterpolant, Cell})

Return the (reciprocal or direct) lattice basis associated with `x`, in Cartesian
coordinates.

Methods in Brillouin will by default return points in the lattice basis, i.e., points are
referred to `basis(x)`. This corresponds to the `setting(x) == LATTICE`.
Coordinates may instead be referred to a Cartesian basis, corresponding to
`setting(x) == CARTESIAN` by using [`cartesianize`](@ref). The result of
`basis(x)`, however, is invariant to this and always refers to the lattice basis in
Cartesian coordinates.
"""
basis(x) = x.basis

export latticize!, latticize,
    cartesianize!, cartesianize,
    basis, setting
# ---------------------------------------------------------------------------------------- #

# MAIN FUNCTIONALITY

include("KPaths.jl")
@reexport using .KPaths

include("WignerSeitz.jl")
@reexport using .WignerSeitz

# ---------------------------------------------------------------------------------------- #

function __init__()
    # plotting extensions on GLMakie load
    @require Makie="ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" begin
        include("requires/makie_wignerseitz.jl")
    end

    # plotting extensions on PlotlyJS load
    @require PlotlyJS="f0f68f2c-4968-5e81-91da-67840de0976a" begin
        include("requires/plotlyjs_wignerseitz.jl")
        include("requires/plotlyjs_kpaths.jl")
        include("requires/plotlyjs_dispersion.jl")
    end

    @require Spglib="f761d5c5-86db-4880-b97f-9680a7cccfb5" begin
        include("requires/spglib.jl")
    end
end

# ---------------------------------------------------------------------------------------- #

end


text:
"""
    IHeightable(nexid)

An interface for widgets whom height can modified.
"""
struct IHeightable <: AbstractINextion
    _nid::NexID

    function IHeightable(nid::NexID)
        new(nid)
    end
end

"""
    obj.height = new_val

Set Nextion object `obj` width to `new_val`.
"""
function setproperty!(obj::IHeightable, property::Symbol, new_val)
    if property == :height
        setnexproperty!(NexID(obj), :hig, new_val)
    else
        setfield!(obj, property, new_val)
    end
end


text:
import LinearAlgebra.cross
cross(a::Pt{2,T}, b::Pt{2,T}) where {T} = a[1]*b[2] - a[2]*b[1]

"""
Compute whether two flat patches of the same dimension overlap or not
"""
function overlap(p1::Simplex{U,D,C,N,T}, p2::Simplex{U,D,C,N,T}) where {U,D,C,N,T}

    # Are the patches in the same D-plane?
    throw(ErrorException("Not implemented yet!"))

end

###
# Version for points: overlap in this case is taken to mean that points approx conincide
###
function overlap(p::Simplex{U,0,C,1,T}, q::Simplex{U,0,C,1,T}) where {T,U,C}

    tol = sqrt(eps(T))
    return norm(p[1]-q[1]) < tol
end

"""
Compute whether two segments in 3D space overlap
"""
function overlap(p::Simplex{U,1,C,2,T}, q::Simplex{U,1,C,2,T}) where {T,U,C}

    tol = sqrt(eps(T))

    # Are these edges on the same line?
    o = q.vertices[2]
    u = p.vertices[1] - o
    t = q.tangents[1]
    t2 = dot(t,t)
    e = u - dot(u,t) / t2 * t

    # if vertex 1 of p is not in the line defined by q, they do not overlap
    norm(e) > tol && return false

    # if the tangents are not linearly dependent return false
    norm(cross(p.tangents[1],q.tangents[1])) > tol && return false

    # we determined p and q are on the same line. Do standard collision testing
    a1 = dot(p.vertices[1]-o, t)
    b1 = dot(p.vertices[2]-o, t)

    a1 < b1 ? (x1 = a1; y1 = b1) : (x1 = b1; y1 = a1)

    y1 <= zero(T) && return false # p to the left of q
    t2 <= x1 && return false      # q to the left of p

    return true
end

"""
Compute whether two triangles in 3D space overlap
"""
function overlap(p::Simplex{3,2,1,3,T}, q::Simplex{3,2,1,3,T}) where T

  tol = sqrt(eps(T))

  # Are the patches in the same plane?
  u1 = q.tangents[1]
  u2 = q.tangents[2]
  v = p.vertices[1] - q.vertices[2]

  # if vertex 1 of p is not in q, return false
  norm(dot(cross(u1,u2),v)) > tol && return false

  # if the two triangles are not coplanar return false
  norm(cross(p.normals[1], q.normals[1])) > tol && return false

  n = p.normals[1]
  for i in 1:3
    a = p.vertices[mod1(i+1,3)]
    b = p.vertices[mod1(i+2,3)]
    c = p.vertices[i]
    t = b - a
    m = cross(t,n)

    sp = zeros(T,3); sp[i] = dot(c-a, m)
    sq = T[dot(q.vertices[j]-a, m) for j in 1:3]

    minp, maxp = extrema(sp)
    minq, maxq = extrema(sq)

    maxq <= minp + tol && return false
    maxp <= minq + tol && return false
  end

  return true
end


text:
module Graph

using AutoHashEquals

abstract Pattern{C}

@auto_hash_equals type Choice{C} <: Pattern{C}
    patterns::Vector{Pattern{C}}
end

@auto_hash_equals type Sequence{C} <: Pattern{C}
    patterns::Vector{Pattern{C}}
end

@auto_hash_equals type Repeat{C} <: Pattern{C}
    pattern::Pattern{C}
    lo::Int
    hi::Int
end

@auto_hash_equals type Range{C} <: Pattern{C}
    lo::C
    hi::C
end

Literal{C}(c::C) = Range{C}(c, c)
Wild{C}(::Type{C}) = Range(typemin(C), typemax(C))

@auto_hash_equals type Group{C} <: Pattern{C}
    index::Int
    pattern::Pattern{C}
end

end


text:
using SpecialFunctions

all_funcs = (
    :~, :conj, :abs, :sin, :cos, :tan, :sinh, :cosh, :tanh, :asin, :acos, :atan,
    :asinh, :acosh, :atanh, :sec, :csc, :cot, :asec, :acsc, :acot, :sech, :csch,
    :coth, :asech, :acsch, :acoth, :sinc, :cosc, :cosd, :cotd, :cscd, :secd,
    :sind, :tand, :acosd, :acotd, :acscd, :asecd, :asind, :atand, :rad2deg,
    :deg2rad, :log, :log2, :log10, :log1p, :exponent, :exp, :exp2, :expm1,
    :cbrt, :sqrt, :ceil, :floor, :trunc, :round, :significand,
    :frexp, :ldexp, :modf, :real, :imag, :!, :identity,
    :zero, :one, :<<, :>>, :abs2, :sign, :sinpi, :cospi, :exp10,
    :iseven, :ispow2, :isfinite, :isinf, :isodd, :isinteger, :isreal,
    :isnan, :isempty, :iszero, :transpose, :copysign, :flipsign, :signbit,
    :+, :-, :*, :/, :\, :^, :(==), :(!=), :<, :(<=), :>, :(>=), :≈, :min, :max,
    :div, :fld, :rem, :mod, :mod1, :cmp, :&, :|, :xor,
    :clamp
)

special_funcs = (
    :airyai,
    :airyaiprime,
    :airybi,
    :airybiprime,
    :airyaix,
    :airyaiprimex,
    :airybix,
    :airybiprimex,
    :besselh,
    :besselhx,
    :besseli,
    :besselix,
    :besselj,
    :besselj0,
    :besselj1,
    :besseljx,
    :besselk,
    :besselkx,
    :bessely,
    :bessely0,
    :bessely1,
    :besselyx,
    :dawson,
    :erf,
    :erfc,
    :erfcinv,
    :erfcx,
    :erfi,
    :erfinv,
    :eta,
    :digamma,
    :invdigamma,
    :polygamma,
    :trigamma,
    :hankelh1,
    :hankelh1x,
    :hankelh2,
    :hankelh2x,
    :zeta
)

function to_abstract(x::Type{T}) where T <: Number
    isa(x, Union) && return to_abstract(reduce(Base.typejoin, Base.uniontypes(x)))
    AbstractNumber
end

function to_abstract(x::T) where T
    x == Any && return Number
    isa(x, TypeVar) && return to_abstract(x.ub)
    return x
end

function isa_number(T::DataType)
    return T <: Number || T == Any
end
function isa_number(tv::TypeVar)
    isa_number(tv.ub)
end

open(joinpath(@__DIR__, "overloads.jl"), "w") do io
    for (mod, funcs) in ((Base, all_funcs), (SpecialFunctions, special_funcs)), f in funcs
        func = getfield(mod, f)
        ms = methods(func)
        sigs = []
        for m in ms
            #TODO infer correct types
            sig = (Base.unwrap_unionall(m.sig).parameters...,)[2:end]
            if all(isa_number, sig)
                # _RT = Base.Core.Inference.return_type(func, Tuple{sig...})
                # isleaftype(_RT) && (RT = _RT)
                push!(sigs, m.nargs - 1)
            end
        end
        sigs = sort(unique(sigs))
        if f in (:+, :-, :*, :|, :&, :xor) && length(sigs) >= 2
            sigs = sigs[1:2]
        end
        for n in sigs
            (n == 3 && f == :cmp) && continue
            boolfunc = f in (
                :(==), :(!=), :<, :(<=), :>, :(>=), :≈
            ) || startswith(string(f), "is")

            argnames = ntuple(i-> Symbol("arg$i"), n)
            args = map(x-> :($x::AbstractNumber), argnames)
            converted = map(x-> :(number($x)), argnames)
            fname = "$mod.$(sprint(show, f))"
            ret = boolfunc ? "tmp" : "like(arg1, tmp)"
            println(io, """
            function $fname($(join(args, ", ")))
                tmp = $fname($(join(converted, ", ")))
                $ret
            end
            """)
            # I think it only makes sense to define all permutations only for binary ops.
            # The rest seems to be too bug ridden
            if n == 2
                ret = boolfunc ? "tmp" : "like(a, tmp)"
                println(io, """
                function $fname(a::AbstractNumber, b::Number)
                    tmp = $fname(number(a), b)
                    $ret
                end
                """)
                ret = boolfunc ? "tmp" : "like(b, tmp)"
                println(io, """
                function $fname(a::Number, b::AbstractNumber)
                    tmp = $fname(a, number(b))
                    $ret
                end
                """)
            end
        end
    end
end


text:
# in principle could be any EBayesShrinkageLocation, as long as it allows sample splitting.

# R: Regression method used
# R: EBayes method that will be applied within each fold

struct EBayesCrossFit{M <: MLJBase.Supervised, I <: Integer} <: AbstractEBayesMethod
    model::M
    num_folds::I
end

EBayesCrossFit(model; num_folds=5) = EBayesCrossFit(model, num_folds)

struct EBayesCrossFitResult{EBS <: NormalSamples,
                        XT,# dataframe
                        KF <: FoldsView, #split into folds
                        MF, #modelfits
                        EBF, #eb_fits
                        T,
                        EBCF}
    ss::EBS
    X::XT
    kfolds::KF
    reg_fits::MF
    eb_fits::EBF
    reg_preds::T
    eb_preds::T
    ebcf_method::EBCF
end

function fit(ebcf::EBayesCrossFit, Xs, ss::NormalSamples; verbosity=0)
    num_folds = ebcf.num_folds

    eb_preds = copy(response(ss))
    reg_preds = copy(response(ss))

    kf = kfolds(shuffleobs(1:length(ss)), num_folds)

    reg_fits = nothing # for now
    eb_fits = []

    # Initialize MLJ machine
    mlj_mach = MLJ.machine(ebcf.model, Xs, response(ss))

    for (train_idx, test_idx) in kf
        fit!(mlj_mach, rows=train_idx, verbosity=verbosity)
        if isa(ebcf.model, MLJBase.Deterministic)
            μs_pred = predict(mlj_mach, rows=test_idx)
        else
            μs_pred = predict_mean(mlj_mach, rows=test_idx)
        end
        eb_fit = fit(Normal(), SURE(), FixedLocation(μs_pred), ss[test_idx])
        push!(eb_fits, eb_fit)

        reg_preds[test_idx] = μs_pred
        eb_preds[test_idx] = predict(eb_fit)
    end

    EBayesCrossFitResult(ss, Xs, kf,
                reg_fits, eb_fits,
                reg_preds, eb_preds,
                ebcf)
end

function predict(ebcf_res::EBayesCrossFitResult)
    ebcf_res.eb_preds
end


text:
using QMTK
using Compat.Test

@testset "Constructors" begin

    space = SiteSpace(Bit, (2, 2); nflips=1)
    @test israndomized(space) == true
    typeof(space.data) <: Sites

    space = SiteSpace(Spin, 2, 2; nflips=1)
    @test israndomized(space) == true
    typeof(space.data) <: Sites

    space = SiteSpace(UnRandomized, space)
    @test israndomized(space) == false

end

@testset "methods" begin

    space = SiteSpace(Spin, 2, 2; nflips=1)
    space = reset!(space)
    @test israndomized(space) == false

    @test acquire(space) == space.data
    @test acquire(space) !== space.data

    @test_throws UnRandomizedError shake!(space)

    space = randomize!(space)
    sample = shake!(space)
    @test sample === space.data

    itr = traverse(space)
    for (i, each) in zip(0 : (1 << length(itr)) - 1, itr)
        @test Int(each) == i
    end

end


text:
"""
    coef(x::EquateMethod)

Return equating coefficient. The input object is equated result under the SG and NEAT design method.
"""
function coef(x::EquateMethod)
    if all(isnothing.(match.(r"estimates", string.(fieldnames(typeof(x))))))
        println("$(typeof(x)) has no field `estimates`")
    else
        x.estimates
    end
end


text:
# Autonomous van der Pol
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot

# parameters and initial values
dvec = [0.1, 5]    # damping parameters
y_0 = [0.1, 0.1]   # initial condition
ttrans = 100   # transient time
tattr = 16     # time on attractor
Δt = 0.01   # sampling time
tvec = 0:Δt:tattr

fig, axs = subplots(1,2; figsize=(0.6figx,0.8figy))

for id = 1:2
    d = dvec[id]
    y_init = y_0
    ds = Systems.vanderpol(; F = 0, μ = d)
    tr = trajectory(ds, tattr, y_0; Δt, Ttr = ttrans)
    x1vec, x2vec = columns(tr)
    axs[1].plot(tvec,x1vec; c = COLORS[id] )
    axs[2].plot(x1vec,x2vec; c = COLORS[id] )
end
axs[1].set_xlabel(L"$t$"; labelpad = -20)
axs[1].set_ylabel(L"$x$", rotation = 0)
axs[2].set_xlabel(L"$x$"; labelpad = -20)
axs[2].set_ylabel(L"$\dot x$", rotation = 0)
axs[1].set_xticks([0,15])
axs[2].set_xticks([-2,0,2])
axs[2].set_xticklabels(["-2", "", "2"])
add_identifiers!(fig)
fig.tight_layout(pad=0.33)

wsave(plotsdir("9", "vanderPol_autonomous"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using Plots

"""
`plot_phimatrix(outputfile::String, burnin::Int64, thin::Int64)`

Plots the mean phi values resulting from pmdi output
## Input
- `outputFile::String` a string referring to the location on disk of pmdi output
- `burnin::Int64` an integer of the number of initial iterations to discard as burn-in
- `thin::Int64` an integer for the rate at which to thin output, `thin = 2` discards
every second iteration

## Output
Outputs a heatmap of pairwise dataset Φ values
"""
function get_phi(outputFile::String, burnin::Int64 = 0, thin::Int64 = 0)
    outputNames = split(readline(outputFile), ',')
    phiColumns = map(outputNames) do str
        occursin(r"phi_", str)
    end
    output = readdlm(outputFile, ',', header = false, skipstart = burnin + 1)

    phiValues = Matrix{Float64}(output[1:thin:end, phiColumns])
    return phiValues
end

function plot_phi_matrix(outputFile::String, burnin::Int64 = 0, thin::Int64 = 1)
    phi_values = get_phi(outputFile, burnin, thin)
    K = Int64(0.5 + sqrt(8 * size(phi_values, 2) + 1) * 0.5)
    @assert K > 1 "Φ not inferred for no. of datasets = 1"
    phi_matrix = Matrix{Float64}(undef, K, K)
    phi_matrix .= NaN

    i = 1
    for k1 = 1:(K - 1)
        for k2 = (k1 + 1):K
            phi_matrix[k1, k2] = phi_matrix[k2, k1] = Statistics.mean(phi_values[:, i])
            i += 1
        end
    end

    Plots.heatmap(["phi (., $i)" for i = 1:K],
                  ["phi ($i, .)" for i = 1:K],
                  phi_matrix,
                  yflip = true,
                  widen = false,
                  left_margin = 20px,
                  c = :viridis,
                  clim = (0, maximum(phi_matrix)),
                  aspect_ratio = 1)
end

"""
`plot_phichain(outputfile::String, burnin::Int64, thin::Int64)`

Plots the phi values at each iteration resulting from pmdi output
## Input
- `outputFile::String` a string referring to the location on disk of pmdi output
- `burnin::Int64` an integer of the number of initial iterations to discard as burn-in
- `thin::Int64` an integer for the rate at which to thin output, `thin = 2` discards
every second iteration

## Output
Outputs a line plot of phi values resulting from pmdi output
"""
function plot_phi_chain(outputFile::String, burnin::Int64 = 0, thin::Int64 = 1)
    phi_values = get_phi(outputFile, burnin, thin)
    n_phis = size(phi_values, 2)
    K = Int64(0.5 + sqrt(8 * size(phi_values, 2) + 1) * 0.5)
    @assert K > 1 "Φ not inferred for no. of datasets = 1"
    phi_names = ["phi ($i, $j)" for i = 1:(K - 1) for j = (i + 1):K]
    Plots.plot(phi_values,
               layout = n_phis,
               legend = false,
               palette = :viridis,
               ylims = [0, maximum(phi_values)],
               title = [phi_names[i] for j = 1:1, i = 1:n_phis],
               xticks = false)
end


text:
module rocRAND

using ..HSA
using ..AMDGPU
using ..AMDGPU: librocrand, mark!, wait!

using CEnum

export rand_logn!, rand_poisson!, rand_logn, rand_poisson

include("librocrand_common.jl")
include("error.jl")
include("librocrand.jl")
function version()
    s = string(ROCRAND_VERSION)
    VersionNumber(join([string(s[1]), s[2:4], s[5:end]], '.'))
end

# stdlib Random integration
include("random.jl")

# TODO: implement thread local RNG like CUDA.jl
const ROCRAND_RNG = Ref{Union{Nothing,RNG}}(nothing)

function default_rng()
    if ROCRAND_RNG[] == nothing
        ROCRAND_RNG[] = RNG()
        Random.seed!(ROCRAND_RNG[])
    end
    return ROCRAND_RNG[]::RNG
end

end


text:
  import LinearAlgebra

  function events!(domain, evals, nframes, nsvt)
    for i in 1:nsvt
      @inbounds for dt in nframes-1:-1:0
        evals[dt + 1, i] = LinearAlgebra.dot(@view(domain[1:nframes-dt,i]),@view(domain[dt+1:nframes,i]))
      end
    end
  end


text:
#=
The following tests check the distribution of `rand(RotType)`.
To visualize the distribution, try the following script.

    using Rotations, Plots
    RotType = RotMatrix3{Float64}

    # Sampling
    n = 10000
    rs = rand(RotType,n)
    q = rand(RotType)
    angles = rotation_angle.([r/q for r in rs])/π

    # Cumulative distribution function
    cdf_sampled(t) = count(<(t), angles)/n
    cdf_true(t) = t-sinpi(t)/π

    # Probability density function and histogram of sampling distribution
    histogram(angles, normed=true, label=false)
    plot!(t->1-cospi(t), 0, 1, label=false)

    # Cumulative distribution function
    plot(cdf_sampled, 0, 1, label=false)
    plot!(cdf_true, 0, 1, label=false)
=#

@testset "Distribution" begin

    # TODO: consider to remove Euler rotations (#155)
    Type_SO3 =  (RotMatrix3, RotMatrix{3}, AngleAxis, RotationVec,
                QuatRotation, MRP, RodriguesParam,
                RotXYZ, RotYZX, RotZXY, RotXZY, RotYXZ, RotZYX,
                RotXYX, RotYZY, RotZXZ, RotXZX, RotYXY, RotZYZ)

    Type_SO2 = (RotMatrix2, RotMatrix{2}, Angle2d, RotX, RotY, RotZ)

    # Number of sampling
    N = 100000

    @testset "Distribution $(RotType)" for RotType in Type_SO3
        # Sampling
        rs = rand(RotType, N)
        q = rand(RotType)
        angles = rotation_angle.([r/q for r in rs])/π

        # Check sampled rotations are rotations
        @test all(isrotation.(rs))
        @test isrotation(q)

        # Cumulative distribution function
        cdf_sampled(t) = count(<(t), angles)/N
        cdf_true(t) = t-sinpi(t)/π

        # Check the CDFs are approximately equal
        ts = 0:0.01:1
        @test norm(cdf_true.(ts) - cdf_sampled.(ts), Inf) < 0.01
    end

    @testset "distribution $(RotType)" for RotType in Type_SO2
        # Sampling
        rs = rand(RotType, N)
        q = rand(RotType)
        angles = rotation_angle.([r/q for r in rs])/2π
        angles = mod.(angles, 1)

        # Check sampled rotations are rotations
        @test all(isrotation.(rs))
        @test isrotation(q)

        # Cumulative distribution function
        cdf_sampled(t) = count(<(t), angles)/N
        cdf_true(t) = t

        # Check the CDFs are approximately equal
        ts = 0:0.01:1
        @test norm(cdf_true.(ts) - cdf_sampled.(ts), Inf) < 0.01
    end

    @testset "distribution $(N)-dimensinal RotMatrix" for N in 2:5
        # Sampling
        RotType = RotMatrix{N}
        n = 1000000
        rs = rand(RotType, n)
        q = rand(RotType)
        norms1 = norm.([(r/q - one(SMatrix{N,N})) for r in rs])
        norms2 = norm.([(r - one(SMatrix{N,N})) for r in rs])

        # Check sampled rotations are rotations
        @test all(isrotation.(rs))
        @test isrotation(q)

        # Cumulative distribution function
        cdf_sampled1(t) = count(<(t), norms1)/n
        cdf_sampled2(t) = count(<(t), norms2)/n

        # Check the CDFs are approximately equal
        ts = 0:0.01:1
        @test norm(cdf_sampled1.(ts) - cdf_sampled2.(ts), Inf) < 0.01
    end
end


text:
module SpikeSorter
using Compat
using DSP
using MultivariateStats
using StatsBase
import StatsBase.transform, StatsBase.fit

include("types.jl")
include("utility.jl")
include("hmm.jl")
include("viterbi.jl")
include("readfiles.jl")
include("features.jl")
include("extraction.jl")
include("io.jl")

#import GUICheck
#if GUICheck.hasgui()
#	include("plot.jl")
#end

export SpikeWaveforms

end # module


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

mutable struct PacketCaptureQueryStatusResult <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    captureStartTime::Any # spec type: Union{ Nothing, DateTime } # spec name: captureStartTime
    packetCaptureStatus::Any # spec type: Union{ Nothing, String } # spec name: packetCaptureStatus
    stopReason::Any # spec type: Union{ Nothing, String } # spec name: stopReason
    packetCaptureError::Any # spec type: Union{ Nothing, Vector{String} } # spec name: packetCaptureError

    function PacketCaptureQueryStatusResult(;name=nothing, id=nothing, captureStartTime=nothing, packetCaptureStatus=nothing, stopReason=nothing, packetCaptureError=nothing)
        o = new()
        validate_property(PacketCaptureQueryStatusResult, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(PacketCaptureQueryStatusResult, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(PacketCaptureQueryStatusResult, Symbol("captureStartTime"), captureStartTime)
        setfield!(o, Symbol("captureStartTime"), captureStartTime)
        validate_property(PacketCaptureQueryStatusResult, Symbol("packetCaptureStatus"), packetCaptureStatus)
        setfield!(o, Symbol("packetCaptureStatus"), packetCaptureStatus)
        validate_property(PacketCaptureQueryStatusResult, Symbol("stopReason"), stopReason)
        setfield!(o, Symbol("stopReason"), stopReason)
        validate_property(PacketCaptureQueryStatusResult, Symbol("packetCaptureError"), packetCaptureError)
        setfield!(o, Symbol("packetCaptureError"), packetCaptureError)
        o
    end
end # type PacketCaptureQueryStatusResult

const _property_map_PacketCaptureQueryStatusResult = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("id")=>Symbol("id"), Symbol("captureStartTime")=>Symbol("captureStartTime"), Symbol("packetCaptureStatus")=>Symbol("packetCaptureStatus"), Symbol("stopReason")=>Symbol("stopReason"), Symbol("packetCaptureError")=>Symbol("packetCaptureError"))
const _property_types_PacketCaptureQueryStatusResult = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("id")=>"String", Symbol("captureStartTime")=>"DateTime", Symbol("packetCaptureStatus")=>"String", Symbol("stopReason")=>"String", Symbol("packetCaptureError")=>"Vector{String}")
Base.propertynames(::Type{ PacketCaptureQueryStatusResult }) = collect(keys(_property_map_PacketCaptureQueryStatusResult))
Swagger.property_type(::Type{ PacketCaptureQueryStatusResult }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_PacketCaptureQueryStatusResult[name]))}
Swagger.field_name(::Type{ PacketCaptureQueryStatusResult }, property_name::Symbol) =  _property_map_PacketCaptureQueryStatusResult[property_name]

const _allowed_PacketCaptureQueryStatusResult_packetCaptureStatus = ["NotStarted", "Running", "Stopped", "Error", "Unknown"]

const _allowed_PacketCaptureQueryStatusResult_packetCaptureError = ["InternalError", "AgentStopped", "CaptureFailed", "LocalFileFailed", "StorageFailed"]

function check_required(o::PacketCaptureQueryStatusResult)
    true
end

function validate_property(::Type{ PacketCaptureQueryStatusResult }, name::Symbol, val)
    if name === Symbol("packetCaptureStatus")
        Swagger.validate_param(name, "PacketCaptureQueryStatusResult", :enum, val, _allowed_PacketCaptureQueryStatusResult_packetCaptureStatus)
    end
    if name === Symbol("packetCaptureError")
        Swagger.validate_param(name, "PacketCaptureQueryStatusResult", :enum, val, _allowed_PacketCaptureQueryStatusResult_packetCaptureError)
    end
end


text:
using GridPotentialFlow
Δx = 0.01
Lx = 2.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx);

v = Vortex(0.0,0.0,1.0);

model = VortexModel(g,vortices=[v]);

s = streamfunction(model);
using Plots
plot(s,g,xlabel="x",ylabel="y")

q = curl(s);
plot(q,g,xlabel="x",ylabel="y")

function ψ_vortex!(ψ::Nodes{Dual,nx,ny},vortex::Vortex,g::PhysicalGrid) where {nx,ny}
    x,y = coordinates(ψ,g)
    for i in 2:nx-1, j in 2:ny-1
        r = sqrt((x[i]-vortex.x)^2+(y[j]-vortex.y)^2)
        ψ[i,j] = ψ_vortex(r,vortex.Γ)
    end
end

function ψ_vortex(r::Real,Γ::Real)
    return -Γ/(2π)*log(r)
end

nv = 4;
vl = Vortex.(-Lx/4 .+ 0.4*Lx*(rand(nv).-0.5),-Lx/4 .+ 0.4*Lx*(rand(nv).-0.5),0.5*rand(nv).+0.5);

grids = [PhysicalGrid(xlim,ylim,Lx/(nx-2)) for nx in [2^p for p in 5:9]];

using LinearAlgebra: norm

errors = []
gridspacings = []
for g in grids
    # LGF
    local model = VortexModel(g,vortices=vl)
    s_lgf = streamfunction(model);

    # Exact solution
    s_exact = Nodes(Dual,size(g))
    s_temp = Nodes(Dual,size(g))
    for v in vl
        ψ_vortex!(s_temp,v,g)
        s_exact += s_temp
    end

    # Bring to same reference level. A constant value can be added/subtracted to any potential flow solution
    s_lgf .-= s_lgf[g.I0[1],2];
    s_exact .-= s_exact[g.I0[1],2];

    error = s_lgf-s_exact
    idx = g.I0[1] + 1 : g.N[1] - 1 # Only look at top right corner
    push!(errors,norm(error[idx,idx])/norm(s_lgf[idx,idx]))
    push!(gridspacings,g.Δx)
end

firstorder=(errors[end]/gridspacings[end].*gridspacings).^1
secondorder=(sqrt(errors[end])/gridspacings[end].*gridspacings).^2
p=plot(gridspacings,errors,xaxis=:log,yaxis=:log,marker=:circle,lab="error",xlabel="dx",legend=:bottomright,title="Error")
plot!(gridspacings,firstorder,lab="1st order",linestyle=:dot,linecolor=:black)
plot!(gridspacings,secondorder,lab="2nd order",linestyle=:dash,linecolor=:black)

d = Lx/2
Γ = 1
v1 = Vortex(d/2,0.0,Γ);
v2 = Vortex(-d/2,0.0,Γ);

Vθ = Γ/(2*π*d); # Analytical tangential velocity
T = π*d/Vθ # Analytical period

model = VortexModel(g,vortices=[v1,v2]);

tspan = (0.0,T);

import OrdinaryDiffEq

#= To construct the `ODEProblem` from this package, we will have to provide a right-hand side function that computes the velocities $Ẋ$ of the vortices, which is the local flow velocity at their positions. These velocities are obtained using `vortexvelocities`, which regularizes the vorticity to the grid, solves the potential flow system, and interpolates the velocities from the grid to the vortex locations as

$\left(U_{q}, V_{q}\right)=\sum_{i, j} \mathsf{v}_{i j} d\left(\frac{\mathsf{x}_{i}-X_{q}}{\Delta x}\right)\left(\frac{\mathsf{y}_{j}-Y_{q}}{\Delta x}\right),$

where $v$ is the velocity field on the nodes.
=#

function rhs(X,model,t)
    setvortexpositions!(model,X)
    Ẋ = vortexvelocities!(model)
    return Ẋ
end

X = getvortexpositions(model)
prob = OrdinaryDiffEq.ODEProblem(rhs,X,tspan,model);
sol = OrdinaryDiffEq.solve(prob,dt=0.1,OrdinaryDiffEq.RK4(),dense=false,adaptive=false);

plot(sol.t,map(s->s.u[1],sol.u),xlabel='t',label="x₁")
plot!(sol.t,map(s->s.v[1],sol.u),label="y₁")

@testset "Corotating point vortices" begin
    @test isapprox(sol.u[end].u[1], 0.5; atol = 1e-2)
    @test isapprox(sol.u[end].u[2], -0.5; atol = 1e-2)
    @test isapprox(sol.u[end].v[1], 0.0; atol = 1e-2)
    @test isapprox(sol.u[end].v[2], 0.0; atol = 1e-2)
end

xlim = (-2,2);
ylim = (-2,2);
Δx = 0.05;
g = PhysicalGrid(xlim,ylim,Δx);

function vortexpatch!(vort,xc,yc,Γ,radius,nring)
    Δr = radius/(nring-1/2)
    dΓ = Γ/(1+8*nring*(nring-1)/2)
    push!(vort,Vortex(xc,yc,dΓ))
    for ir in 1:nring-1
        nθ = 8*ir
        for j = 0:nθ-1
            push!(vort,Vortex(xc + ir*Δr*cos(2π*j/nθ),yc + ir*Δr*sin(2π*j/nθ),dΓ))
        end
    end
    return vort
end

vortexpatch(xc,yc,Γ,radius,nring) = vortexpatch!(Vortex[],xc,yc,Γ,radius,nring)

d0 = d
anims = []
for r0 in [0.2*d0, 0.4*d0]
    vortices = vcat(vortexpatch(0.0,0.0+d0/2,Γ,r0,10),vortexpatch(0.0,0.0-d0/2,Γ,r0,10));
    vortexcolors = vcat(fill(:blue,length(vortices)÷2),fill(:red,length(vortices)÷2));
    local model = VortexModel(g,vortices=vortices);
    local X = getvortexpositions(model)
    local prob = OrdinaryDiffEq.ODEProblem(rhs,X,tspan,model);
    local sol = OrdinaryDiffEq.solve(prob,dt=0.1,OrdinaryDiffEq.RK4(),dense=false,adaptive=false);
    anim = @animate for i=1:length(sol.t)-1
        plot(xlims=xlim,ylims=ylim,ratio=:equal,legend=:none,title="r0/d0 = $(r0/d0)")
        scatter!(sol.u[i].u,sol.u[i].v,markerstrokewidth=0,markersize=3,color=vortexcolors)
    end
    push!(anims,anim)
end

gif(anims[1])

gif(anims[2])

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using BenchmarkTools
globalsuite = BenchmarkGroup()

include("profile_ylm.jl")

# judge(minimum(results), minimum(results1)) |> display


text:
using ONNX, Flux, ProtoBuf
include("ops_tests.jl")

#test add:
main_test("$ONNX_TEST_PATH/test_add",
    read_output("$ONNX_TEST_PATH/test_add"),
        read_input("$ONNX_TEST_PATH/test_add")[1],
            read_input("$ONNX_TEST_PATH/test_add")[2])

#test add bcast
main_test("$ONNX_TEST_PATH/test_add_bcast",
          read_output("$ONNX_TEST_PATH/test_add_bcast"),
          read_input("$ONNX_TEST_PATH/test_add_bcast")[1],
          read_input("$ONNX_TEST_PATH/test_add_bcast")[2])

#test mul
main_test("$ONNX_TEST_PATH/test_mul",
    read_output("$ONNX_TEST_PATH/test_mul"),
        read_input("$ONNX_TEST_PATH/test_mul")[1],
            read_input("$ONNX_TEST_PATH/test_mul")[2])

#test mul bcast
main_test("$ONNX_TEST_PATH/test_mul_bcast",
    read_output("$ONNX_TEST_PATH/test_mul_bcast"),
        read_input("$ONNX_TEST_PATH/test_mul_bcast")[1],
            read_input("$ONNX_TEST_PATH/test_mul_bcast")[2])

#test sub
main_test("$ONNX_TEST_PATH/test_sub",
    read_output("$ONNX_TEST_PATH/test_sub"),
        read_input("$ONNX_TEST_PATH/test_sub")[1],
            read_input("$ONNX_TEST_PATH/test_sub")[2])

#test sub bcast
main_test("$ONNX_TEST_PATH/test_sub_bcast",
    read_output("$ONNX_TEST_PATH/test_sub_bcast"),
        read_input("$ONNX_TEST_PATH/test_sub_bcast")[1],
            read_input("$ONNX_TEST_PATH/test_sub_bcast")[2])

#test div
main_test("$ONNX_TEST_PATH/test_div",
    read_output("$ONNX_TEST_PATH/test_div"),
        read_input("$ONNX_TEST_PATH/test_div")[1],
            read_input("$ONNX_TEST_PATH/test_div")[2])

#test div bcast
main_test("$ONNX_TEST_PATH/test_div_bcast",
    read_output("$ONNX_TEST_PATH/test_div_bcast"),
        read_input("$ONNX_TEST_PATH/test_div_bcast")[1],
            read_input("$ONNX_TEST_PATH/test_div_bcast")[2])

#test matmul 2d
main_test("$ONNX_TEST_PATH/test_matmul_2d",
    read_output("$ONNX_TEST_PATH/test_matmul_2d"),
        read_input("$ONNX_TEST_PATH/test_matmul_2d")[1],
            read_input("$ONNX_TEST_PATH/test_matmul_2d")[2])
#test exp
main_test("$ONNX_TEST_PATH/test_exp",
    read_output("$ONNX_TEST_PATH/test_exp"),
        read_input("$ONNX_TEST_PATH/test_exp")[1])

#test reciprocal
main_test("$ONNX_TEST_PATH/test_reciprocal",
    read_output("$ONNX_TEST_PATH/test_reciprocal"),
        read_input("$ONNX_TEST_PATH/test_reciprocal")[1])

#test reciprocal example
main_test("$ONNX_TEST_PATH/test_reciprocal_example",
    read_output("$ONNX_TEST_PATH/test_reciprocal_example"),
        read_input("$ONNX_TEST_PATH/test_reciprocal_example")[1])

#test floor
main_test("$ONNX_TEST_PATH/test_floor",
    read_output("$ONNX_TEST_PATH/test_floor"),
        read_input("$ONNX_TEST_PATH/test_floor")[1])

#test ceil
main_test("$ONNX_TEST_PATH/test_ceil",
    read_output("$ONNX_TEST_PATH/test_ceil"),
        read_input("$ONNX_TEST_PATH/test_ceil")[1])

#test log
main_test("$ONNX_TEST_PATH/test_log",
    read_output("$ONNX_TEST_PATH/test_log"),
        read_input("$ONNX_TEST_PATH/test_log")[1])

#test pow
main_test("$ONNX_TEST_PATH/test_pow",
    read_output("$ONNX_TEST_PATH/test_pow"),
        read_input("$ONNX_TEST_PATH/test_pow")[1],
            read_input("$ONNX_TEST_PATH/test_pow")[2])

#test pow bcast
main_test("$ONNX_TEST_PATH/test_pow_bcast_array",
    read_output("$ONNX_TEST_PATH/test_pow_bcast_array"),
        read_input("$ONNX_TEST_PATH/test_pow_bcast_array")[1],
            read_input("$ONNX_TEST_PATH/test_pow_bcast_array")[2])

#test pow bcast scalar
main_test("$ONNX_TEST_PATH/test_pow_bcast_scalar",
    read_output("$ONNX_TEST_PATH/test_pow_bcast_scalar"),
        read_input("$ONNX_TEST_PATH/test_pow_bcast_scalar")[1],
            read_input("$ONNX_TEST_PATH/test_pow_bcast_scalar")[2])

#test pow example
main_test("$ONNX_TEST_PATH/test_pow_example",
    read_output("$ONNX_TEST_PATH/test_pow_example"),
        read_input("$ONNX_TEST_PATH/test_pow_example")[1],
            read_input("$ONNX_TEST_PATH/test_pow_example")[2])

#test relu
main_test("$ONNX_TEST_PATH/test_relu",
    read_output("$ONNX_TEST_PATH/test_relu"),
        read_input("$ONNX_TEST_PATH/test_relu")[1])

#Test sum one input
main_test("$ONNX_TEST_PATH/test_sum_one_input",
          read_output("$ONNX_TEST_PATH/test_sum_one_input"),
          read_input( "$ONNX_TEST_PATH/test_sum_one_input")[1])

#Test sum two inputs
main_test("$ONNX_TEST_PATH/test_sum_two_inputs",
          read_output("$ONNX_TEST_PATH/test_sum_two_inputs"),
          read_input("$ONNX_TEST_PATH/test_sum_two_inputs")[1],
          read_input("$ONNX_TEST_PATH/test_sum_two_inputs")[2])

# Test Prelu example
main_test("$ONNX_TEST_PATH/test_prelu_example", read_output("$ONNX_TEST_PATH/test_prelu_example"),
    read_input("$ONNX_TEST_PATH/test_prelu_example")[1],
        read_input("$ONNX_TEST_PATH/test_prelu_example")[2])

# Test PRelu broadcast
main_test("$ONNX_TEST_PATH/test_prelu_broadcast", read_output("$ONNX_TEST_PATH/test_prelu_broadcast"),
    read_input("$ONNX_TEST_PATH/test_prelu_broadcast")[1],
        read_input("$ONNX_TEST_PATH/test_prelu_broadcast")[2])

## Trigonometric ops

#Test sin
main_test("$ONNX_TEST_PATH/test_sin",
          read_output("$ONNX_TEST_PATH/test_sin"),
          read_input("$ONNX_TEST_PATH/test_sin")[1])

main_test("$ONNX_TEST_PATH/test_sin_example",
          read_output("$ONNX_TEST_PATH/test_sin_example"),
          read_input("$ONNX_TEST_PATH/test_sin_example")[1])
#Test cos
main_test("$ONNX_TEST_PATH/test_cos",
          read_output("$ONNX_TEST_PATH/test_cos"),
          read_input("$ONNX_TEST_PATH/test_cos")[1])

main_test("$ONNX_TEST_PATH/test_cos_example",
          read_output("$ONNX_TEST_PATH/test_cos_example"),
          read_input("$ONNX_TEST_PATH/test_cos_example")[1])

#Test tan
main_test("$ONNX_TEST_PATH/test_tan",
          read_output("$ONNX_TEST_PATH/test_tan"),
          read_input("$ONNX_TEST_PATH/test_tan")[1])

main_test("$ONNX_TEST_PATH/test_tan_example",
          read_output("$ONNX_TEST_PATH/test_tan_example"),
          read_input("$ONNX_TEST_PATH/test_tan_example")[1])

#test asin
main_test("$ONNX_TEST_PATH/test_asin",
          read_output("$ONNX_TEST_PATH/test_asin"),
          read_input("$ONNX_TEST_PATH/test_asin")[1])

main_test("$ONNX_TEST_PATH/test_asin_example",
          read_output("$ONNX_TEST_PATH/test_asin_example"),
          read_input("$ONNX_TEST_PATH/test_asin_example")[1])

#test acos
main_test("$ONNX_TEST_PATH/test_acos",
          read_output("$ONNX_TEST_PATH/test_acos"),
          read_input("$ONNX_TEST_PATH/test_acos")[1])

main_test("$ONNX_TEST_PATH/test_acos_example",
          read_output("$ONNX_TEST_PATH/test_acos_example"),
          read_input("$ONNX_TEST_PATH/test_acos_example")[1])

#test atan
main_test("$ONNX_TEST_PATH/test_atan",
          read_output("$ONNX_TEST_PATH/test_atan"),
          read_input("$ONNX_TEST_PATH/test_atan")[1])

main_test("$ONNX_TEST_PATH/test_atan_example",
          read_output("$ONNX_TEST_PATH/test_atan_example"),
          read_input("$ONNX_TEST_PATH/test_atan_example")[1])

# Flatten axis 0
main_test("$ONNX_TEST_PATH/test_flatten_axis0",
    read_output("$ONNX_TEST_PATH/test_flatten_axis0"),
        read_input("$ONNX_TEST_PATH/test_flatten_axis0")[1])

# Flatten axis 1
main_test("$ONNX_TEST_PATH/test_flatten_axis1",
    read_output("$ONNX_TEST_PATH/test_flatten_axis1"),
        read_input("$ONNX_TEST_PATH/test_flatten_axis1")[1])

# Flatten axis 2
main_test("$ONNX_TEST_PATH/test_flatten_axis2",
    read_output("$ONNX_TEST_PATH/test_flatten_axis2"),
        read_input("$ONNX_TEST_PATH/test_flatten_axis2")[1])

# Flatten axis 3
main_test("$ONNX_TEST_PATH/test_flatten_axis3",
    read_output("$ONNX_TEST_PATH/test_flatten_axis3"),
        read_input("$ONNX_TEST_PATH/test_flatten_axis3")[1])

# Flatten default axis
main_test("$ONNX_TEST_PATH/test_flatten_default_axis",
    read_output("$ONNX_TEST_PATH/test_flatten_default_axis"),
        read_input("$ONNX_TEST_PATH/test_flatten_default_axis")[1])

# test gemm broadcast
main_test("$ONNX_TEST_PATH/test_gemm_broadcast", read_output("$ONNX_TEST_PATH/test_gemm_broadcast"),
     read_input("$ONNX_TEST_PATH/test_gemm_broadcast")[1],
     read_input("$ONNX_TEST_PATH/test_gemm_broadcast")[2],
     read_input("$ONNX_TEST_PATH/test_gemm_broadcast")[3])

# test gemm nobroadcast
main_test("$ONNX_TEST_PATH/test_gemm_nobroadcast", read_output("$ONNX_TEST_PATH/test_gemm_nobroadcast"),
     read_input("$ONNX_TEST_PATH/test_gemm_nobroadcast")[1],
     read_input("$ONNX_TEST_PATH/test_gemm_nobroadcast")[2],
     read_input("$ONNX_TEST_PATH/test_gemm_nobroadcast")[3])

# test unsqueeze
main_test("$ONNX_TEST_PATH/test_unsqueeze", read_output("$ONNX_TEST_PATH/test_unsqueeze"),
     read_input("$ONNX_TEST_PATH/test_unsqueeze")[1])

# test abs
main_test("$ONNX_TEST_PATH/test_abs", read_output("$ONNX_TEST_PATH/test_abs"),
        read_input("$ONNX_TEST_PATH/test_abs")[1])

# test clip
main_test("$ONNX_TEST_PATH/test_clip", read_output("$ONNX_TEST_PATH/test_clip"),
    read_input("$ONNX_TEST_PATH/test_clip")[1])

# test clip default inbounds
main_test("$ONNX_TEST_PATH/test_clip_default_inbounds",
    read_output("$ONNX_TEST_PATH/test_clip_default_inbounds"),
        read_input("$ONNX_TEST_PATH/test_clip_default_inbounds")[1])

# test clip default max
main_test("$ONNX_TEST_PATH/test_clip_default_max",
    read_output("$ONNX_TEST_PATH/test_clip_default_max"),
        read_input("$ONNX_TEST_PATH/test_clip_default_max")[1])

# test clip default min
main_test("$ONNX_TEST_PATH/test_clip_default_min",
    read_output("$ONNX_TEST_PATH/test_clip_default_min"),
        read_input("$ONNX_TEST_PATH/test_clip_default_min")[1])

# test clip example
main_test("$ONNX_TEST_PATH/test_clip_example",
    read_output("$ONNX_TEST_PATH/test_clip_example"),
        read_input("$ONNX_TEST_PATH/test_clip_example")[1])

# test clip inbounds
main_test("$ONNX_TEST_PATH/test_clip_inbounds",
    read_output("$ONNX_TEST_PATH/test_clip_inbounds"),
        read_input("$ONNX_TEST_PATH/test_clip_inbounds")[1])

# test clip outbounds
main_test("$ONNX_TEST_PATH/test_clip_outbounds",
    read_output("$ONNX_TEST_PATH/test_clip_outbounds"),
        read_input("$ONNX_TEST_PATH/test_clip_outbounds")[1])

# test clip splitbounds
main_test("$ONNX_TEST_PATH/test_clip_splitbounds",
    read_output("$ONNX_TEST_PATH/test_clip_splitbounds"),
        read_input("$ONNX_TEST_PATH/test_clip_splitbounds")[1])


text:
export saveModel, loadModel
using JSON

"""
Function to save a model instance to a jsonFile

    saveModel(ModelName, \"jsonFilenameAsString.json\")
"""
function saveModel(M::Model, jsonFile::AbstractString)
    # Save a Dict as json file
    f = open(jsonFile, "w")
    write(f, JSON.json(M))
    close(f)
end

"""
Function to load a model instance from a json file

    M = loadModel(\"jsonFilenameAsString.json\")
"""
function loadModel(jsonFile::AbstractString)
    # Save a Dict as json file
    f = open(jsonFile, "r")
    D = JSON.parse(f)
    Sims = Dict([k => SimulationData(k, D["sims"][k]["I"],D["sims"][k]["P"], D["sims"][k]["D"]) for k in keys(D["sims"])])
    close(f)
    M = Model(D["odes"], D["init"], D["pars"]; name =  D["name"], aux = D["aux"], spec = D["spec"], vars = D["vars"], sims = Sims)
    return(M)
end

max_stars_repo_name: novakgroupoxford/XPPjl

max_issues_repo_name: novakgroupoxford/XPPjl

max_forks_repo_name: novakgroupoxford/XPPjl


text:
for N in 2:11
    facts("Hypercube in $N dimensions") do
        m = Model()
        @variable(m, 0 ≤ x[1:N] ≤ 1)

        v, r = get_extrema(m)

        @fact length(v) --> 2^N
        @fact length(r) --> 0

        for k in 0:N
            for p in combinations(1:N, k)
                ex = zeros(N)
                ex[p] = ones(length(p))
                @fact ex in v --> true
            end
        end
    end
end


text:
# insert at the given index, "123" inserted '_' at 2 would end up "1_23"
function insert_char(og_string, index, character)
    return string(og_string[1:index - 1], character, og_string[index:length(og_string)])
end

# if given an input sequence and a corresponding chromosome representation
# rebuild the string by inserting gaps into the gap positions labeled in the chrom_rep
function build_strings(input_sequence, chrom_rep)
    full_str = []
    for i = 1:length(input_sequence)
        my_string = input_sequence[i]

        for j = 1:length(chrom_rep[i])
           my_string = insert_char(my_string, chrom_rep[i][j], '_')
        end
        push!(full_str, my_string)
    end

    return full_str
end

# take the input sequence and growth fraction and calculate the extension to the longest string in the sequence
# then create a list of the number of gaps needed to reach that length for each string in the sequence
# return gap list and length
function get_gap_count(input_sequence, growth)

    # set our max to negative, will be overwritten by a length, which is always positive
    mx = -1

    # store the lengths of each sequence to determine number of gaps later
    gap_counts = []

    # go through each input sequence and overwrite mx if the length is larger than what we have seen already
    for i = 1:length(input_sequence)
        # put the length of each sequence into gap counts to be updated once we know the final length
        push!(gap_counts, length(input_sequence[i]))
        mx = max(mx, length(input_sequence[i]))
    end

    # this length will be what we put the actual size of all of our strings to, including the gaps
    len = ceil(growth * mx)

    # change the lengths to be the number of gaps needed to reach the calculated length
    for i = 1:length(gap_counts)
        gap_counts[i] = Int(len - gap_counts[i])
    end

    # return gap count for each string needed to extend it to calculated length
    return gap_counts, Int(len)
end

function msa_from_file(file)
    score_list = []
    input_sequence = []
    open(file, "r") do io
        lines = readlines(io)
        push!(score_list, parse(Int64, lines[1]))
        push!(score_list, parse(Int64, lines[2]))
        push!(score_list, parse(Int64, lines[3]))

        range = 4 + parse(Int64, lines[4]) * 2
        for i = 5:2:range
            push!(input_sequence, [lines[i], lines[i + 1]])
        end
    end;

    println(score_list)
    println(input_sequence)

    return score_list, input_sequence
end

function pairwise_score(reference_sequence, input_sequence)
    sequence_length = length(reference_sequence[1])
    score = 0

    match = 0
    mismatch = 0
    total = 0

    ref_list = []
    index_list = fill(1, length(reference_sequence))

    # for every sequence but the last
    for i = 1:length(reference_sequence) - 1
        # go through all characters (columns)
        for pos = 1:sequence_length
            # only check the characters
            if reference_sequence[i][pos] == '_'
                continue
            end

            for j = i + 1:length(reference_sequence)
                if reference_sequence[j][pos] == '_'
                    continue
                end

                count = 0
                for k = 1:pos
                    if reference_sequence[j][k] == '_'
                        continue
                    end

                    count += 1
                end

                push!(ref_list, [ [i, index_list[i]], [j, count] ])
                total += 1
            end

            index_list[i] += 1
        end
    end

    ref_list2 = []
    index_list2 = fill(1, length(input_sequence))

    # for every sequence but the last
    for i = 1:length(input_sequence) - 1
        # go through all characters (columns)
        for pos = 1:sequence_length
            # only check the characters
            if input_sequence[i][pos] == '_'
                continue
            end

            for j = i + 1:length(input_sequence)
                if input_sequence[j][pos] == '_'
                    continue
                end

                count = 0
                for k = 1:pos
                    if input_sequence[j][k] == '_'
                        continue
                    end

                    count += 1
                end

                push!(ref_list2, [ [i, index_list2[i]], [j, count] ])
                total += 1
            end

            index_list2[i] += 1
        end
    end

    for i = 1:length(ref_list)
        if ref_list[i] in ref_list2
            println(ref_list[i])
            match += 1
        end
    end

    score = (match) / total

    #println(ref_list)
    println(score)

    return score
end

export insert_char
export build_strings
export get_gap_count
export msa_from_file
export pairwise_score


text:
"""
This enum is used to select the distance calculation type.
We converted the name into a string, and the magic happens! :D

# DistanceType Possibilities:

    Euclidean, Cityblock, TotalVariation, Chebyshev,
    Jaccard, BrayCurtis, SpanNormDist

# Examples:

    Euclidean::DistanceType
    Cityblock::DistanceType
"""
@enum DistanceType begin

    Euclidean = 1
    Cityblock = 2
    TotalVariation = 3
    Chebyshev = 4
    Jaccard = 5
    BrayCurtis = 6
    SpanNormDist = 7

end

#=

    This dictionary is used with the enum DistanceType to
    select the distance calculation.
    After convert the enum in a string we applied in this dictionary. :)

=#

#add
# Creating the dictionary

distanceFormulas = Dict()

# Setting Euclidean
temp_euc(x, y) = sqrt(sum((x - y) .^ 2))
distanceFormulas["Euclidean"] = temp_euc

# Setting Cityblock
temp_cb(x, y) = sum(abs.(x - y))
distanceFormulas["Cityblock"] = temp_cb

# Setting TotalVariation
temp_tv(x, y) = sum(abs.(x - y)) / 2
distanceFormulas["TotalVariation"] = temp_tv

# Setting Chebyshev
temp_chby(x, y) = maximum(abs.(x - y))
distanceFormulas["Chebyshev"] = temp_chby

# Setting Jaccard
temp_jcc(x, y) = 1 - sum(min(x, y)) / sum(max(x, y))
distanceFormulas["Jaccard"] = temp_jcc

# Setting BrayCurtis
temp_bc(x, y) = sum(abs.(x - y)) / sum(abs.(x + y))
distanceFormulas["BrayCurtis"] = temp_bc

# Setting SpanNormDist
temp_snd(x, y) = maximum(x - y) - minimum(x - y)
distanceFormulas["SpanNormDist"] = temp_snd

"""
    GetDistance(distanceType::DistanceType,
    v1::AbstractArray{<: Number}, v2::AbstractArray{<: Number})

This function compute distance between two numeric vectors. You can use
diferents distance calculations selecting with "DistanceType".

# Examples:

    GetDistance(Euclidean::DistanceType, [2,2], [2,4])
    GetDistance(Cityblock::DistanceType, [2,2], [2,4])
    GetDistance(TotalVariation::DistanceType, [2,2], [2,4])
"""
function GetDistance(distanceType::DistanceType, v1::AbstractArray{<: Number}, v2::AbstractArray{<: Number})
    if length(v1) != length(v2)
        throw(DimensionMismatch("The dimensions of the vectors are different. Therefore, impossible to calculate!\n
        Vector1: $(length(v1)) // Vector2: $(length(v2))."))
    end

    return distanceFormulas[string(distanceType)](v1, v2)
end

# Formulas from: https://github.com/JuliaStats/Distances.jl


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Hanabi_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Hanabi")
JLLWrappers.@generate_main_file("Hanabi", UUID("37516f07-4740-538e-b7cb-eec320f4a184"))
end  # module Hanabi_jll


text:
# domain[state => f] shorthand for evaluating `f` within a domain and state
Base.getindex(domain::Domain, statevar::Pair{<:State, <:Term}) =
    evaluate(domain, first(statevar), last(statevar))
Base.getindex(domain::Domain, statevar::Pair{<:State, Symbol}) =
    evaluate(domain, first(statevar), Const(last(statevar)))
Base.getindex(domain::Domain, statevar::Pair{<:State, String}) =
    evaluate(domain, first(statevar), Parser.parse_formula(last(statevar)))

get_objects(domain::Domain, state::State) = get_objects(state)

function get_objects(domain::Domain, state::State, type::Symbol)
    return Const[o for (o, ty) in get_objtypes(state)
                 if ty == type || ty in get_all_subtypes(domain, type)]
end

"Returns all (recursive) subtypes of `type` in a domain."
function get_all_subtypes(domain::Domain, type::Symbol)
    subtypes = get_subtypes(domain, type)
    if isempty(subtypes) return subtypes end
    return reduce(vcat, [get_all_subtypes(domain, ty) for ty in subtypes],
                  init=subtypes)
end