text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

struct NullException <: Exception
end

"""
    Nullable(x, hasvalue::Bool=true)

Wrap value `x` in an object of type `Nullable`, which indicates whether a value is present.
`Nullable(x)` yields a non-empty wrapper and `Nullable{T}()` yields an empty instance of a
wrapper that might contain a value of type `T`.

`Nullable(x, false)` yields `Nullable{typeof(x)}()` with `x` stored in the result's `value`
field.

# Examples

```jldoctest
julia> Nullable(1)
Nullable{Int64}(1)

julia> Nullable{Int64}()
Nullable{Int64}()

julia> Nullable(1, false)
Nullable{Int64}()

julia> dump(Nullable(1, false))
Nullable{Int64}
  hasvalue: Bool false
  value: Int64 1
```
"""
Nullable(value::T, hasvalue::Bool=true) where {T} = Nullable{T}(value, hasvalue)
Nullable() = Nullable{Union{}}()

eltype(::Type{Nullable{T}}) where {T} = T

convert(::Type{Nullable{T}}, x::Nullable{T}) where {T} = x
convert(::Type{Nullable   }, x::Nullable   ) = x

convert(t::Type{Nullable{T}}, x::Any) where {T} = convert(t, convert(T, x))

function convert(::Type{Nullable{T}}, x::Nullable) where T
    return isnull(x) ? Nullable{T}() : Nullable{T}(convert(T, get(x)))
end

convert(::Type{Nullable{T}}, x::T) where {T<:Nullable} = Nullable{T}(x)
convert(::Type{Nullable{T}}, x::T) where {T} = Nullable{T}(x)
convert(::Type{Nullable   }, x::T) where {T} = Nullable{T}(x)

convert(::Type{Nullable{T}}, ::Void) where {T} = Nullable{T}()
convert(::Type{Nullable   }, ::Void) = Nullable{Union{}}()

promote_rule(::Type{Nullable{S}}, ::Type{T}) where {S,T} = Nullable{promote_type(S, T)}
promote_rule(::Type{Nullable{S}}, ::Type{Nullable{T}}) where {S,T} = Nullable{promote_type(S, T)}
promote_op(op::Any, ::Type{Nullable{S}}, ::Type{Nullable{T}}) where {S,T} = Nullable{promote_op(op, S, T)}
promote_op(op::Type, ::Type{Nullable{S}}, ::Type{Nullable{T}}) where {S,T} = Nullable{promote_op(op, S, T)}

function show(io::IO, x::Nullable)
    if get(io, :compact, false)
        if isnull(x)
            print(io, "#NULL")
        else
            show(io, x.value)
        end
    else
        print(io, "Nullable{")
        showcompact(io, eltype(x))
        print(io, "}(")
        if !isnull(x)
            showcompact(io, x.value)
        end
        print(io, ')')
    end
end

"""
    get(x::Nullable[, y])

Attempt to access the value of `x`. Returns the value if it is present;
otherwise, returns `y` if provided, or throws a `NullException` if not.
"""
@inline function get(x::Nullable{T}, y) where T
    if isbits(T)
        ifelse(isnull(x), y, x.value)
    else
        isnull(x) ? y : x.value
    end
end

get(x::Nullable) = isnull(x) ? throw(NullException()) : x.value

"""
    unsafe_get(x)

Return the value of `x` for [`Nullable`](@ref) `x`; return `x` for
all other `x`.

This method does not check whether or not `x` is null before attempting to
access the value of `x` for `x::Nullable` (hence "unsafe").

```jldoctest
julia> x = Nullable(1)
Nullable{Int64}(1)

julia> unsafe_get(x)
1

julia> x = Nullable{String}()
Nullable{String}()

julia> unsafe_get(x)
ERROR: UndefRefError: access to undefined reference
Stacktrace:
 [1] unsafe_get(::Nullable{String}) at ./nullable.jl:125

julia> x = 1
1

julia> unsafe_get(x)
1
```
"""
unsafe_get(x::Nullable) = x.value
unsafe_get(x) = x

"""
    isnull(x)

Return whether or not `x` is null for [`Nullable`](@ref) `x`; return
`false` for all other `x`.

# Examples

```jldoctest
julia> x = Nullable(1, false)
Nullable{Int64}()

julia> isnull(x)
true

julia> x = Nullable(1, true)
Nullable{Int64}(1)

julia> isnull(x)
false

julia> x = 1
1

julia> isnull(x)
false
```
"""
isnull(x::Nullable) = !x.hasvalue
isnull(x) = false

## Operators

"""
    null_safe_op(f::Any, ::Type, ::Type...)::Bool

Returns whether an operation `f` can safely be applied to any value of the passed type(s).
Returns `false` by default.

Custom types should implement methods for some or all operations `f` when applicable:
returning `true` means that the operation may be called on any bit pattern without
throwing an error (though returning invalid or nonsensical results is not a problem).
In particular, this means that the operation can be applied on the whole domain of the
type *and on uninitialized objects*. As a general rule, these properties are only true for
safe operations on `isbits` types.

Types declared as safe can benefit from higher performance for operations on nullable: by
always computing the result even for null values, a branch is avoided, which helps
vectorization.
"""
null_safe_op(f::Any, ::Type, ::Type...) = false

const NullSafeSignedInts = Union{Type{Int128}, Type{Int16}, Type{Int32},
                                 Type{Int64}, Type{Int8}}
const NullSafeUnsignedInts = Union{Type{Bool}, Type{UInt128}, Type{UInt16},
                                   Type{UInt32}, Type{UInt64}, Type{UInt8}}
const NullSafeInts = Union{NullSafeSignedInts, NullSafeUnsignedInts}
const NullSafeFloats = Union{Type{Float16}, Type{Float32}, Type{Float64}}
const NullSafeTypes = Union{NullSafeInts, NullSafeFloats}
const EqualOrLess = Union{typeof(isequal), typeof(isless)}

null_safe_op(::typeof(identity), ::Type{T}) where {T} = isbits(T)

null_safe_op(f::EqualOrLess, ::NullSafeTypes, ::NullSafeTypes) = true
null_safe_op(f::EqualOrLess, ::Type{Rational{S}}, ::Type{T}) where {S,T} =
    null_safe_op(f, T, S)
# complex numbers can be compared for equality but not in general ordered
null_safe_op(::typeof(isequal), ::Type{Complex{S}}, ::Type{T}) where {S,T} =
    null_safe_op(isequal, T, S)

"""
    isequal(x::Nullable, y::Nullable)

If neither `x` nor `y` is null, compare them according to their values
(i.e. `isequal(get(x), get(y))`). Else, return `true` if both arguments are null,
and `false` if one is null but not the other: nulls are considered equal.
"""
@inline function isequal(x::Nullable{S}, y::Nullable{T}) where {S,T}
    if null_safe_op(isequal, S, T)
        (isnull(x) & isnull(y)) | (!isnull(x) & !isnull(y) & isequal(x.value, y.value))
    else
        (isnull(x) & isnull(y)) || (!isnull(x) & !isnull(y) && isequal(x.value, y.value))
    end
end

isequal(x::Nullable{Union{}}, y::Nullable{Union{}}) = true
isequal(x::Nullable{Union{}}, y::Nullable) = isnull(y)
isequal(x::Nullable, y::Nullable{Union{}}) = isnull(x)

"""
    isless(x::Nullable, y::Nullable)

If neither `x` nor `y` is null, compare them according to their values
(i.e. `isless(get(x), get(y))`). Else, return `true` if only `y` is null, and `false`
otherwise: nulls are always considered greater than non-nulls, but not greater than
another null.
"""
@inline function isless(x::Nullable{S}, y::Nullable{T}) where {S,T}
    # NULL values are sorted last
    if null_safe_op(isless, S, T)
        (!isnull(x) & isnull(y)) | (!isnull(x) & !isnull(y) & isless(x.value, y.value))
    else
        (!isnull(x) & isnull(y)) || (!isnull(x) & !isnull(y) && isless(x.value, y.value))
    end
end

isless(x::Nullable{Union{}}, y::Nullable{Union{}}) = false
isless(x::Nullable{Union{}}, y::Nullable) = false
isless(x::Nullable, y::Nullable{Union{}}) = !isnull(x)

==(x::Nullable, y::Nullable) = throw(NullException())

const nullablehash_seed = UInt === UInt64 ? 0x932e0143e51d0171 : 0xe51d0171

function hash(x::Nullable, h::UInt)
    if isnull(x)
        return h + nullablehash_seed
    else
        return hash(x.value, h + nullablehash_seed)
    end
end

# higher-order functions
"""
    filter(p, x::Nullable)

Return null if either `x` is null or `p(get(x))` is false, and `x` otherwise.
"""
function filter(p, x::Nullable{T}) where T
    if isbits(T)
        val = unsafe_get(x)
        Nullable{T}(val, !isnull(x) && p(val))
    else
        isnull(x) || p(unsafe_get(x)) ? x : Nullable{T}()
    end
end

"""
Return the given type if it is concrete, and `Union{}` otherwise.
"""
nullable_returntype(::Type{T}) where {T} = isleaftype(T) ? T : Union{}

"""
    map(f, x::Nullable)

Return `f` applied to the value of `x` if it has one, as a `Nullable`. If `x`
is null, then return a null value of type `Nullable{S}`. `S` is guaranteed to
be either `Union{}` or a concrete type. Whichever of these is chosen is an
implementation detail, but typically the choice that maximizes performance
would be used. If `x` has a value, then the return type is guaranteed to be of
type `Nullable{typeof(f(x))}`.
"""
function map(f, x::Nullable{T}) where T
    S = promote_op(f, T)
    if isleaftype(S) && null_safe_op(f, T)
        Nullable(f(unsafe_get(x)), !isnull(x))
    else
        if isnull(x)
            Nullable{nullable_returntype(S)}()
        else
            Nullable(f(unsafe_get(x)))
        end
    end
end

# We need the following function and specializations because LLVM cannot
# optimize !any(isnull, t) without further guidance.
hasvalue(x::Nullable) = x.hasvalue
hasvalue(x) = true
all(f::typeof(hasvalue), t::Tuple) = f(t[1]) & all(f, tail(t))
all(f::typeof(hasvalue), t::Tuple{}) = true

# Overloads of null_safe_op
# Unary operators

# Note this list does not include sqrt since it can raise a DomainError
for op in (+, -, abs, abs2)
    null_safe_op(::typeof(op), ::NullSafeTypes) = true
    null_safe_op(::typeof(op), ::Type{Complex{S}}) where {S} = null_safe_op(op, S)
    null_safe_op(::typeof(op), ::Type{Rational{S}}) where {S} = null_safe_op(op, S)
end

null_safe_op(::typeof(~), ::NullSafeInts) = true
null_safe_op(::typeof(!), ::Type{Bool}) = true

# Binary operators

# Note this list does not include ^, ÷ and %
# Operations between signed and unsigned types are not safe: promotion to unsigned
# gives an InexactError for negative numbers
for op in (+, -, *, /, &, |, <<, >>, >>>,
           scalarmin, scalarmax)
    # to fix ambiguities
    null_safe_op(::typeof(op), ::NullSafeFloats, ::NullSafeFloats) = true
    null_safe_op(::typeof(op), ::NullSafeSignedInts, ::NullSafeSignedInts) = true
    null_safe_op(::typeof(op), ::NullSafeUnsignedInts, ::NullSafeUnsignedInts) = true
end
for op in (+, -, *, /)
    null_safe_op(::typeof(op), ::Type{Complex{S}}, ::Type{T}) where {S,T} =
        null_safe_op(op, T, S)
    null_safe_op(::typeof(op), ::Type{Rational{S}}, ::Type{T}) where {S,T} =
        null_safe_op(op, T, S)
end

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
module IfPartialBase

include("./../../constraints/geometric/constants.jl")
using JuMP
using .Constants

export if_pb

function if_pb(m)

    #if (flbz_front_left_bot_z[i] == flbz_front_left_bot_z[j] + hob_height_of_box[j] &&
    #    flbz_front_left_bot_z[i] - hob_height_of_box[j] >= 0 &&
    #    flbx_front_left_bot_x[j] > flbx_front_left_bot_x[i] - box_i_x_axis_cover &&
    #    flbx_front_left_bot_x[j] < flbx_front_left_bot_x[i] + box_i_x_axis_cover &&
    #    flby_front_left_bot_y[j] > flby_front_left_bot_y[i] - box_i_y_axis_cover &&
    #    flby_front_left_bot_y[j] < flby_front_left_bot_y[i] + box_i_y_axis_cover)

    # IF 1
    @variable(m, if11_pbs[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @variable(m, if12_pbs[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @variable(m, if1_pbs[i=1:length(boxes), j=1:length(boxes); i!=j] >= 0, Bin)

    #Hz1 >= x-b+e
    #H(1-z1) >= b-x-e
    @constraint(m, c1_if11_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        H * m[:if11_pbs][i, j] >=
        m[:flbz_front_left_bot_z][i] - (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j]) + 0.001
    )
    @constraint(m, c2_if11_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        H * (1 - m[:if11_pbs][i, j]) >=
        (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j]) - m[:flbz_front_left_bot_z][i] - 0.001
    )

    #Hz2 >= b-x+e
    #H(1-z2) >= x-b-e
    @constraint(m, c1_if12_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        H * m[:if12_pbs][i, j] >=
        (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j]) - m[:flbz_front_left_bot_z][i] + 0.001
    )
    @constraint(m, c2_if12_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        H * (1 - m[:if12_pbs][i, j]) >=
        m[:flbz_front_left_bot_z][i] - (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j]) - 0.001
    )

    #If if11, if12 = 1, then box i and j have z-bot/top value within interval -e <= z <= +e, hence, expression = 1
    @constraint(m, c_if1_pbs[i=1:length(boxes), j=1:length(boxes); i!=j],
        m[:if11_pbs][i, j] + m[:if12_pbs][i, j] >=
        2 * m[:if1_pbs][i, j]
    )

    return m
end

end

ext: jl


text:
importall Bukdu

type WelcomeController <: ApplicationController
    stack::Vector
    WelcomeController() = new([])
end

show(c::WelcomeController) = c.stack
create(c::WelcomeController) = c.stack
secret(c::WelcomeController) = c.stack
prepend(c::WelcomeController) = push!(c.stack, (:p,c[:action]))

function plugins(c::WelcomeController)
   if c[:action] in [show, create, secret]
       plug(prepend, c)
   end
end

function before(c::WelcomeController)
    push!(c.stack, (:b,c[:action]))
end

Router() do
    get("/", WelcomeController, show)
    post("/", WelcomeController, create)
    patch("/", WelcomeController, secret)
end

using Base.Test
conn = (Router)(get, "/")
@test 200 == conn.status
@test [(:p,show),(:b,show)] == conn.resp_body

conn = (Router)(post, "/")
@test 200 == conn.status
@test [(:p,create),(:b,create)] == conn.resp_body

conn = (Router)(patch, "/")
@test 200 == conn.status
@test [(:p,secret),(:b,secret)] == conn.resp_body


text:
# The missing values are `-`. Replace with `?`

dataset_path = joinpath(dirname(@__FILE__), "trains.data.orig")

f = open(dataset_path)
data = readlines(f)
close(f)

f = open("$dataset_path.custom", "w")
for line in data
  line = replace(line, "-", "?")
  write(f, line)
end
close(f)

max_stars_repo_name: CoventryResearch/uci-data

max_issues_repo_name: CoventryResearch/uci-data

max_forks_repo_name: CoventryResearch/uci-data


text:
############################################################################
#
#  Airframe.jl
#
#  Purpose:    Single file that contains all functions needed to simulate a
#              6-DOF Airframe.
#
#  Author:     Edward Daughtery
#
#  Date:       05 August 2016
#
#
############################################################################

# Initialize Position
function Airframe5DOF( time, latitude, longitude, altitude )

    ## Transformation Matrices ECEF, ECI, NED
    Tr_ei = TR_EI(time);
    Tr_ie = Tr_ei';
    Tr_ne = TR_NE(latitude, longitude);

    ## Position Vectors
    rm_ecef = LLH2ECEF(latitude, longitude, altitude);
    rm_eci = Tr_ie*rm_ecef;

    ## Populating Position Structure
    position = WGS_Pos(time, latitude, longitude, altitude, rm_ecef, rm_eci, Tr_ei, Tr_ne);
    return position;
end

############################################################################
# Initialize Velocity
function InitVel( vmag::Float64, gamma::Float64, azimuth::Float64, position::WGS_Pos )

    time = position.time;
    Tr_ei = position.Tr_ei;
    Tr_ie = Tr_ei';
    Tr_ne = position.Tr_ne;
    Tr_en = Tr_ne';

    wb_ei_eci = [0.0;0.0;WEII3];

    vm_ned = [
        vmag*cos(gamma)*cos(azimuth)
        vmag*cos(gamma)*sin(azimuth)
        -vmag*sin(gamma)
    ];

    vm_ecef = Tr_en*vm_ned;
    vm_eci = Tr_ie*vm_ecef + cross(wb_ei_eci, position.rm_eci);

    velocity = Velocity_States( time, vmag, gamma, azimuth, vm_ned, vm_ecef, vm_eci );
    return velocity;
end

############################################################################
# Initialize Euler
function InitEuler( Roll::Float64, Pitch::Float64, Heading::Float64 )
    Tr_bn = TR_BN( Roll, Pitch, Heading);

    euler = Euler_Angles(Roll, Pitch, Heading, Tr_bn);

    return euler;
end

############################################################################
# Initialize Rates
function InitAngleRates(wb_bn_body::Array{Float64,1}, position::WGS_Pos, euler::Euler_Angles)

    ## Transformation Matrices
    Tr_ei = position.Tr_ei;
    Tr_ie = Tr_ei';
    Tr_ne = position.Tr_ne;
    Tr_en = Tr_ne';
    Tr_bn = euler.Tr_bn;
    Tr_nb = Tr_bn';
    Tr_bi = Tr_bn*Tr_ne*Tr_ei;
    Tr_ib = Tr_bi';

    ## Earth Rotation Rate
    wb_ne_body = [0.0;0.0;0.0];
    wb_ne_eci = Tr_ib*wb_ne_body;
    wb_ei_eci = [0.0;0.0;WEII3];
    wb_ei_ecef = Tr_ei*wb_ei_eci;
    wb_bi_body = wb_bn_body + wb_ne_body + Tr_bi*wb_ei_eci;

    ## Transport Rate
    angle_rates = Angle_Rates( wb_bi_body, wb_bn_body, wb_ei_eci, wb_ei_ecef, wb_ne_eci )
    return angle_rates
end

############################################################################
# Initialize
function Airframe5_Initialize( tsim::Float64,
                     time::Float64, dt::Float64,
                     position::WGS_Pos, velocity::Velocity_States,
                     euler::Euler_Angles, angle_rates::Angle_Rates, mass::Float64)

    ## Time Computations
    dt_eci = time - tsim;

    ## Transformation Matrices
    Tr_ei = position.Tr_ei;
    Tr_ie = Tr_ei';
    Tr_ne = position.Tr_ne;
    Tr_en = Tr_ne';
    Tr_bn = euler.Tr_bn;
    Tr_nb = Tr_bn';
    Tr_bi = Tr_bn*Tr_ne*Tr_ei;
    Tr_ib = Tr_bi';

    quat = QuatInit( Tr_bi );
    quat_bi = Quaternions(quat[1], quat[2], quat[3], quat[4] );

    zz = [0.0;0.0;0.0];

    xo = [ position.rm_eci; velocity.vm_eci];

    ## Quad Rotor Intialization
    airframe_init = Aero5DOF(time, dt, position, velocity, euler, angle_rates, zz,  mass, quat_bi, xo);

    return airframe_init
end

############################################################################
# Update
function Airframe5_Update( airframe::Aero5DOF, ftot_body::Array{Float64,1}
			   ,alpha::Float64, alphadot::Float64, phi::Float64, phidot::Float64)
    dt = airframe.dt;
    airframe_out = airframe;

    xo = airframe.xo;

    rm_eci = airframe.position.rm_eci;
    vm_eci = airframe.velocity.vm_eci;

    ## Time
    time = airframe.time;
    wb_ei_eci = [0.0;0.0;WEII3];
    wb_ei_ecef = wb_ei_eci;

    ## Transformation Matrices
    Tr_ei = TR_EI(time);
    Tr_ie = Tr_ei';
    Tr_ne = TR_NE(airframe.position.latitude, airframe.position.longitude)
    Tr_en = Tr_ne';
    Tr_bn = airframe.euler.Tr_bn;
    Tr_nb = Tr_bn';

    ## Mass Properties
    mass = airframe.mass;

    ## Airframe Velocity
    vel = airframe.velocity;
    vm_eci = airframe.velocity.vm_eci;
    vm_ecef = Tr_ei*vm_eci - cross(wb_ei_ecef, airframe.position.rm_ecef);
    vm_ned = Tr_ne*vm_ecef;
    vm_body = Tr_bn*vm_ned;
    vmag = vel.vmag;
    gamma = vel.gamma;
    azimuth = vel.azimuth;

    ## Finishing Computation of Transformation Matrices
    Tr_vn = TR_VN(gamma,azimuth);
    Tr_nv = Tr_vn';
    Tr_bv = TR_BV(alpha, phi);
    Tr_vb = Tr_bv';
    Tr_bi = Tr_bv*Tr_vn*Tr_ne*Tr_ei;
    Tr_ib = Tr_bi';
    Tr_bn = Tr_bv*Tr_vn;
    Tr_nb = Tr_bn';

    ## Airframe Position
    pos = airframe.position;

    ## Quaternions
    quat_bi = airframe.quat_bi;
    quat_bi_vec = [quat_bi.q0; quat_bi.q1; quat_bi.q2; quat_bi.q3];

    ## Gravity
    grav_ecef = WGS84_GRAVITY(airframe.position.rm_ecef);
    grav_eci = Tr_ie*grav_ecef - cross( wb_ei_eci, cross(wb_ei_eci, airframe.position.rm_eci) );
    grav_ecef = Tr_ei*grav_eci;
    grav_body = Tr_bi*grav_eci;
    grav_ned = Tr_nb*grav_body;
    fgrav_body = mass*grav_body;

    ## Total Forces in Body Frame
    acc_meas_body = ftot_body/mass;
    ftot_eci = Tr_ib*ftot_body;
    airframe_out.accel_meas = acc_meas_body;

    ## Total Acceleration
    acc_tot_body = acc_meas_body + grav_body;
    acc_tot_ned = Tr_nb*acc_tot_body;
    acc_tot_vel = Tr_vn*acc_tot_ned;

    ## Gamma Dot and Azimuth Dot
    gamdot  = -acc_tot_vel[3]/vmag;
    azdot   = acc_tot_vel[2]/(vmag*cos(gamma));

    ## Rotation Rate of Velocity Frame wrt to NED frame
    wvn_vel = zeros(3);
    wvn_vel[1] = -azdot*sin(gamma);
    wvn_vel[2] = gamdot;
    wvn_vel[3] = azdot*cos(gamma);
    wvn_body = Tr_bv*wvn_vel;

    ## Transport Rate
    #wne_ned[1] =  vm_ned[2]/RC.Lon;
    #wne_ned[2] = -vm_ned[1]/RC.Lat;
    #wne_ned[3] = -( vm_ned[2]*wgs_Pos.Latitude.Tan() )/RC.Lon;
    wne_ned = zeros(3);
    wne_body = Tr_bn*wne_ned;
    wne_eci = Tr_ib*wne_body;

    ## Velocity Angle rates
    wbv_body = zeros(3);
    wbv_body[1] = phidot*cos(alpha);
    wbv_body[2] = alphadot;
    wbv_body[3] = phidot*sin(alpha);

    ## Converting Earth Rotation Rate Vector to Body Frame
    wei_body = Tr_bi*wb_ei_eci;

    ## Computing Rate of body frame wrt inertial frame
    wbi_body = wbv_body + wvn_body + wne_body + wei_body;
    wbn_body = wbv_body + wvn_body;
    pp = wbn_body[1];
    qq = wbn_body[2];
    rr = wbn_body[3];

    ## Computing State Vector
    xo = [
            airframe.position.rm_eci;
            airframe.velocity.vm_eci;
         ];

    ## Computing Derivatives
    dx_pos = vm_eci;
    dx_vel = Tr_ib*(1.0/mass*( ftot_body + fgrav_body - cross( wbi_body, mass*vm_body ) ) );

    ## Derivative Vector
    dx = [dx_pos; dx_vel];

    ## Euler Integration
    xo = xo + dt*dx;
    time = time + dt;

    ## Updating States
    rm_eci = xo[1:3];
    vm_eci = xo[4:6];

    ## Position Update
    rm_ecef = Tr_ei*rm_eci;
    llh = ECEF2LLH(rm_ecef);
    lat = llh[1];
    lon = llh[2];
    alt = llh[3];
    Tr_ne = TR_NE(lat,lon);
    Tr_ei = TR_EI(time);
    Tr_ie = Tr_ei';
    pos = WGS_Pos(time, lat, lon, alt, rm_ecef, rm_eci, Tr_ei, Tr_ne);

    ## Velocity Update
    vm_ecef = Tr_ei*vm_eci - cross(wb_ei_ecef, rm_ecef);
    vm_ned = Tr_ne*vm_ecef;
    vmag = norm(vm_ned);
    vn = vm_ned[1];
    ve = vm_ned[2];
    vd = vm_ned[3];
    gamma = atan( -vd, sqrt( vn*vn + ve*ve ) );
    azimuth = atan( ve, vn );
    Tr_vn = TR_VN(gamma,azimuth);
    vel = Velocity_States(time, vmag, gamma, azimuth, vm_ned, vm_ecef, vm_eci);

    ## Updating Transformations
    Tr_nv = Tr_vn';
    Tr_en = Tr_ne';
    Tr_bn = Tr_bv*Tr_vn;
    Tr_nb = Tr_bn';
    Tr_eb = Tr_en*Tr_nb;
    Tr_be = Tr_eb';
    Tr_bi = Tr_be*Tr_ei;
    Tr_ib = Tr_bi';

    ## Euler Angle Updates
    eulervec = ComputeEuler( Tr_bn )
    roll = eulervec[1];
    pitch = eulervec[2];
    heading = eulervec[3];
    euler = Euler_Angles(roll,  pitch, heading, Tr_bn);

    ## Updating Quaternions
    quat_bi_body = QuatInit(Tr_bi);
    quat = Quaternions(quat_bi_body[1],
                       quat_bi_body[2],
                       quat_bi_body[3],
                       quat_bi_body[4] );

    ## Update Angle Rates
    angle_rates = Angle_Rates(wbi_body, wbn_body, wb_ei_eci, wb_ei_ecef, wne_eci)

    ## Updating the airframe structure
    airframe_out.time = time;
    airframe_out.position = pos;
    airframe_out.velocity = vel;
    airframe_out.euler = euler;
    airframe_out.angle_rates = angle_rates;
    airframe_out.quat_bi = quat;
    airframe_out.xo = xo;
    airframe_out.accel_meas = acc_meas_body;

    return airframe_out;
end


text:
function Base.show(io::IO, x::EXPR, offset = 0, d = 0, er = false)
    T = headof(x)
    c =  T === :errortoken || er ? :red : :normal
    # Print span as 1-based range of the source string. This presentation is
    # simple to understand when strings are presented to CSTParser.parse().
    print(io, lpad(offset + 1, 3), ":", rpad(offset + x.fullspan, 3), " ")
    if isidentifier(x)
        printstyled(io, " "^d, headof(x) == :NONSTDIDENTIFIER ? valof(x.args[2]) : valof(x), color = :yellow)
        x.meta !== nothing && show(io, x.meta)
        print(io)
    elseif isoperator(x)
        printstyled(io, " "^d, "OP: ", valof(x), color = c)
    elseif iskeyword(x)
        printstyled(io, " "^d, headof(x), color = :magenta)
    elseif ispunctuation(x)
        printstyled(io, " "^d, punctuationprinting[headof(x)], color = c)
    elseif isliteral(x)
        printstyled(io, " "^d, "$(headof(x)): ", valof(x) === nothing ? "nothing" : valof(x), color = c)
    else
        printstyled(io, " "^d, T, color=c)
        if x.meta !== nothing
            print(io, "( ")
            show(io, x.meta)
            print(io, ")")
        end
        x.args === nothing && return
        for a in x.args
            println(io)
            show(io, a, offset, d + 1, er)
            offset += a.fullspan
        end
    end
end

const punctuationprinting = Dict(
    :COMMA => ",",
    :LPAREN => "(",
    :RPAREN => ")",
    :LSQUARE => "[",
    :RSQUARE => "]",
    :LBRACE => "{",
    :RBRACE => "}",
    :ATSIGN => "@",
    :DOT => "."
)

struct CSTInfiniteLoop <: Exception
    msg::AbstractString
end

function Base.showerror(io::IO, ex::CSTInfiniteLoop)
    print(io, ex.msg)
end


text:
const NO_ENTITY_SUBSTITUTION = Set([:script, :style])
const EMPTY_TAGS = Set([
    :area,
    :base,
    :basefont,
    :bgsound,
    :br,
    :command,
    :col,
    :embed,
    Symbol("event-source"),
    :frame,
    :hr,
    :image,
    :img,
    :input,
    :keygen,
    :link,
    :menuitem,
    :meta,
    :param,
    :source,
    :spacer,
    :track,
    :wbr
])
const RELEVANT_WHITESPACE = Set([:pre, :textarea, :script, :style])

function substitute_text_entities(str)
    str = replace(str, "&" => "&")
    str = replace(str, "<" => "<")
    str = replace(str, ">" => ">")

    return str
end

function substitute_attribute_entities(str)
    str = substitute_text_entities(str)
    str = replace(str, "\""=> """)
    str = replace(str, "'" => "'")

    return str
end

function Base.print(io::IO, elem::HTMLElement{T}; pretty = false, depth = 0, substitution = true) where {T}
    empty_tag = T in EMPTY_TAGS
    ws_relevant = T in RELEVANT_WHITESPACE
    has_children = !isempty(elem.children)

    pretty_children = pretty && !ws_relevant

    pretty && print(io, ' '^(2*depth))
    print(io, '<', T)
    for (name, value) in sort(collect(elem.attributes), by = first)
        print(io, ' ', name, "=\"", substitute_attribute_entities(value), '"')
    end
    if empty_tag
        print(io, '/')
    end
    print(io, '>')
    pretty_children && has_children && print(io, '\n')

    if !empty_tag
        for child in elem.children
            print(io, child; pretty = pretty_children, depth = depth + 1, substitution = substitution && !in(T, NO_ENTITY_SUBSTITUTION))
        end

        pretty && has_children && print(io, ' '^(2*depth))
        print(io, "</", T, '>')
    end
    pretty && print(io, '\n')

    return nothing
end

function Base.print(io::IO, node::HTMLText; pretty = false, depth = 0, substitution = true)
    substitutor = substitution ? substitute_text_entities : identity

    if !pretty
        print(io, substitutor(node.text))
        return nothing
    end

    for line in strip.(split(node.text, '\n'))
        isempty(line) && continue
        print(io, ' '^(2*depth), substitutor(line), '\n')
    end
end

function Base.print(io::IO, doc::HTMLDocument; pretty = false)
    write(io, "<!DOCTYPE ", doc.doctype, ">")
    Base.print(io, doc.root, pretty = pretty)
end

prettyprint(io::IO, doc::HTMLDocument) = Base.print(io, doc, pretty = true)
prettyprint(doc::HTMLDocument) = prettyprint(stdout, doc)
prettyprint(io::IO, elem::HTMLElement) = print(io, elem, pretty = true)
prettyprint(elem::HTMLElement) = print(stdout, elem, pretty = true)

function Base.show(io::IO, elem::HTMLElement)
    write(io, summary(elem), ":")
    if get(io, :compact, false)
        return
    elseif get(io, :limit, false)
        buf = IOBuffer()
        print(buf, elem, pretty = true)
        for (i, line) in enumerate(split(String(take!(buf)), '\n'))
            if i > 20
                println(io, "...")
                return
            end

            println(io, line)
        end
    else
        print(io, elem, pretty=true)
    end
end

function Base.show(io::IO, t::HTMLText)
    write(io,"HTML Text: `", t.text, '`')
end

function Base.show(io::IO, doc::HTMLDocument)
    write(io, "HTML Document:\n")
    write(io, "<!DOCTYPE ", doc.doctype, ">\n")
    Base.show(io, doc.root)
end


text:
function mccormick_envelope(model::Model, w::Variable, u::Variable, v::Variable, urange::ClosedInterval, vrange::ClosedInterval)
    # convex approximation of w == u * v on u ∈ [umin, umax], v ∈ [vmin, vmax]
    umin, umax = minimum(urange), maximum(urange)
    vmin, vmax = minimum(vrange), maximum(vrange)
    vertices = [umin vmin umin * vmin; umax vmax umax * vmax; umax vmin umax * vmin; umin vmax umin * vmax]
    SimpleHRepresentation(hrep(polyhedron(SimpleVRepresentation(vertices), CDDLibrary())))
end

function piecewise_mccormick_envelope_constraints(model::Model, u::Variable, v::Variable, w::Variable,
    urange::Range, vrange::ClosedInterval)
    vars = [u; v; w]
    if length(urange) == 2
        envelope = mccormick_envelope(model, u, v, w, first(urange)..last(urange), vrange)
        polyhedron_constraints(model, envelope, vars)
        nothing
    else
        N = length(urange) - 1
        envelopes = Vector{SimpleHRepresentation}(N)
        for i = 1 : N
            envelopes[i] = mccormick_envelope(model, u, v, w, urange[i]..urange[i + 1], vrange)
        end
        hull_reformulation(model, envelopes, vars)
        nothing
    end
end


text:
export BrickletPiezoBuzzerIdentity
struct BrickletPiezoBuzzerIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletPiezoBuzzer
"""
Creates 1kHz beep
"""
mutable struct BrickletPiezoBuzzer <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletPiezoBuzzer(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_piezo_buzzer")
        deviceInternal = package.BrickletPiezoBuzzer(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export beep
"""
    $(SIGNATURES)

Beeps for the given duration.
"""
function beep(device::BrickletPiezoBuzzer, duration)
    device.deviceInternal.beep(duration)
end

export morse_code
"""
    $(SIGNATURES)

Sets morse code that will be played by the piezo buzzer. The morse code
is given as a string consisting of "." (dot), "-" (minus) and " " (space)
for *dits*, *dahs* and *pauses*. Every other character is ignored.

For example: If you set the string "...---...", the piezo buzzer will beep
nine times with the durations "short short short long long long short
short short".
"""
function morse_code(device::BrickletPiezoBuzzer, morse)
    device.deviceInternal.morse_code(morse)
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletPiezoBuzzer)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletPiezoBuzzer, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using Test
using PrettyTables

data = Any[1    false      1.0     0x01 ;
           2     true      2.0     0x02 ;
           3    false      3.0     0x03 ;
           4     true      4.0     0x04 ;
           5    false      5.0     0x05 ;
           6     true      6.0     0x06 ;]

# Default
# ==============================================================================

@testset "Default" begin

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = sprint(pretty_table, data)
    @test result == expected
end

# Alignments
# ==============================================================================

@testset "Alignments" begin
    # Left
    # ==========================================================================
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│ 1      │ false  │ 1.0    │ 1      │
│ 2      │ true   │ 2.0    │ 2      │
│ 3      │ false  │ 3.0    │ 3      │
│ 4      │ true   │ 4.0    │ 4      │
│ 5      │ false  │ 5.0    │ 5      │
│ 6      │ true   │ 6.0    │ 6      │
└────────┴────────┴────────┴────────┘
"""
    result = sprint((io, data)->pretty_table(io, data; alignment = :l), data)
    @test result == expected

    # Center
    # ==========================================================================
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│   1    │ false  │  1.0   │   1    │
│   2    │  true  │  2.0   │   2    │
│   3    │ false  │  3.0   │   3    │
│   4    │  true  │  4.0   │   4    │
│   5    │ false  │  5.0   │   5    │
│   6    │  true  │  6.0   │   6    │
└────────┴────────┴────────┴────────┘
"""
    result = sprint((io, data)->pretty_table(io, data; alignment = :c), data)
    @test result == expected

    # Per column configuration
    # ==========================================================================

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│ 1      │  false │  1.0   │      1 │
│ 2      │   true │  2.0   │      2 │
│ 3      │  false │  3.0   │      3 │
│ 4      │   true │  4.0   │      4 │
│ 5      │  false │  5.0   │      5 │
│ 6      │   true │  6.0   │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = sprint((io, data)->pretty_table(io, data;
                                             alignment = [:l,:r,:c,:r]),
                    data)
    @test result == expected
end

# Formatter
# ==============================================================================

@testset "Formatter" begin
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │      1 │      1 │
│      0 │   true │      0 │      0 │
│      3 │  false │      3 │      3 │
│      0 │   true │      0 │      0 │
│      5 │  false │      5 │      5 │
│      0 │   true │      0 │      0 │
└────────┴────────┴────────┴────────┘
"""
    formatter = Dict(0 => (v,i) -> isodd(i) ? i : 0,
                     2 => (v,i) -> v)
    result = sprint((io, data)->pretty_table(io, data;
                                             formatter = formatter), data)
    @test result == expected
end

# Show row number
# ==============================================================================

@testset "Show row number" begin
    expected = """
┌─────┬────────┬────────┬────────┬────────┐
│ Row │ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├─────┼────────┼────────┼────────┼────────┤
│   1 │ 1      │  false │  1.0   │      1 │
│   2 │ 2      │   true │  2.0   │      2 │
│   3 │ 3      │  false │  3.0   │      3 │
│   4 │ 4      │   true │  4.0   │      4 │
│   5 │ 5      │  false │  5.0   │      5 │
│   6 │ 6      │   true │  6.0   │      6 │
└─────┴────────┴────────┴────────┴────────┘
"""
    result = sprint((io, data)->pretty_table(io, data;
                                             alignment       = [:l,:r,:c,:r],
                                             show_row_number = true),
                    data)
    @test result == expected
end

# Pre-defined formats
# ==============================================================================

@testset "Pre-defined formats" begin

    # ascii_dots
    # ==========================================================================
    expected = """
.....................................
: Col. 1 : Col. 2 : Col. 3 : Col. 4 :
:........:........:........:........:
:      1 :  false :    1.0 :      1 :
:      2 :   true :    2.0 :      2 :
:      3 :  false :    3.0 :      3 :
:      4 :   true :    4.0 :      4 :
:      5 :  false :    5.0 :      5 :
:      6 :   true :    6.0 :      6 :
:........:........:........:........:
"""
    result = sprint(pretty_table, data, ascii_dots)
    @test result == expected

    # ascii_rounded
    # ==========================================================================
    expected = """
.--------.--------.--------.--------.
| Col. 1 | Col. 2 | Col. 3 | Col. 4 |
:--------+--------+--------+--------:
|      1 |  false |    1.0 |      1 |
|      2 |   true |    2.0 |      2 |
|      3 |  false |    3.0 |      3 |
|      4 |   true |    4.0 |      4 |
|      5 |  false |    5.0 |      5 |
|      6 |   true |    6.0 |      6 |
'--------'--------'--------'--------'
"""
    result = sprint(pretty_table, data, ascii_rounded)
    @test result == expected

    # borderless
    # ==========================================================================
    expected = """
  Col. 1   Col. 2   Col. 3   Col. 4

       1    false      1.0        1
       2     true      2.0        2
       3    false      3.0        3
       4     true      4.0        4
       5    false      5.0        5
       6     true      6.0        6
"""
    result = sprint(pretty_table, data, borderless)
    @test result == expected

    # compact
    # ==========================================================================
    expected = """
 -------- -------- -------- --------
  Col. 1   Col. 2   Col. 3   Col. 4
 -------- -------- -------- --------
       1    false      1.0        1
       2     true      2.0        2
       3    false      3.0        3
       4     true      4.0        4
       5    false      5.0        5
       6     true      6.0        6
 -------- -------- -------- --------
"""
    result = sprint(pretty_table, data, compact)
    @test result == expected

    # markdown
    # ==========================================================================
    expected = """
| Col. 1 | Col. 2 | Col. 3 | Col. 4 |
|--------|--------|--------|--------|
|      1 |  false |    1.0 |      1 |
|      2 |   true |    2.0 |      2 |
|      3 |  false |    3.0 |      3 |
|      4 |   true |    4.0 |      4 |
|      5 |  false |    5.0 |      5 |
|      6 |   true |    6.0 |      6 |
"""
    result = sprint(pretty_table, data, markdown)
    @test result == expected

    # mysql
    # ==========================================================================
    expected = """
+--------+--------+--------+--------+
| Col. 1 | Col. 2 | Col. 3 | Col. 4 |
+--------+--------+--------+--------+
|      1 |  false |    1.0 |      1 |
|      2 |   true |    2.0 |      2 |
|      3 |  false |    3.0 |      3 |
|      4 |   true |    4.0 |      4 |
|      5 |  false |    5.0 |      5 |
|      6 |   true |    6.0 |      6 |
+--------+--------+--------+--------+
"""
    result = sprint(pretty_table, data, mysql)
    @test result == expected

    # simple
    # ==========================================================================
    expected = """
========= ======== ======== =========
  Col. 1   Col. 2   Col. 3   Col. 4
========= ======== ======== =========
       1    false      1.0        1
       2     true      2.0        2
       3    false      3.0        3
       4     true      4.0        4
       5    false      5.0        5
       6     true      6.0        6
========= ======== ======== =========
"""
    result = sprint(pretty_table, data, simple)
    @test result == expected

    # unicode_rounded
    # ==========================================================================
    expected = """
╭────────┬────────┬────────┬────────╮
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
╰────────┴────────┴────────┴────────╯
"""
    result = sprint(pretty_table, data, unicode_rounded)
    @test result == expected
end

# Pre-defined formatters
# ==============================================================================

@testset "Pre-defined formatters" begin

    # ft_round
    # ==========================================================================
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│    1.0 │    0.0 │    1.0 │    1.0 │
│    2.0 │    1.0 │    2.0 │    2.0 │
│    3.0 │    0.0 │    3.0 │    3.0 │
│    4.0 │    1.0 │    4.0 │    4.0 │
│    5.0 │    0.0 │    5.0 │    5.0 │
│    6.0 │    1.0 │    6.0 │    6.0 │
└────────┴────────┴────────┴────────┘
"""
    result = sprint((io,data)->pretty_table(io,data;
                                            formatter = ft_round(1)),
                    data)
    @test result == expected

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│    1.0 │  false │    1.0 │      1 │
│    2.0 │   true │    2.0 │      2 │
│    3.0 │  false │    3.0 │      3 │
│    4.0 │   true │    4.0 │      4 │
│    5.0 │  false │    5.0 │      5 │
│    6.0 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = sprint((io,data)->pretty_table(io,data;
                                            formatter = ft_round(1,[3,1])),
                    data)
    @test result == expected

    # ft_printf
    # ==========================================================================

    expected = """
┌──────────┬──────────┬──────────┬──────────┐
│   Col. 1 │   Col. 2 │   Col. 3 │   Col. 4 │
├──────────┼──────────┼──────────┼──────────┤
│    1.000 │    0.000 │    1.000 │    1.000 │
│    2.000 │    1.000 │    2.000 │    2.000 │
│    3.000 │    0.000 │    3.000 │    3.000 │
│    4.000 │    1.000 │    4.000 │    4.000 │
│    5.000 │    0.000 │    5.000 │    5.000 │
│    6.000 │    1.000 │    6.000 │    6.000 │
└──────────┴──────────┴──────────┴──────────┘
"""
    result = sprint((io,data)->pretty_table(io,data;
                                            formatter = ft_printf("%8.3f")),
                    data)
    @test result == expected

    expected = """
┌──────────┬────────┬────────┬──────────┐
│   Col. 1 │ Col. 2 │ Col. 3 │   Col. 4 │
├──────────┼────────┼────────┼──────────┤
│    1.000 │  false │    1.0 │    1.000 │
│    2.000 │   true │    2.0 │    2.000 │
│    3.000 │  false │    3.0 │    3.000 │
│    4.000 │   true │    4.0 │    4.000 │
│    5.000 │  false │    5.0 │    5.000 │
│    6.000 │   true │    6.0 │    6.000 │
└──────────┴────────┴────────┴──────────┘
"""
    result = sprint((io,data)->
                    pretty_table(io,data; formatter = ft_printf("%8.3f",[1,4])),
                    data)
    @test result == expected

    expected = """
┌──────────┬────────┬────────┬──────────┐
│   Col. 1 │ Col. 2 │ Col. 3 │   Col. 4 │
├──────────┼────────┼────────┼──────────┤
│     1.00 │  false │    1.0 │   1.0000 │
│     2.00 │   true │    2.0 │   2.0000 │
│     3.00 │  false │    3.0 │   3.0000 │
│     4.00 │   true │    4.0 │   4.0000 │
│     5.00 │  false │    5.0 │   5.0000 │
│     6.00 │   true │    6.0 │   6.0000 │
└──────────┴────────┴────────┴──────────┘
"""
    result = sprint((io,data)->
                    pretty_table(io,data;
                                 formatter = ft_printf(["%8.2f","%8.4f"],[1,4])),
                    data)
    @test result == expected
end

# Sub-headers
# ==============================================================================

@testset "Sub-headers" begin
    expected = """
┌───┬───────┬─────┬───┐
│ 1 │     2 │   3 │ 4 │
├───┼───────┼─────┼───┤
│ 1 │ false │ 1.0 │ 1 │
│ 2 │  true │ 2.0 │ 2 │
│ 3 │ false │ 3.0 │ 3 │
│ 4 │  true │ 4.0 │ 4 │
│ 5 │ false │ 5.0 │ 5 │
│ 6 │  true │ 6.0 │ 6 │
└───┴───────┴─────┴───┘
"""
    result = sprint(pretty_table, data, [1;2;3;4])
    @test result == expected

    result = sprint(pretty_table, data, [1 2 3 4])
    @test result == expected

    expected = """
┌───┬───────┬─────┬───┐
│ 1 │     2 │   3 │ 4 │
│ 5 │     6 │   7 │ 8 │
├───┼───────┼─────┼───┤
│ 1 │ false │ 1.0 │ 1 │
│ 2 │  true │ 2.0 │ 2 │
│ 3 │ false │ 3.0 │ 3 │
│ 4 │  true │ 4.0 │ 4 │
│ 5 │ false │ 5.0 │ 5 │
│ 6 │  true │ 6.0 │ 6 │
└───┴───────┴─────┴───┘
"""
    result = sprint(pretty_table, data, [1 2 3 4; 5 6 7 8])
    @test result == expected
end

# Horizontal lines
# ==============================================================================

@testset "Horizontal lines" begin
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
├────────┼────────┼────────┼────────┤
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
├────────┼────────┼────────┼────────┤
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = sprint((io,data)->
                    pretty_table(io,data; hlines = findall(x->x == true, data[:,2])),
                    data)
    @test result == expected
end

# Line breaks inside cells
# ==============================================================================

@testset "Line breaks inside cells" begin
    data = ["This line contains\nthe velocity [m/s]" 10.0;
            "This line contains\nthe acceleration [m/s^2]" 1.0;
            "This line contains\nthe time from the\nbeginning of the simulation" 10;]

    header = ["Information", "Value"]

    expected = """
┌─────────────────────────────┬───────┐
│                 Information │ Value │
├─────────────────────────────┼───────┤
│          This line contains │  10.0 │
│          the velocity [m/s] │       │
│          This line contains │   1.0 │
│    the acceleration [m/s^2] │       │
│          This line contains │    10 │
│           the time from the │       │
│ beginning of the simulation │       │
└─────────────────────────────┴───────┘
"""

    result = sprint((io,data)->pretty_table(io,data,header; linebreaks = true),
                    data)
    @test result == expected

    expected = """
┌─────────────────────────────┬───────┐
│         Information         │ Value │
├─────────────────────────────┼───────┤
│     This line contains      │ 10.0  │
│     the velocity [m/s]      │       │
│     This line contains      │  1.0  │
│  the acceleration [m/s^2]   │       │
│     This line contains      │  10   │
│      the time from the      │       │
│ beginning of the simulation │       │
└─────────────────────────────┴───────┘
"""

    result = sprint((io,data)->pretty_table(io,data,header;
                                            alignment = :c,
                                            linebreaks = true),
                    data)
    @test result == expected

    expected = """
┌─────────────────────────────┬───────┐
│ Information                 │ Value │
├─────────────────────────────┼───────┤
│ This line contains          │ 10.0  │
│ the velocity [m/s]          │       │
│ This line contains          │ 1.0   │
│ the acceleration [m/s^2]    │       │
│ This line contains          │ 10    │
│ the time from the           │       │
│ beginning of the simulation │       │
└─────────────────────────────┴───────┘
"""

    result = sprint((io,data)->pretty_table(io,data,header;
                                            alignment = :l,
                                            linebreaks = true),
                    data)
    @test result == expected

    expected = """
┌────────────────────────────────────────────────────────────────────┬───────┐
│                                                        Information │ Value │
├────────────────────────────────────────────────────────────────────┼───────┤
│                             This line contains\\nthe velocity [m/s] │  10.0 │
│                       This line contains\\nthe acceleration [m/s^2] │   1.0 │
│ This line contains\\nthe time from the\\nbeginning of the simulation │    10 │
└────────────────────────────────────────────────────────────────────┴───────┘
"""

    result = sprint(pretty_table, data, header)
    @test result == expected
end

# Line breaks inside cells
# ==============================================================================

@testset "Hide header" begin

    expected = """
┌───┬───────┬─────┬───┐
│ 1 │ false │ 1.0 │ 1 │
│ 2 │  true │ 2.0 │ 2 │
│ 3 │ false │ 3.0 │ 3 │
│ 4 │  true │ 4.0 │ 4 │
│ 5 │ false │ 5.0 │ 5 │
│ 6 │  true │ 6.0 │ 6 │
└───┴───────┴─────┴───┘
"""

    result = sprint((io, data)->pretty_table(io, data; noheader = true), data)
    @test result == expected

    result = sprint((io, data)->pretty_table(io, data, [1 2]; noheader = true),
                    data)
    @test result == expected

    expected = """
┌───┬───────┬─────┬───┐
│ 1 │ false │ 1.0 │ 1 │
│ 2 │ true  │ 2.0 │ 2 │
│ 3 │ false │ 3.0 │ 3 │
│ 4 │ true  │ 4.0 │ 4 │
│ 5 │ false │ 5.0 │ 5 │
│ 6 │ true  │ 6.0 │ 6 │
└───┴───────┴─────┴───┘
"""
    result = sprint((io, data)->pretty_table(io, data, [1 2]; alignment = :l,
                                             noheader = true), data)
    @test result == expected
end

# Issue #4
# ==============================================================================

@testset "Issue #4 - Header with escape sequences" begin

    expected = """
┌───┬───────┬─────┬───┐
│ 1 │   2\\n │   3 │ 4 │
├───┼───────┼─────┼───┤
│ 1 │ false │ 1.0 │ 1 │
│ 2 │  true │ 2.0 │ 2 │
│ 3 │ false │ 3.0 │ 3 │
│ 4 │  true │ 4.0 │ 4 │
│ 5 │ false │ 5.0 │ 5 │
│ 6 │  true │ 6.0 │ 6 │
└───┴───────┴─────┴───┘
"""

    result = sprint(pretty_table, data, ["1" "2\n" "3" "4"])
    @test result == expected
end


text:
#i64o-typedefs.jl

#a helper function which generates the appropriate number of cells for any given FSS
__cell_length(FSS) = 1 << (FSS - 6)

# a helper function that checks arraynum inputs.
function __check_ArrayNum(FSS, a::Array{UInt64,1})
  FSS < 7 && throw(ArgumentError("invalid FSS == $FSS < 7"))
  _al = __cell_length(FSS)
  length(a) < _al && throw(ArgumentError("invalid array length, should be at least $_al > $(length(a))"))
end

doc"""
`Unums.ArrayNum` is a variadic type which maps an `FSS` variable to an `Int64`
array of a size corresponding to `FSS`.
"""
@dev_check type ArrayNum{FSS}
  a::Array{UInt64,1}
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
"""
    Poles(z) = [x->1/(x - z) for z in z]
"""
struct Poles{Z} <: AbstractBasis
    poles::Z
end

Base.length(B::Poles) = length(B.poles)

evaluate_basis(B::Poles, i, x) = inv(x - B.poles[i]*I)


text:
## Allows to iterate over kernels
Base.length(::Kernel) = 1
Base.iterate(k::Kernel) = (k,nothing)
Base.iterate(k::Kernel, ::Any) = nothing

# default fallback for evaluating a kernel with two arguments (such as vectors etc)
kappa(κ::Kernel, x, y) = kappa(κ, evaluate(metric(κ), x, y))
kappa(κ::TransformedKernel, x, y) = kappa(kernel(κ), apply(κ.transform,x), apply(κ.transform,y))
kappa(κ::TransformedKernel{<:BaseKernel,<:ScaleTransform}, x, y) = kappa(κ, _scale(κ.transform, metric(κ), x, y))
_scale(t::ScaleTransform, metric::Euclidean, x, y) =  first(t.s) * evaluate(metric, x, y)
_scale(t::ScaleTransform, metric::Union{SqEuclidean,DotProduct}, x, y) =  first(t.s)^2 * evaluate(metric, x, y)
_scale(t::ScaleTransform, metric, x, y) = evaluate(metric, apply(t, x), apply(t, y))

printshifted(io::IO,κ::Kernel,shift::Int) = print(io,"$κ")
Base.show(io::IO,κ::Kernel) = print(io,nameof(typeof(κ)))

### Syntactic sugar for creating matrices and using kernel functions
for k in subtypes(BaseKernel)
    @eval begin
        @inline (κ::$k)(d::Real) = kappa(κ,d) #TODO Add test
        @inline (κ::$k)(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}) = kappa(κ, x, y)
        @inline (κ::$k)(X::AbstractMatrix{T}, Y::AbstractMatrix{T}; obsdim::Integer=defaultobs) where {T} = kernelmatrix(κ, X, Y, obsdim=obsdim)
        @inline (κ::$k)(X::AbstractMatrix{T}; obsdim::Integer=defaultobs) where {T} = kernelmatrix(κ, X, obsdim=obsdim)
    end
end

for k in nameof.(subtypes(BaseKernel))
    @eval begin
        @deprecate($k(ρ::Real;args...),transform($k(args...),ρ))
        @deprecate($k(ρ::AbstractVector{<:Real};args...),transform($k(args...),ρ))
    end
end


text:
"""
    struct Density{M,B}
        μ::M
        base::B
    end

For measures μ and ν with μ≪ν, the density of μ with respect to ν (also called
the Radon-Nikodym derivative dμ/dν) is a function f defined on the support of ν
with the property that for any measurable a ⊂ supp(ν), μ(a) = ∫ₐ f dν.

Because this function is often difficult to express in closed form, there are
many different ways of computing it. We therefore provide a formal
representation to allow comptuational flexibilty.
"""
struct Density{M,B,L}
    μ::M
    base::B
    log::L
end

export 𝒹

export log𝒹

log𝒹(μ, base) = Density(μ, base, Val{true}())

"""
    𝒹(μ::AbstractMeasure, base::AbstractMeasure; log=false)

Compute the Radom-Nikodym derivative (or its log, if `log=false`) of μ with
respect to `base`.
"""
function 𝒹(μ::AbstractMeasure, base::AbstractMeasure; log=false)
    return Density(μ, base, Val(log))
end

(f::Density{M,B,Val{true}})(x) where {M,B} = logdensity(f.μ, f.base, x)

(f::Density{M,B,Val{false}})(x) where {M,B} = density(f.μ, f.base, x)

"""
    struct DensityMeasure{F,B} <: AbstractMeasure
        density :: F
        base    :: B
    end

A `DensityMeasure` is a measure defined by a density with respect to some other
"base" measure
"""
struct DensityMeasure{F,B,L} <: AbstractMeasure
    f    :: F
    base :: B
    log  :: L
end

function Base.show(io::IO, μ::DensityMeasure{F,B,Val{L}}) where {F,B,L}
    io = IOContext(io, :compact => true)
    print(io, "DensityMeasure ")
    print(io, "∫(", μ.f)
    print(io, ", ", μ.base)
    print(io, "; log = ", L, ")")
end

function Base.rand(rng::AbstractRNG, T::Type, d::DensityMeasure)
    x = rand(rng, T, d.base)
    WeightedMeasure(d.f(x), Dirac(x))
end

basemeasure(μ::DensityMeasure) = μ.base

logdensity(μ::DensityMeasure{F,B,Val{true}}, x) where {F,B} = μ.f(x)

density(μ::DensityMeasure{F,B,Val{false}}, x) where {F,B} = μ.f(x)

logdensity(μ::DensityMeasure{F,B,Val{false}}, x) where {F,B} = log(density(μ,x))

export ∫

"""
    ∫(f, base::AbstractMeasure)

Define a new measure in terms of a density `f` over some measure `base`.
"""
∫(f, base::AbstractMeasure) = DensityMeasure(f, base, Val(false))

∫(μ::AbstractMeasure, base::AbstractMeasure) = ∫exp(log𝒹(μ, base), base)

export ∫exp

"""
    ∫exp(f, base::AbstractMeasure; log=false)

Define a new measure in terms of a density `f` over some measure `base`.
"""
∫exp(f,μ) = DensityMeasure(f,μ,Val{true}())

# TODO: `density` and `logdensity` functions for `DensityMeasure`

function logdensity(μ::T, ν::T, x) where {T<:AbstractMeasure}
    μ==ν && return 0.0
    invoke(logdensity, Tuple{AbstractMeasure, AbstractMeasure, typeof(x)}, μ, ν, x)
end

function logdensity(μ::AbstractMeasure, ν::AbstractMeasure, x)
    α = basemeasure(μ)
    β = basemeasure(ν)

    ℓμ = @trysupport logdensity(μ, x)
    ℓν = @trysupport logdensity(ν, x)

    if α==μ && β==ν
        @warn """
        No method found for logdensity(μ, ν, x) where
        μ == $μ
        ν == $ν

        Assuming logdensity of 0.0
        """
        return ℓμ - ℓν
    end

    ℓμ - ℓν + logdensity(α, β, x)
end

logdensity(::Lebesgue, ::Lebesgue, x) = 0.0

# logdensity(::Lebesgue{ℝ}, ::Lebesgue{ℝ}, x) = zero(x)

export density

density(μ, ν::AbstractMeasure, x) = exp(logdensity(μ, ν, x))

density(μ, x) = exp(logdensity(μ, x))


text:
using LinearAlgebra, ForwardDiff

function row_col!(row,col,r,c)
    for cc in c
        for rr in r
            push!(row,convert(Int,rr))
            push!(col,convert(Int,cc))
        end
    end
    return row, col
end

v = [5.0; 5.0; 0.]
n = 1.0
μ = 0.5

P = [1. 0. 0.;
	 0. 1. 0.]

_norm(x) = sqrt(x'*x)
function vec_norm(x)
	if norm(x) == 0.0
		return ones(length(x))./norm(ones(length(x)))
	else
		x./_norm(x)
	end
end

function d_vec_norm(x)
	if norm(x) == 0.0
		y = 1.0*ones(length(x))#./norm(ones(length(x)))
		return (I - y*y'/(y'*y))/_norm(y)
	else
		(I - x*x'/(x'*x))/_norm(x)
	end
end
b = 1.0e-3*ones(2)
b = zeros(2)
d_vec_norm(b)
_norm(b)
vec_norm(b)
d_vec_norm(b)
norm(vec(d_vec_norm(b)) - vec(ForwardDiff.jacobian(vec_norm,norm(b) == 0.0 ? 1.0*ones(2) : b)))

ψ = 0.1
P*v + ψ*vec_norm(b)

##
include(joinpath(pwd(),"src/Ipopt.jl"))

nx = 4
nc = 2 + 1 + 1
function obj(x)
	b = x[1:2]
	ψ = x[3]
	s = x[4]

	return λ_global[1]*s + 0.5/mu_prev*s^2
end

function con!(c,x)
	b = x[1:2]
	ψ = x[3]
	s = x[4]

	c[1:2] = P*v + ψ*vec_norm(b)
	c[3] = μ*n - _norm(b)
	c[4] = s - min(ψ,(μ*n - _norm(b)))

	return nothing
end

function ∇con!(C,x)
	b = x[1:2]
	ψ = x[3]

	C[1:2,1:2] = ψ*d_vec_norm(b)
	C[1:2,3] = vec_norm(b)

	C[3,1:2] = -1.0*vec_norm(b)
	C[3,3] = 0.0

	if ψ < μ*n - _norm(b)
		C[4,3] = -1.0
	end

	if ψ > μ*n - _norm(b)
		C[4,1:2] = vec_norm(b)
	end

	C[4,4] = 1.0
	return nothing
end

_x = [1.0e-3; 1.0e-3; 1.0e-3; 1.0e-3]
obj(_x)
_c = zeros(nc)
con!(_c,_x)
_c
_C = zeros(nc,nx)
∇con!(_C,_x)

function eval_f(x::Vector{Float64})
	x_global .= x
  return obj(x)
end

function eval_g(x::Vector{Float64}, g::Vector{Float64})
	con!(g,x)

	return g
end

function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64})
	grad_f .= ForwardDiff.gradient(obj,x)
end

function eval_jac_g(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, values::Vector{Float64})
	if mode == :Structure
		row = []
		col = []
		row_col!(row,col,(1:nc),(1:nx))
		rows .= vec(row)
		cols .= vec(col)
	else
		∇con!(reshape(values,nc,nx),x)
	end
end

function eval_h(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, obj_factor::Float64, lambda::Vector{Float64}, values::Vector{Float64})
	if mode == :Structure
		nothing
	else
		nothing
	end
end

# This tests callbacks.
function intermediate(alg_mod::Int, iter_count::Int,
  obj_value::Float64, inf_pr::Float64, inf_du::Float64, mu::Float64,
  d_norm::Float64, regularization_size::Float64, alpha_du::Float64, alpha_pr::Float64,
  ls_trials::Int)

	println("CALLBACK")
  println("   mu_prev: $(mu_prev)")
  println("   mu: $(mu)")
  println("   x: $(x_global)")

  if mu != mu_prev
    println("-mu difference-")
    println("-AL update-")
    λ_global[1] += 1/mu_prev*x_global[4]
    global mu_prev = mu
  end

  return true#!(norm(x_global[9:11],Inf) < 1.0e-3 && norm(eval_g(x_global,zeros(7))) < 1.0e-3) #iter_count < 1  # Interrupts after one iteration.
end

x_L = [-1.0e8; -1.0e8; 0.0; 0.0]
x_U = 1.0e8*ones(nx)

g_L = zeros(nc)
g_U = [zeros(2);1.0e8;0.0]

prob = Ipopt.createProblem(nx, x_L, x_U, nc, g_L, g_U, nc*nx, nx*nx,
                     eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h)
Ipopt.addOption(prob,"hessian_approximation","limited-memory")
# Ipopt.addOption(prob,"mu_strategy","monotone")

Ipopt.addOption(prob,"constr_viol_tol",1.0e-3)
Ipopt.addOption(prob,"tol",1.0e-3)
Ipopt.setIntermediateCallback(prob, intermediate)

global x_global = zeros(nx)
global λ_global = zeros(1)
global mu_prev = 0.1

prob.x .= 1.0e-3*rand(nx)
x_global .= prob.x
@time solvestat = Ipopt.solveProblem(prob)

b_sol = prob.x[1:2]
ψ_sol = prob.x[3]
s_sol = prob.x[4]
λ_global
mu_prev

norm(v[1:2])
vec_norm(v[1:2])'*vec_norm(b_sol)

con!(_c,prob.x)
_c
_C = zeros(nc,nx)
∇con!(_C,prob.x)
_C


text:
using Documenter, JuliaDBMeta

makedocs(
    format = :html,
    sitename = "JuliaDBMeta",
    authors = "Pietro Vertechi",
    pages = [
        "Introduction" => "index.md",
        "Getting Started" => "getting_started.md",
        "Row-wise macros" => "row_macros.md",
        "Column-wise macros" => "column_macros.md",
        "Selection" => "selection.md",
        "Pipeline macros" => "pipeline_macros.md",
        "Grouping operations" => "grouping.md",
        "Out-of-core support" => "out_of_core.md",
        "Plotting" => "plotting.md",
        "Tutorial" => "tutorial.md",
    ]
)

deploydocs(
    repo = "github.com/piever/JuliaDBMeta.jl.git",
    target = "build",
)


text:
"""
    allele_freq(locus::T) where T<:GenoArray
Return a `Dict` of allele frequencies of a single locus in a `PopData`
object.
"""
@inline function allele_freq(locus::GenoArray)
    d = Dict{Int16,Float32}()
    flat_alleles = alleles(locus)
    len = length(flat_alleles)
    len == 0 && return d
    @inbounds @simd for allele in flat_alleles
        @inbounds d[allele] = get!(d, allele, 0) + 1/len
    end
    return d
end

"""
    allele_freq_vec(locus::T) where T<:GenoArray
Return a Vector of allele frequencies of a single locus in a `PopData`
object. Similar to `allele_freq()`, except it returns only the frequencies,
without the allele names, meaning they can be in any order. This is useful
for getting the expected genotype frequencies.
"""
@inline function allele_freq_vec(locus::GenoArray)
    flat_alleles = alleles(locus)
    len = length(flat_alleles)
    d = [count(i -> i == j, flat_alleles) for j in unique(flat_alleles)]
    return d ./ len
end

@inline function allele_freq_vec(::Missing)
    return missing
end

#TODO add to docs (API/AlleleFreq)
"""
    avg_allele_freq(allele_dicts::AbstractVector{T}) where T<:Dict{Int16,Float32}
Takes a Vector of Dicts generated by `allele_freq` and returns a Dict of the average
allele frequencies. This is typically done to calculate average allele frequencies
across populations.

***Example***
```
cats = nancycats();
alleles_df = DataFrames.combine(
    groupby(cats.loci, [:locus, :population]),
    :genotype => allele_freq => :alleles
);
DataFrames.combine(
    groupby(alleles_df, :locus),
    :alleles => (i -> sum(avg_allele_freq(i).^2)) => :avg_freq
)
```
"""
function avg_allele_freq(allele_dicts::AbstractVector{T}) where T<:Dict{Int16,Float32}
   sum_dict = Dict{Int16, Tuple{Float32, Int}}()
   # remove any dicts with no entries (i.e. from a group without that locus)
   allele_dicts = allele_dicts[length.(allele_dicts) .> 0]
   # create a list of all the alleles
   all_alleles = keys.(allele_dicts) |> Base.Iterators.flatten |> collect |> unique
   # populate the sum dict with allele frequency and n for each allele
   @inbounds for allele in all_alleles
       for allele_dict in allele_dicts
           sum_dict[allele] = get!(sum_dict, allele, (0., 0)) .+ (get!(allele_dict, allele, 0.), 1)
       end
   end
   avg_dict = Dict{Int16, Float32}()
   # collapse the sum dict into a dict of averages
   @inbounds for (key, value) in sum_dict
       freq_sum, n = value
       avg = freq_sum / n
       # drop zeroes
       if !iszero(avg)
           @inbounds avg_dict[key] = avg
       end
   end
   return avg_dict
end

"""
    allele_freq(data::PopData, locus::String; population::Bool = false)
Return a `Dict` of allele frequencies of a single locus in a `PopData`
object. Use `population = true` to return a table of allele frequencies
of that locus per population.
### Example
```
cats = nancycats()
allele_freq(cats, "fca8")
allele_freq(cats, "fca8", population = true)
```
"""
function allele_freq(data::PopData, locus::String, population::Bool=false)
    if !population
        data.loci[data.loci.locus .== locus, :genotype] |> allele_freq
    else
        tmp = groupby(data.loci[data.loci.locus .== locus, :], :population)
        DataFrames.combine(tmp, :genotype => allele_freq => :frequency)
    end
end

#TODO add to docs (API)
"""
    allele_freq(allele::Int, genos::T) where T<:GenoArray
Return the frequency of an `allele` from a vector of `genotypes`

### Example
```
using DataFramesMeta
ncats = nancycats();
ncats_sub @where(ncats.loci, :locus .== "fca8", :genotype .!== missing)
pop_grp = groupby(ncats_sub, :population)
DataFrames.combine(pop_grp, :genotype => (geno,) -> allele_freq(146, geno))
```
"""
function allele_freq(allele::Int, genos::T) where T<:GenoArray
    ploidy = (length.(genos) |> unique)[1]
    count(allele .== Base.Iterators.flatten(genos))/(ploidy*nonmissing(genos))
end

"""
    geno_count_observed(locus::T) where T<:GenoArray
Return a `Dict` of genotype counts of a single locus in a
`PopData` object.
"""
@inline function geno_count_observed(locus::T) where T<:GenoArray
    # conditional testing if all genos are missing
    all(ismissing.(locus)) && return missing
    d = Dict{Tuple, Float32}()
    @inbounds for genotype in skipmissing(locus)
        # sum up non-missing genotypes
        d[genotype] = get!(d, genotype, 0) + 1
    end
    return d
end

"""
    geno_count_expected(locus::T) where T<:GenoArray
Return a `Dict` of the expected genotype counts of a single locus in a
`PopData` object. Expected counts are calculated as the product of observed
allele frequencies multiplied by the number of non-missing genotypes.
"""
function geno_count_expected(locus::T) where T<:GenoArray
    #count number of non-missing genotypes in the locus
    n = nonmissing(locus)

    # Get expected number of genotypes in a locus
    ## get the observed allele frequencies
    allele_dict = allele_freq(locus)

    ## split the appropriate pairs into their own vectors
    alle, freq = collect(keys(allele_dict)), collect(values(allele_dict))
    #return alle, freq
    ## calculate expected genotype frequencies by multiplying all-by-all x n
    expected_genotype_freq = vec(freq * freq' .* n)

    # reform genotype frequencies with same all-by-all approach
    genos = reverse.(Base.Iterators.product(alle, alle) |> collect |> vec)
    expected = Dict{Tuple, Float64}()
    for (geno, freq) in zip(genos, expected_genotype_freq)
        expected[geno] = get!(expected, geno, 0.0) + freq
    end

    return expected
end

"""
    geno_freq(locus::T) where T<:GenoArray
Return a `Dict` of genotype frequencies of a single locus in a
`PopData` object.
"""
@inline function geno_freq(locus::T) where T<:GenoArray
    # conditional testing if all genos are missing
    all(ismissing.(locus)) && return missing
    d = Dict{Tuple, Float32}()
    @inbounds for genotype in skipmissing(locus)
        # sum up non-missing genotypes
        d[genotype] = get!(d, genotype, 0) + 1
    end
    total = values(d) |> sum    # sum of all non-missing genotypes
    [d[i] = d[i] / total for i in keys(d)] # genotype count/total
    return d
end

"""
    geno_freq(data::PopData, locus::String; population::Bool = false)
Return a `Dict` of genotype frequencies of a single locus in a `PopData`
object. Use `population = true` to return a table of genotype frequencies
of that locus per population.
### Example
```
cats = nancycats()
geno_freq(cats, "fca8")
geno_freq(cats, "fca8", population = true)
```
"""
function geno_freq(data::PopData, locus::String, population::Bool=false)
    if !population
        tmp = data.loci[data.loci.locus .== locus, :genotype]
        geno_freq(tmp)
    else
        tmp = data.loci[data.loci.locus .== locus, :]
        DataFrames.combine(groupby(tmp, :population), :genotype => geno_freq => :freq)
    end
end

#TODO replace allele creation with Base.Iterators.product
"""
    geno_freq_expected(locus::T) where T<:GenoArray
Return a `Dict` of the expected genotype frequencies of a single locus in a
`PopData` object. Expected frequencies are calculated as the product of
observed allele frequencies.
"""
function geno_freq_expected(locus::T) where T<:GenoArray
    # Get expected number of genotypes in a locus
    ## get the observed allele frequencies
    allele_dict = allele_freq(locus)

    ## split the appropriate pairs into their own vectors
    alle, freq = collect(keys(allele_dict)), collect(values(allele_dict))

    ## calculate expected genotype frequencies by multiplying all-by-all
    expected_genotype_freq = vec(freq * freq')

    # reform genotype frequencies with same all-by-all approach
    genos = reverse.(Base.Iterators.product(alle, alle) |> collect |> vec)

    # reform genotype frequencies into a Dict
    expected = Dict{Tuple, Float64}()
    for (geno, freq) in zip(genos, expected_genotype_freq)
        expected[geno] = get!(expected, geno, 0) + freq
    end

    return expected
end

"""
    geno_freq_expected(data::PopData, locus::String; population::Bool = false)
Return a `Dict` of expected genotype frequencies of a single locus in a
`PopData` object. Use `population = true` to return a table of expected genotype
frequencies of that locus per population.
### Example
```
cats = nancycats()
geno_freq_expeced(cats, "fca8")
geno_freq_expected(cats, "fca8", population = true)
```
"""
function geno_freq_expected(data::PopData, locus::String, population::Bool=false)
    if !population
        tmp = data.loci[data.loci.locus .== locus, :genotype]
        geno_freq_expected(tmp)
    else
        tmp = data.loci[data.loci.locus .== locus, :]
        DataFrames.combine(groupby(tmp, :population), :genotype => geno_freq_expected => :freq)
    end
end


text:
"""
A `GenericCompositeOperator` contains a list of operators that are applied
consecutively to any input. It is a generic operator.
"""
struct GenericCompositeOperator <: AbstractOperator
    operators
end

# This constructor takes several operators as arguments and ensures that
# the spaces are compatible. This step can by sidestepped by calling the
# default inner constructor directly with a vector or tuple of operators instead.
function GenericCompositeOperator(operators::AbstractOperator...)
    for i in 1:length(operators)-1
        @assert coefficienttype(dest_space(operators[i])) == coefficienttype(src_space(operators[i+1]))
    end
    # Pass the tuple of operators to the inner constructor
    GenericCompositeOperator(operators)
end

src_space(op::GenericCompositeOperator) = src_space(op.operators[1])
dest_space(op::GenericCompositeOperator) = dest_space(op.operators[end])

# Generic functions for composite types:
components(op::GenericCompositeOperator) = op.operators
component(op::GenericCompositeOperator, j::Int) = op.operators[j]

# If the GenericCompositeOperator happens to be a composite operator of
# product operators, then the result is a product operator too. You can select
# its components with the factor routine.
function factor(op::GenericCompositeOperator, j::Int)
    GenericCompositeOperator(map(t -> factor(t, j), components(op))...)
end
factors(op::GenericCompositeOperator) = tuple([factor(op, j) for j in 1:nfactors(op)]...)
nfactors(op::GenericCompositeOperator) = nfactors(component(op,1))

function apply(comp::GenericCompositeOperator, fun; options...)
    output = fun
    for op in components(comp)
        input = output
        # TODO: clean this op. The apply of a generic operator can accept options,
        # but the apply of a DictionaryOperator can not.
        if op isa DictionaryOperator
            output = apply(op, input)
        else
            output = apply(op, input; options...)
        end
    end
    output
end

"Can an operator allocate storage for its expected output?"
can_allocate_output(op::DictionaryOperator) = true
# The answer is no in general, but yes if the output is the span of a dictionary
# - Yes for all dictionary operators

# - AbstractOperator: it depends on the output, which we determine using dispatch
can_allocate_output(op::AbstractOperator) = _can_allocate_output(op, dest_space(op))
_can_allocate_output(op, span::Span) = true
_can_allocate_output(op, space::FunctionSpace) = false

allocate_output(op::DictionaryOperator{T}) where {T} = zeros(T, dest(op))

Display.combinationsymbol(op::GenericCompositeOperator) = Display.Symbol('∘')
Display.displaystencil(op::GenericCompositeOperator) =
    composite_displaystencil(op; reversecomponents=true)
show(io::IO, mime::MIME"text/plain", op::GenericCompositeOperator) =
    composite_show(io, mime, op)

"""
A `CompositeOperator` consists of a sequence of operators that are applied
consecutively.

Whenever possible, scratch space is allocated to hold intermediate results.
"""
struct CompositeOperator{T} <: DictionaryOperator{T}
    # We explicitly store src and dest, because that information may be lost
    # when the list of operators is optimized (for example, an Identity mapping
    # between two spaces could disappear).
    src     ::  Dictionary
    dest    ::  Dictionary
    "The list of operators"
    operators
    "Scratch space for the result of each operator, except the last one"
    # We don't need it for the last one, because the final result goes to coef_dest
    scratch
end

# Generic functions for composite types:
components(op::CompositeOperator) = op.operators
component(op::CompositeOperator, j::Int) = op.operators[j]

function factor(op::CompositeOperator, j::Int)
    compose(map(t -> factor(t, j), components(op))...)
end
factors(op::CompositeOperator) = tuple([factor(op, j) for j in 1:nfactors(op)]...)
nfactors(op::CompositeOperator) = nfactors(component(op,1))

isinplace(op::CompositeOperator) = all(map(isinplace, op.operators))
isdiag(op::CompositeOperator) = all(map(isdiag, op.operators))

function compose_and_simplify(composite_src::Dictionary, composite_dest::Dictionary, operators::DictionaryOperator...; simplify = true)
    # Check operator compatibility
    for i in 1:length(operators)-1
        @assert size(dest(operators[i])) == size(src(operators[i+1]))
#       TODO: at one point we should enable strict checking again as follows:
#        @assert dest(operators[i]) == src(operators[i+1])
    end

    # Checks pass, now apply some simplifications
    if simplify
        # Flatten away nested CompositeOperators
        operators = flatten(CompositeOperator, operators...)
        operators = compose_and_simplify(operators...)
    end

    L = length(operators)
    if L == 1
        return wrap_operator(composite_src, composite_dest, operators[1])
    end
    T = promote_type(map(eltype, operators)...)
    # We are going to reserve scratch space, but only for operators that are not
    # in-place. We do reserve scratch space for the first operator, even if it
    # is in-place, because we may want to call the composite operator out of place.
    # In that case we need a place to store the result of the first operator.
    scratch_array = Any[zeros(T, dest(operators[1]))]
    for m = 2:L-1
        if ~isinplace(operators[m])
            push!(scratch_array, zeros(T, dest(operators[m])))
        end
    end
    scratch = tuple(scratch_array...)
    CompositeOperator{T}(composite_src, composite_dest, operators, scratch)
end

unsafe_wrap_operator(src, dest, op::CompositeOperator{T}) where T =
    CompositeOperator{T}(src, dest, op.operators, op.scratch)

apply_inplace!(op::CompositeOperator, coef_srcdest) =
    apply_inplace_composite!(op, coef_srcdest, op.operators)

apply!(op::CompositeOperator, coef_dest, coef_src) =
    apply_composite!(op, coef_dest, coef_src, op.operators, op.scratch)

function apply_composite!(op::CompositeOperator, coef_dest, coef_src, operators, scratch)
    L = length(operators)
    apply!(operators[1], scratch[1], coef_src)
    l = 1
    for i in 2:L-1
        if isinplace(operators[i])
            apply_inplace!(operators[i], scratch[l])
        else
            apply!(operators[i], scratch[l+1], scratch[l])
            l += 1
        end
    end
    # We loose a little bit of efficiency if the last operator was in-place
    apply!(operators[L], coef_dest, scratch[l])
    # Return coef_dest even though the last apply! should also do that, to
    # help type-inference continue afterwards
    coef_dest
end

# Below is the ideal scenario for lengths 2 and 3, written explicitly.
# function apply_composite!(op::CompositeOperator, coef_dest, coef_src, operators::NTuple{2}, scratch)
#     if isinplace(operators[2])
#         apply!(operators[1], coef_dest, coef_src)
#         apply!(operators[2], coef_dest)
#     else
#         apply!(operators[1], scratch[1], coef_src)
#         apply!(operators[2], coef_dest, scratch[1])
#     end
#     coef_dest
# end
#
# function apply_composite!(op::CompositeOperator, coef_dest, coef_src, operators::NTuple{3}, scratch)
#     ip2 = isinplace(operators[2])
#     ip3 = isinplace(operators[3])
#     if ip2
#         if ip3
#             # 2 and 3 are in-place
#             apply!(operators[1], coef_dest, coef_src)
#             apply!(operators[2], coef_dest)
#             apply!(operators[3], coef_dest)
#         else
#             # 2 is in place, 3 is not
#             apply!(operators[1], scratch[1], coef_src)
#             apply!(operators[2], scratch[1])
#             apply!(operators[3], coef_dest, scratch[1])
#         end
#     else
#         if ip3
#             # 2 is not in place, but 3 is
#             apply!(operators[1], scratch[1], coef_src)
#             apply!(operators[2], coef_dest, scratch[1])
#             apply!(operators[3], coef_dest)
#         else
#             # noone is in place
#             apply!(operators[1], scratch[1], coef_src)
#             apply!(operators[2], scratch[2], scratch[1])
#             apply!(operators[3], coef_dest, scratch[2])
#         end
#     end
#     coef_dest
# end

function apply_inplace_composite!(op::CompositeOperator, coef_srcdest, operators)
    for operator in operators
        apply!(operator, coef_srcdest)
    end
    coef_srcdest
end

inv(op::CompositeOperator) = unsafe_wrap_operator(dest(op), src(op), (*)(map(inv, op.operators)...))
adjoint(op::CompositeOperator) = unsafe_wrap_operator(dest(op), src(op), (*)(map(adjoint, op.operators)...))

conj(op::CompositeOperator{T}) where T  = CompositeOperator{T}(src(op), dest(op), map(conj,op.operators), op.scratch)
conj(op::CompositeOperator{T}) where {T<:Real} = op

sqrt(op::CompositeOperator{T}) where T  = CompositeOperator{T}(src(op), dest(op), map(sqrt,op.operators), op.scratch)

apply(op1::AbstractOperator, op2::AbstractOperator) = compose(op2,op1)
apply(op1::DictionaryOperator, op2::AbstractOperator) = compose(op2,op1)

# Don't do anything if we have just one operator
compose(op::AbstractOperator) = op

compose(ops::DictionaryOperator...) = compose_and_simplify(src(ops[1]), dest(ops[end]), ops...)
compose(ops::AbstractOperator...) = GenericCompositeOperator(flatten(GenericCompositeOperator, ops...)...)

sparse(op::CompositeOperator; options...) = *([sparse(opi; options...) for opi in components(op)[end:-1:1]]...)

CompositeOperators = Union{CompositeOperator,GenericCompositeOperator}

function stencilarray(op::CompositeOperators)
    A = Any[]
    push!(A,component(op,ncomponents(op)))
    for i in ncomponents(op)-1:-1:1
        push!(A, " * ")
        push!(A, component(op,i))
    end
    A
end

stencil_parentheses(op::CompositeOperators) = true
object_parentheses(op::CompositeOperators) = true


text:
#####
##### pyMem.jl
##### Script to save the memory in a python
##### usable format for experimentation
#####

using AlphaZero
using Serialization: deserialize
using NPZ

include("games.jl")
GAME = get(ENV, "GAME", "double-dummy")
SelectedGame = GAME_MODULE[GAME]
using .SelectedGame: Game

memory_path = "sessions/double-dummy/mem.data"
memory = deserialize(memory_path)

values = [sample.z for sample in memory]
states = zeros(UInt8, size(memory[1].s.board)[1], size(memory[1].s.board)[2], size(memory[1].s.board)[3], length(memory))
policies = zeros(Float64, 52, length(memory))

for i in 1:length(memory)
	for j in 1:size(memory[1].s.board)[3]
		for k in 1:size(memory[1].s.board)[2]
			for l in 1:size(memory[1].s.board)[1]
				states[l,k,j,i] = memory[i].s.board[l,k,j]
			end
		end
	end

	amask = GI.actions_mask(Game(memory[i].s))
	actions = zeros(size(amask))
	actions[amask] = memory[i].π
	policies[:,i] = actions
end

npzwrite("pyMem/states.npz", states)
npzwrite("pyMem/values.npz", values)
npzwrite("pyMem/states.policies", policies)

max_stars_repo_name: btheodorou/Bridge-NN-Research

max_issues_repo_name: btheodorou/Bridge-NN-Research

max_forks_repo_name: btheodorou/Bridge-NN-Research


text:
module TestQuery

  using GBIF
  using DataFrames
  using Query
  using Test

  t = taxon("Mammalia", strict=false)
  set = occurrences(t)
  [occurrences!(set) for i in 1:10]

  tdf = view(set) |>
    @filter(_.rank == "SPECIES") |>
    @map({_.key, _.taxon.name, _.country}) |>
    DataFrame

  @test typeof(tdf) <: DataFrame

end


text:
# list of integrators (integrator, tableau, sim_id)

function get_tableau_list_erk(T=Float64)
    tableau_list = (
        # (IntegratorERK,               getTableauExplicitEuler(),     "euler"  ),
        # (IntegratorERK,               getTableauExplicitMidpoint(),  "emp"    ),
        # (IntegratorERK,               getTableauHeun(),              "heun"   ),
        # (IntegratorERK,               getTableauKutta(),             "kutta"  ),
        (IntegratorERK,               getTableauERK4(),              "erk4_16"),
        (IntegratorERK,               getTableauERK438(),            "erk4_38"),
    )
end

function get_tableau_list_firk(T=Float64)
    tableau_list = (
        (IntegratorFIRK,              getTableauGLRK(1),      "glrk1"      ),
        (IntegratorFIRK,              getTableauGLRK(2),      "glrk2"      ),
        (IntegratorFIRK,              getTableauGLRK(3),      "glrk3"      ),
        (IntegratorFIRK,              getTableauGLRK(4),      "glrk4"      ),
        (IntegratorFIRK,              getTableauGLRK(5),      "glrk5"      ),
        (IntegratorFIRK,              getTableauGLRK(6),      "glrk6"      ),
        (IntegratorFIRK,              getTableauLobIIIA2(),   "lobIIIA2"   ),
        (IntegratorFIRK,              getTableauLobIIIA3(),   "lobIIIA3"   ),
        (IntegratorFIRK,              getTableauLobIIIA4(),   "lobIIIA4"   ),
        (IntegratorFIRK,              getTableauLobIIIB2(),   "lobIIIB2"   ),
        (IntegratorFIRK,              getTableauLobIIIB3(),   "lobIIIB3"   ),
        (IntegratorFIRK,              getTableauLobIIIB4(),   "lobIIIB4"   ),
        (IntegratorFIRK,              getTableauLobIIIC2(),   "lobIIIC2"   ),
        (IntegratorFIRK,              getTableauLobIIIC3(),   "lobIIIC3"   ),
        (IntegratorFIRK,              getTableauLobIIIC4(),   "lobIIIC4"   ),
        (IntegratorFIRK,              getTableauLobIIID2(),   "lobIIID2"   ),
        (IntegratorFIRK,              getTableauLobIIID3(),   "lobIIID3"   ),
        (IntegratorFIRK,              getTableauLobIIID4(),   "lobIIID4"   ),
        (IntegratorFIRK,              getTableauLobIIIE2(),   "lobIIIE2"   ),
        (IntegratorFIRK,              getTableauLobIIIE3(),   "lobIIIE3"   ),
        (IntegratorFIRK,              getTableauLobIIIE4(),   "lobIIIE4"   ),
        (IntegratorFIRK,              getTableauLobIIIF2(),   "lobIIIF2"   ),
        (IntegratorFIRK,              getTableauLobIIIF3(),   "lobIIIF3"   ),
        (IntegratorFIRK,              getTableauLobIIIF4(),   "lobIIIF4"   ),
        (IntegratorFIRK,              getTableauLobIIIG2(),   "lobIIIG2"   ),
        (IntegratorFIRK,              getTableauLobIIIG3(),   "lobIIIG3"   ),
        (IntegratorFIRK,              getTableauLobIIIG4(),   "lobIIIG4"   ),
    )
end

function get_tableau_list_vprk(T=Float64)
    tableau_list = (
        (IntegratorVPRKpNone,         getTableauVPGLRK(1, T=T), "vprk_glrk1"   ),
        (IntegratorVPRKpNone,         getTableauVPGLRK(2, T=T), "vprk_glrk2"   ),
        (IntegratorVPRKpNone,         getTableauVPGLRK(3, T=T), "vprk_glrk3"   ),
        (IntegratorVPRKpNone,         getTableauVPGLRK(4, T=T), "vprk_glrk4"   ),
        (IntegratorVPRKpNone,         getTableauVPGLRK(5, T=T), "vprk_glrk5"   ),
        (IntegratorVPRKpNone,         getTableauVPGLRK(6, T=T), "vprk_glrk6"   ),
        (IntegratorVPRKpNone,         getTableauVPSRK3(),       "vprk_srk3"    ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIA2(),   "vprk_lobIIIA2"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIA3(),   "vprk_lobIIIA3"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIA4(),   "vprk_lobIIIA4"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIB2(),   "vprk_lobIIIB2"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIB3(),   "vprk_lobIIIB3"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIB4(),   "vprk_lobIIIB4"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIC2(),   "vprk_lobIIIC2"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIC3(),   "vprk_lobIIIC3"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIC4(),   "vprk_lobIIIC4"),
        (IntegratorVPRKpNone,         getTableauVPLobIIID2(),   "vprk_lobIIID2"),
        (IntegratorVPRKpNone,         getTableauVPLobIIID3(),   "vprk_lobIIID3"),
        (IntegratorVPRKpNone,         getTableauVPLobIIID4(),   "vprk_lobIIID4"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIE2(),   "vprk_lobIIIE2"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIE3(),   "vprk_lobIIIE3"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIE4(),   "vprk_lobIIIE4"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIE2(),   "vprk_lobIIIF2"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIE3(),   "vprk_lobIIIF3"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIE4(),   "vprk_lobIIIF4"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIE2(),   "vprk_lobIIIG2"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIE3(),   "vprk_lobIIIG3"),
        (IntegratorVPRKpNone,         getTableauVPLobIIIE4(),   "vprk_lobIIIG4"),
    )
end

function get_tableau_list_vprk_projection(T=Float64)
    tableau_list = (
        (IntegratorVPRKpNone,         getTableauVPGLRK(1, T=T), "vprk_glrk1_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPGLRK(1, T=T), "vprk_glrk1_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPGLRK(1, T=T), "vprk_glrk1_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPGLRK(1, T=T), "vprk_glrk1_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPGLRK(1, T=T), "vprk_glrk1_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPGLRK(2, T=T), "vprk_glrk2_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPGLRK(2, T=T), "vprk_glrk2_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPGLRK(2, T=T), "vprk_glrk2_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPGLRK(2, T=T), "vprk_glrk2_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPGLRK(2, T=T), "vprk_glrk2_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPGLRK(3, T=T), "vprk_glrk3_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPGLRK(3, T=T), "vprk_glrk3_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPGLRK(3, T=T), "vprk_glrk3_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPGLRK(3, T=T), "vprk_glrk3_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPGLRK(3, T=T), "vprk_glrk3_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPGLRK(4, T=T), "vprk_glrk4_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPGLRK(4, T=T), "vprk_glrk4_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPGLRK(4, T=T), "vprk_glrk4_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPGLRK(4, T=T), "vprk_glrk4_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPGLRK(4, T=T), "vprk_glrk4_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPGLRK(5, T=T), "vprk_glrk5_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPGLRK(5, T=T), "vprk_glrk5_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPGLRK(5, T=T), "vprk_glrk5_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPGLRK(5, T=T), "vprk_glrk5_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPGLRK(5, T=T), "vprk_glrk5_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPGLRK(6, T=T), "vprk_glrk6_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPGLRK(6, T=T), "vprk_glrk6_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPGLRK(6, T=T), "vprk_glrk6_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPGLRK(6, T=T), "vprk_glrk6_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPGLRK(6, T=T), "vprk_glrk6_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPSRK3(),       "vprk_srk3_pnone"       ),
        (IntegratorVPRKpStandard,     getTableauVPSRK3(),       "vprk_srk3_pstandard"   ),
        (IntegratorVPRKpSymmetric,    getTableauVPSRK3(),       "vprk_srk3_psymmetric"  ),
        (IntegratorVPRKpSymplectic,   getTableauVPSRK3(),       "vprk_srk3_psymplectic" ),
        (IntegratorVPRKpMidpoint,     getTableauVPSRK3(),       "vprk_srk3_pmidpoint"   ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIA2(),   "vprk_lobIIIA2_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIA2(),   "vprk_lobIIIA2_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIA2(),   "vprk_lobIIIA2_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIA2(),   "vprk_lobIIIA2_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIA2(),   "vprk_lobIIIA2_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIA3(),   "vprk_lobIIIA3_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIA3(),   "vprk_lobIIIA3_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIA3(),   "vprk_lobIIIA3_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIA3(),   "vprk_lobIIIA3_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIA3(),   "vprk_lobIIIA3_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIA4(),   "vprk_lobIIIA4_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIA4(),   "vprk_lobIIIA4_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIA4(),   "vprk_lobIIIA4_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIA4(),   "vprk_lobIIIA4_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIA4(),   "vprk_lobIIIA4_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIB2(),   "vprk_lobIIIB2_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIB2(),   "vprk_lobIIIB2_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIB2(),   "vprk_lobIIIB2_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIB2(),   "vprk_lobIIIB2_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIB2(),   "vprk_lobIIIB2_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIB3(),   "vprk_lobIIIB3_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIB3(),   "vprk_lobIIIB3_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIB3(),   "vprk_lobIIIB3_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIB3(),   "vprk_lobIIIB3_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIB3(),   "vprk_lobIIIB3_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIB4(),   "vprk_lobIIIB4_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIB4(),   "vprk_lobIIIB4_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIB4(),   "vprk_lobIIIB4_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIB4(),   "vprk_lobIIIB4_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIB4(),   "vprk_lobIIIB4_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIC2(),   "vprk_lobIIIC2_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIC2(),   "vprk_lobIIIC2_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIC2(),   "vprk_lobIIIC2_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIC2(),   "vprk_lobIIIC2_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIC2(),   "vprk_lobIIIC2_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIC3(),   "vprk_lobIIIC3_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIC3(),   "vprk_lobIIIC3_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIC3(),   "vprk_lobIIIC3_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIC3(),   "vprk_lobIIIC3_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIC3(),   "vprk_lobIIIC3_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIC4(),   "vprk_lobIIIC4_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIC4(),   "vprk_lobIIIC4_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIC4(),   "vprk_lobIIIC4_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIC4(),   "vprk_lobIIIC4_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIC4(),   "vprk_lobIIIC4_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIID2(),   "vprk_lobIIID2_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIID2(),   "vprk_lobIIID2_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIID2(),   "vprk_lobIIID2_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIID2(),   "vprk_lobIIID2_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIID2(),   "vprk_lobIIID2_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIID3(),   "vprk_lobIIID3_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIID3(),   "vprk_lobIIID3_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIID3(),   "vprk_lobIIID3_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIID3(),   "vprk_lobIIID3_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIID3(),   "vprk_lobIIID3_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIID4(),   "vprk_lobIIID4_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIID4(),   "vprk_lobIIID4_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIID4(),   "vprk_lobIIID4_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIID4(),   "vprk_lobIIID4_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIID4(),   "vprk_lobIIID4_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIE2(),   "vprk_lobIIIE2_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIE2(),   "vprk_lobIIIE2_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIE2(),   "vprk_lobIIIE2_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIE2(),   "vprk_lobIIIE2_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIE2(),   "vprk_lobIIIE2_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIE3(),   "vprk_lobIIIE3_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIE3(),   "vprk_lobIIIE3_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIE3(),   "vprk_lobIIIE3_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIE3(),   "vprk_lobIIIE3_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIE3(),   "vprk_lobIIIE3_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIE4(),   "vprk_lobIIIE4_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIE4(),   "vprk_lobIIIE4_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIE4(),   "vprk_lobIIIE4_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIE4(),   "vprk_lobIIIE4_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIE4(),   "vprk_lobIIIE4_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIF2(),   "vprk_lobIIIF2_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIF2(),   "vprk_lobIIIF2_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIF2(),   "vprk_lobIIIF2_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIF2(),   "vprk_lobIIIF2_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIF2(),   "vprk_lobIIIF2_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIF3(),   "vprk_lobIIIF3_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIF3(),   "vprk_lobIIIF3_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIF3(),   "vprk_lobIIIF3_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIF3(),   "vprk_lobIIIF3_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIF3(),   "vprk_lobIIIF3_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIF4(),   "vprk_lobIIIF4_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIF4(),   "vprk_lobIIIF4_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIF4(),   "vprk_lobIIIF4_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIF4(),   "vprk_lobIIIF4_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIF4(),   "vprk_lobIIIF4_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIG2(),   "vprk_lobIIIG2_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIG2(),   "vprk_lobIIIG2_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIG2(),   "vprk_lobIIIG2_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIG2(),   "vprk_lobIIIG2_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIG2(),   "vprk_lobIIIG2_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIG3(),   "vprk_lobIIIG3_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIG3(),   "vprk_lobIIIG3_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIG3(),   "vprk_lobIIIG3_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIG3(),   "vprk_lobIIIG3_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIG3(),   "vprk_lobIIIG3_pmidpoint"  ),
        (IntegratorVPRKpNone,         getTableauVPLobIIIG4(),   "vprk_lobIIIG4_pnone"      ),
        (IntegratorVPRKpStandard,     getTableauVPLobIIIG4(),   "vprk_lobIIIG4_pstandard"  ),
        (IntegratorVPRKpSymmetric,    getTableauVPLobIIIG4(),   "vprk_lobIIIG4_psymmetric" ),
        (IntegratorVPRKpSymplectic,   getTableauVPLobIIIG4(),   "vprk_lobIIIG4_psymplectic"),
        (IntegratorVPRKpMidpoint,     getTableauVPLobIIIG4(),   "vprk_lobIIIG4_pmidpoint"  ),
        )
    end

function get_tableau_list_vprk_radau(T=Float64)
    tableau_list = (
        (IntegratorVPRK,              getTableauVPRadIIAIIA2(), "vprk_radIIA2"),
        (IntegratorVPRK,              getTableauVPRadIIAIIA3(), "vprk_radIIA3"),
    )
end

function get_tableau_list_splitting(T=Float64)
    tableau_list = (
        (IntegratorSplitting,  getTableauLieA(),             "LieA"           ),
        (IntegratorSplitting,  getTableauLieB(),             "LieB"           ),
        (IntegratorSplitting,  getTableauStrang(),           "Strang"         ),
        (IntegratorSplitting,  getTableauMcLachlan2(),       "McLachlan2"     ),
        (IntegratorSplitting,  getTableauMcLachlan4(),       "McLachlan4"     ),
        (IntegratorSplitting,  getTableauTripleJump(),       "TripleJump"     ),
        (IntegratorSplitting,  getTableauSuzukiFractal(),    "SuzukiFractal"  ),
    )
end