text:
#=
    Sorting networks that process 17,..25 values are given.

    These sorting networks should minimize conditional swaps.
    There may exist one or more comparable alternative solutions;
      those may perform marginally better for some kinds of hardware.

    Functions are written with groups of internally paralleli1zable
    statements given together and sequential steps are separated.

    `swapsort(x1::T, .., xn::T)::NTuple{n,T}`
    `swapsort(tup::NTuple{n,T})::NTuple{n,T}`

    networks were selected using software by John Gamble
        http://pages.ripco.net/~jgamble/nw.html
        http://search.cpan.org/dist/Algorithm-Networksort/
=#

const MORE_ARGS_MIN = 17
const MORE_ARGS_MAX = 25

#=
   sort NTuples of 17..25 values
=#

for N in collect(MORE_ARGS_MIN:MORE_ARGS_MAX)
    @eval swapsort(x::NTuple{$N, T}) where T = swapsort(x...)
end

#=
    sorting networks for 17..25 values
    using known compact realizations
=#

#    sort 17 values with 74 minmaxs in 12 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
                f::T, g::T, h::T, i::T, j::T,
                k::T, l::T, m::T, n::T, o::T,
                p::T, q::T) where {T}

    a, q = min_max(a, q)
    b, j = min_max(b, j)
    c, k = min_max(c, k)
    d, l = min_max(d, l)
    e, m = min_max(e, m)
    f, n = min_max(f, n)
    g, o = min_max(g, o)
    h, p = min_max(h, p)

    a, i = min_max(a, i)
    b, f = min_max(b, f)
    c, g = min_max(c, g)
    d, h = min_max(d, h)
    j, n = min_max(j, n)
    k, o = min_max(k, o)
    l, p = min_max(l, p)

    i, q = min_max(i, q)
    a, e = min_max(a, e)
    f, j = min_max(f, j)
    g, k = min_max(g, k)
    h, l = min_max(h, l)
    b, d = min_max(b, d)
    n, p = min_max(n, p)

    i, m = min_max(i, m)
    e, q = min_max(e, q)
    a, c = min_max(a, c)
    f, h = min_max(f, h)
    j, l = min_max(j, l)

    e, i = min_max(e, i)
    m, q = min_max(m, q)
    d, j = min_max(d, j)
    h, n = min_max(h, n)
    a, b = min_max(a, b)

    e, g = min_max(e, g)
    i, k = min_max(i, k)
    m, o = min_max(m, o)
    c, q = min_max(c, q)
    d, f = min_max(d, f)
    h, j = min_max(h, j)
    l, n = min_max(l, n)

    c, i = min_max(c, i)
    g, m = min_max(g, m)
    k, q = min_max(k, q)

    c, e = min_max(c, e)
    g, i = min_max(g, i)
    k, m = min_max(k, m)
    o, q = min_max(o, q)

    c, d = min_max(c, d)
    e, f = min_max(e, f)
    g, h = min_max(g, h)
    i, j = min_max(i, j)
    k, l = min_max(k, l)
    m, n = min_max(m, n)
    o, p = min_max(o, p)
    b, q = min_max(b, q)

    b, i = min_max(b, i)
    d, k = min_max(d, k)
    f, m = min_max(f, m)
    h, o = min_max(h, o)
    j, q = min_max(j, q)

    b, e = min_max(b, e)
    d, g = min_max(d, g)
    f, i = min_max(f, i)
    h, k = min_max(h, k)
    j, m = min_max(j, m)
    l, o = min_max(l, o)
    n, q = min_max(n, q)

    b, c = min_max(b, c)
    d, e = min_max(d, e)
    f, g = min_max(f, g)
    h, i = min_max(h, i)
    j, k = min_max(j, k)
    l, m = min_max(l, m)
    n, o = min_max(n, o)
    p, q = min_max(p, q)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q
end

#    sort 18 values with 82 minmaxs in 13 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T, q::T, r::T) where {T}

    a, q = min_max(a, q)
    b, r = min_max(b, r)
    c, k = min_max(c, k)
    d, l = min_max(d, l)
    e, m = min_max(e, m)
    f, n = min_max(f, n)
    g, o = min_max(g, o)
    h, p = min_max(h, p)

    a, i = min_max(a, i)
    b, j = min_max(b, j)
    c, g = min_max(c, g)
    d, h = min_max(d, h)
    k, o = min_max(k, o)
    l, p = min_max(l, p)

    i, q = min_max(i, q)
    j, r = min_max(j, r)
    a, e = min_max(a, e)
    b, f = min_max(b, f)
    g, k = min_max(g, k)
    h, l = min_max(h, l)

    i, m = min_max(i, m)
    j, n = min_max(j, n)
    e, q = min_max(e, q)
    f, r = min_max(f, r)
    a, c = min_max(a, c)
    b, d = min_max(b, d)

    e, i = min_max(e, i)
    f, j = min_max(f, j)
    m, q = min_max(m, q)
    n, r = min_max(n, r)
    a, b = min_max(a, b)

    e, g = min_max(e, g)
    f, h = min_max(f, h)
    i, k = min_max(i, k)
    j, l = min_max(j, l)
    m, o = min_max(m, o)
    n, p = min_max(n, p)
    c, q = min_max(c, q)
    d, r = min_max(d, r)

    c, i = min_max(c, i)
    d, j = min_max(d, j)
    g, m = min_max(g, m)
    h, n = min_max(h, n)
    k, q = min_max(k, q)
    l, r = min_max(l, r)

    c, e = min_max(c, e)
    d, f = min_max(d, f)
    g, i = min_max(g, i)
    h, j = min_max(h, j)
    k, m = min_max(k, m)
    l, n = min_max(l, n)
    o, q = min_max(o, q)
    p, r = min_max(p, r)

    c, d = min_max(c, d)
    e, f = min_max(e, f)
    g, h = min_max(g, h)
    i, j = min_max(i, j)
    k, l = min_max(k, l)
    m, n = min_max(m, n)
    o, p = min_max(o, p)
    q, r = min_max(q, r)

    b, q = min_max(b, q)
    d, k = min_max(d, k)
    f, m = min_max(f, m)
    h, o = min_max(h, o)

    b, i = min_max(b, i)
    j, q = min_max(j, q)
    d, g = min_max(d, g)
    h, k = min_max(h, k)
    l, o = min_max(l, o)

    b, e = min_max(b, e)
    f, i = min_max(f, i)
    j, m = min_max(j, m)
    n, q = min_max(n, q)

    b, c = min_max(b, c)
    d, e = min_max(d, e)
    f, g = min_max(f, g)
    h, i = min_max(h, i)
    j, k = min_max(j, k)
    l, m = min_max(l, m)
    n, o = min_max(n, o)
    p, q = min_max(p, q)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r
end

#    sort 19 values with 91 minmaxs in 14 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T, q::T, r::T, s::T) where {T}

    a, q = min_max(a, q)
    b, r = min_max(b, r)
    c, s = min_max(c, s)
    d, l = min_max(d, l)
    e, m = min_max(e, m)
    f, n = min_max(f, n)
    g, o = min_max(g, o)
    h, p = min_max(h, p)

    a, i = min_max(a, i)
    b, j = min_max(b, j)
    c, k = min_max(c, k)
    d, h = min_max(d, h)
    l, p = min_max(l, p)

    i, q = min_max(i, q)
    j, r = min_max(j, r)
    k, s = min_max(k, s)
    a, e = min_max(a, e)
    b, f = min_max(b, f)
    c, g = min_max(c, g)
    h, l = min_max(h, l)

    i, m = min_max(i, m)
    j, n = min_max(j, n)
    k, o = min_max(k, o)
    e, q = min_max(e, q)
    f, r = min_max(f, r)
    g, s = min_max(g, s)
    a, c = min_max(a, c)
    b, d = min_max(b, d)

    e, i = min_max(e, i)
    f, j = min_max(f, j)
    g, k = min_max(g, k)
    m, q = min_max(m, q)
    n, r = min_max(n, r)
    o, s = min_max(o, s)
    a, b = min_max(a, b)

    e, g = min_max(e, g)
    f, h = min_max(f, h)
    i, k = min_max(i, k)
    j, l = min_max(j, l)
    m, o = min_max(m, o)
    n, p = min_max(n, p)
    q, s = min_max(q, s)
    d, r = min_max(d, r)

    c, q = min_max(c, q)
    d, j = min_max(d, j)
    g, m = min_max(g, m)
    h, n = min_max(h, n)
    l, r = min_max(l, r)

    c, i = min_max(c, i)
    k, q = min_max(k, q)
    d, f = min_max(d, f)
    h, j = min_max(h, j)
    l, n = min_max(l, n)
    p, r = min_max(p, r)

    c, e = min_max(c, e)
    g, i = min_max(g, i)
    k, m = min_max(k, m)
    o, q = min_max(o, q)

    c, d = min_max(c, d)
    e, f = min_max(e, f)
    g, h = min_max(g, h)
    i, j = min_max(i, j)
    k, l = min_max(k, l)
    m, n = min_max(m, n)
    o, p = min_max(o, p)
    q, r = min_max(q, r)

    b, q = min_max(b, q)
    d, s = min_max(d, s)
    f, m = min_max(f, m)
    h, o = min_max(h, o)

    b, i = min_max(b, i)
    d, k = min_max(d, k)
    j, q = min_max(j, q)
    l, s = min_max(l, s)

    b, e = min_max(b, e)
    d, g = min_max(d, g)
    f, i = min_max(f, i)
    h, k = min_max(h, k)
    j, m = min_max(j, m)
    l, o = min_max(l, o)
    n, q = min_max(n, q)
    p, s = min_max(p, s)

    b, c = min_max(b, c)
    d, e = min_max(d, e)
    f, g = min_max(f, g)
    h, i = min_max(h, i)
    j, k = min_max(j, k)
    l, m = min_max(l, m)
    n, o = min_max(n, o)
    p, q = min_max(p, q)
    r, s = min_max(r, s)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s
end

#    sort 20 values with 97 minmaxs in 14 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T, q::T, r::T, s::T, t::T) where {T}

    a, q = min_max(a, q)
    b, r = min_max(b, r)
    c, s = min_max(c, s)
    d, t = min_max(d, t)
    e, m = min_max(e, m)
    f, n = min_max(f, n)
    g, o = min_max(g, o)
    h, p = min_max(h, p)

    a, i = min_max(a, i)
    b, j = min_max(b, j)
    c, k = min_max(c, k)
    d, l = min_max(d, l)

    i, q = min_max(i, q)
    j, r = min_max(j, r)
    k, s = min_max(k, s)
    l, t = min_max(l, t)
    a, e = min_max(a, e)
    b, f = min_max(b, f)
    c, g = min_max(c, g)
    d, h = min_max(d, h)

    i, m = min_max(i, m)
    j, n = min_max(j, n)
    k, o = min_max(k, o)
    l, p = min_max(l, p)
    e, q = min_max(e, q)
    f, r = min_max(f, r)
    g, s = min_max(g, s)
    h, t = min_max(h, t)
    a, c = min_max(a, c)
    b, d = min_max(b, d)

    e, i = min_max(e, i)
    f, j = min_max(f, j)
    g, k = min_max(g, k)
    h, l = min_max(h, l)
    m, q = min_max(m, q)
    n, r = min_max(n, r)
    o, s = min_max(o, s)
    p, t = min_max(p, t)
    a, b = min_max(a, b)

    e, g = min_max(e, g)
    f, h = min_max(f, h)
    i, k = min_max(i, k)
    j, l = min_max(j, l)
    m, o = min_max(m, o)
    n, p = min_max(n, p)
    q, s = min_max(q, s)
    r, t = min_max(r, t)

    c, q = min_max(c, q)
    d, r = min_max(d, r)
    g, m = min_max(g, m)
    h, n = min_max(h, n)
    s, t = min_max(s, t)

    c, i = min_max(c, i)
    d, j = min_max(d, j)
    k, q = min_max(k, q)
    l, r = min_max(l, r)

    c, e = min_max(c, e)
    d, f = min_max(d, f)
    g, i = min_max(g, i)
    h, j = min_max(h, j)
    k, m = min_max(k, m)
    l, n = min_max(l, n)
    o, q = min_max(o, q)
    p, r = min_max(p, r)

    c, d = min_max(c, d)
    e, f = min_max(e, f)
    g, h = min_max(g, h)
    i, j = min_max(i, j)
    k, l = min_max(k, l)
    m, n = min_max(m, n)
    o, p = min_max(o, p)
    q, r = min_max(q, r)

    b, q = min_max(b, q)
    d, s = min_max(d, s)
    f, m = min_max(f, m)
    h, o = min_max(h, o)

    b, i = min_max(b, i)
    d, k = min_max(d, k)
    j, q = min_max(j, q)
    l, s = min_max(l, s)

    b, e = min_max(b, e)
    d, g = min_max(d, g)
    f, i = min_max(f, i)
    h, k = min_max(h, k)
    j, m = min_max(j, m)
    l, o = min_max(l, o)
    n, q = min_max(n, q)
    p, s = min_max(p, s)

    b, c = min_max(b, c)
    d, e = min_max(d, e)
    f, g = min_max(f, g)
    h, i = min_max(h, i)
    j, k = min_max(j, k)
    l, m = min_max(l, m)
    n, o = min_max(n, o)
    p, q = min_max(p, q)
    r, s = min_max(r, s)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t
end

#    sort 21 values with 107 minmaxs in 15 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T, q::T, r::T, s::T, t::T,
              u::T) where {T}

    a, q = min_max(a, q)
    b, r = min_max(b, r)
    c, s = min_max(c, s)
    d, t = min_max(d, t)
    e, u = min_max(e, u)
    f, n = min_max(f, n)
    g, o = min_max(g, o)
    h, p = min_max(h, p)

    a, i = min_max(a, i)
    b, j = min_max(b, j)
    c, k = min_max(c, k)
    d, l = min_max(d, l)
    e, m = min_max(e, m)

    i, q = min_max(i, q)
    j, r = min_max(j, r)
    k, s = min_max(k, s)
    l, t = min_max(l, t)
    m, u = min_max(m, u)
    a, e = min_max(a, e)
    b, f = min_max(b, f)
    c, g = min_max(c, g)
    d, h = min_max(d, h)

    i, m = min_max(i, m)
    j, n = min_max(j, n)
    k, o = min_max(k, o)
    l, p = min_max(l, p)
    q, u = min_max(q, u)
    f, r = min_max(f, r)
    g, s = min_max(g, s)
    h, t = min_max(h, t)
    a, c = min_max(a, c)
    b, d = min_max(b, d)

    e, q = min_max(e, q)
    f, j = min_max(f, j)
    g, k = min_max(g, k)
    h, l = min_max(h, l)
    n, r = min_max(n, r)
    o, s = min_max(o, s)
    p, t = min_max(p, t)
    a, b = min_max(a, b)

    e, i = min_max(e, i)
    m, q = min_max(m, q)
    f, h = min_max(f, h)
    j, l = min_max(j, l)
    n, p = min_max(n, p)
    r, t = min_max(r, t)

    e, g = min_max(e, g)
    i, k = min_max(i, k)
    m, o = min_max(m, o)
    q, s = min_max(q, s)
    d, r = min_max(d, r)
    h, n = min_max(h, n)

    c, q = min_max(c, q)
    g, u = min_max(g, u)
    d, j = min_max(d, j)
    l, r = min_max(l, r)

    c, i = min_max(c, i)
    g, m = min_max(g, m)
    k, q = min_max(k, q)
    o, u = min_max(o, u)
    d, f = min_max(d, f)
    h, j = min_max(h, j)
    l, n = min_max(l, n)
    p, r = min_max(p, r)

    c, e = min_max(c, e)
    g, i = min_max(g, i)
    k, m = min_max(k, m)
    o, q = min_max(o, q)
    s, u = min_max(s, u)

    c, d = min_max(c, d)
    e, f = min_max(e, f)
    g, h = min_max(g, h)
    i, j = min_max(i, j)
    k, l = min_max(k, l)
    m, n = min_max(m, n)
    o, p = min_max(o, p)
    q, r = min_max(q, r)
    s, t = min_max(s, t)

    b, q = min_max(b, q)
    d, s = min_max(d, s)
    f, u = min_max(f, u)
    h, o = min_max(h, o)

    b, i = min_max(b, i)
    d, k = min_max(d, k)
    f, m = min_max(f, m)
    j, q = min_max(j, q)
    l, s = min_max(l, s)
    n, u = min_max(n, u)

    b, e = min_max(b, e)
    d, g = min_max(d, g)
    f, i = min_max(f, i)
    h, k = min_max(h, k)
    j, m = min_max(j, m)
    l, o = min_max(l, o)
    n, q = min_max(n, q)
    p, s = min_max(p, s)
    r, u = min_max(r, u)

    b, c = min_max(b, c)
    d, e = min_max(d, e)
    f, g = min_max(f, g)
    h, i = min_max(h, i)
    j, k = min_max(j, k)
    l, m = min_max(l, m)
    n, o = min_max(n, o)
    p, q = min_max(p, q)
    r, s = min_max(r, s)
    t, u = min_max(t, u)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u
end

#    sort 22 values with 114 minmaxs in 15 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T, q::T, r::T, s::T, t::T,
              u::T, v::T) where {T}

    a, q = min_max(a, q)
    b, r = min_max(b, r)
    c, s = min_max(c, s)
    d, t = min_max(d, t)
    e, u = min_max(e, u)
    f, v = min_max(f, v)
    g, o = min_max(g, o)
    h, p = min_max(h, p)

    a, i = min_max(a, i)
    b, j = min_max(b, j)
    c, k = min_max(c, k)
    d, l = min_max(d, l)
    e, m = min_max(e, m)
    f, n = min_max(f, n)

    i, q = min_max(i, q)
    j, r = min_max(j, r)
    k, s = min_max(k, s)
    l, t = min_max(l, t)
    m, u = min_max(m, u)
    n, v = min_max(n, v)
    a, e = min_max(a, e)
    b, f = min_max(b, f)
    c, g = min_max(c, g)
    d, h = min_max(d, h)

    i, m = min_max(i, m)
    j, n = min_max(j, n)
    k, o = min_max(k, o)
    l, p = min_max(l, p)
    q, u = min_max(q, u)
    r, v = min_max(r, v)
    g, s = min_max(g, s)
    h, t = min_max(h, t)
    a, c = min_max(a, c)
    b, d = min_max(b, d)

    e, q = min_max(e, q)
    f, r = min_max(f, r)
    g, k = min_max(g, k)
    h, l = min_max(h, l)
    o, s = min_max(o, s)
    p, t = min_max(p, t)
    a, b = min_max(a, b)

    e, i = min_max(e, i)
    f, j = min_max(f, j)
    m, q = min_max(m, q)
    n, r = min_max(n, r)

    e, g = min_max(e, g)
    f, h = min_max(f, h)
    i, k = min_max(i, k)
    j, l = min_max(j, l)
    m, o = min_max(m, o)
    n, p = min_max(n, p)
    q, s = min_max(q, s)
    r, t = min_max(r, t)

    c, q = min_max(c, q)
    d, r = min_max(d, r)
    g, u = min_max(g, u)
    h, v = min_max(h, v)

    c, i = min_max(c, i)
    d, j = min_max(d, j)
    g, m = min_max(g, m)
    h, n = min_max(h, n)
    k, q = min_max(k, q)
    l, r = min_max(l, r)
    o, u = min_max(o, u)
    p, v = min_max(p, v)

    c, e = min_max(c, e)
    d, f = min_max(d, f)
    g, i = min_max(g, i)
    h, j = min_max(h, j)
    k, m = min_max(k, m)
    l, n = min_max(l, n)
    o, q = min_max(o, q)
    p, r = min_max(p, r)
    s, u = min_max(s, u)
    t, v = min_max(t, v)

    c, d = min_max(c, d)
    e, f = min_max(e, f)
    g, h = min_max(g, h)
    i, j = min_max(i, j)
    k, l = min_max(k, l)
    m, n = min_max(m, n)
    o, p = min_max(o, p)
    q, r = min_max(q, r)
    s, t = min_max(s, t)
    u, v = min_max(u, v)

    b, q = min_max(b, q)
    d, s = min_max(d, s)
    f, u = min_max(f, u)
    h, o = min_max(h, o)

    b, i = min_max(b, i)
    d, k = min_max(d, k)
    f, m = min_max(f, m)
    j, q = min_max(j, q)
    l, s = min_max(l, s)
    n, u = min_max(n, u)

    b, e = min_max(b, e)
    d, g = min_max(d, g)
    f, i = min_max(f, i)
    h, k = min_max(h, k)
    j, m = min_max(j, m)
    l, o = min_max(l, o)
    n, q = min_max(n, q)
    p, s = min_max(p, s)
    r, u = min_max(r, u)

    b, c = min_max(b, c)
    d, e = min_max(d, e)
    f, g = min_max(f, g)
    h, i = min_max(h, i)
    j, k = min_max(j, k)
    l, m = min_max(l, m)
    n, o = min_max(n, o)
    p, q = min_max(p, q)
    r, s = min_max(r, s)
    t, u = min_max(t, u)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v
end

#    sort 23 values with 122 minmaxs in 15 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T, q::T, r::T, s::T, t::T,
              u::T, v::T, w::T) where {T}

    a, q = min_max(a, q)
    b, r = min_max(b, r)
    c, s = min_max(c, s)
    d, t = min_max(d, t)
    e, u = min_max(e, u)
    f, v = min_max(f, v)
    g, w = min_max(g, w)
    h, p = min_max(h, p)

    a, i = min_max(a, i)
    b, j = min_max(b, j)
    c, k = min_max(c, k)
    d, l = min_max(d, l)
    e, m = min_max(e, m)
    f, n = min_max(f, n)
    g, o = min_max(g, o)

    i, q = min_max(i, q)
    j, r = min_max(j, r)
    k, s = min_max(k, s)
    l, t = min_max(l, t)
    m, u = min_max(m, u)
    n, v = min_max(n, v)
    o, w = min_max(o, w)
    a, e = min_max(a, e)
    b, f = min_max(b, f)
    c, g = min_max(c, g)
    d, h = min_max(d, h)

    i, m = min_max(i, m)
    j, n = min_max(j, n)
    k, o = min_max(k, o)
    l, p = min_max(l, p)
    q, u = min_max(q, u)
    r, v = min_max(r, v)
    s, w = min_max(s, w)
    h, t = min_max(h, t)
    a, c = min_max(a, c)
    b, d = min_max(b, d)

    e, q = min_max(e, q)
    f, r = min_max(f, r)
    g, s = min_max(g, s)
    h, l = min_max(h, l)
    p, t = min_max(p, t)
    u, w = min_max(u, w)
    a, b = min_max(a, b)

    e, i = min_max(e, i)
    f, j = min_max(f, j)
    g, k = min_max(g, k)
    m, q = min_max(m, q)
    n, r = min_max(n, r)
    o, s = min_max(o, s)

    e, g = min_max(e, g)
    f, h = min_max(f, h)
    i, k = min_max(i, k)
    j, l = min_max(j, l)
    m, o = min_max(m, o)
    n, p = min_max(n, p)
    q, s = min_max(q, s)
    r, t = min_max(r, t)

    c, q = min_max(c, q)
    d, r = min_max(d, r)
    g, u = min_max(g, u)
    h, v = min_max(h, v)

    c, i = min_max(c, i)
    d, j = min_max(d, j)
    g, m = min_max(g, m)
    h, n = min_max(h, n)
    k, q = min_max(k, q)
    l, r = min_max(l, r)
    o, u = min_max(o, u)
    p, v = min_max(p, v)

    c, e = min_max(c, e)
    d, f = min_max(d, f)
    g, i = min_max(g, i)
    h, j = min_max(h, j)
    k, m = min_max(k, m)
    l, n = min_max(l, n)
    o, q = min_max(o, q)
    p, r = min_max(p, r)
    s, u = min_max(s, u)
    t, v = min_max(t, v)

    c, d = min_max(c, d)
    e, f = min_max(e, f)
    g, h = min_max(g, h)
    i, j = min_max(i, j)
    k, l = min_max(k, l)
    m, n = min_max(m, n)
    o, p = min_max(o, p)
    q, r = min_max(q, r)
    s, t = min_max(s, t)
    u, v = min_max(u, v)

    b, q = min_max(b, q)
    d, s = min_max(d, s)
    f, u = min_max(f, u)
    h, w = min_max(h, w)

    b, i = min_max(b, i)
    d, k = min_max(d, k)
    f, m = min_max(f, m)
    h, o = min_max(h, o)
    j, q = min_max(j, q)
    l, s = min_max(l, s)
    n, u = min_max(n, u)
    p, w = min_max(p, w)

    b, e = min_max(b, e)
    d, g = min_max(d, g)
    f, i = min_max(f, i)
    h, k = min_max(h, k)
    j, m = min_max(j, m)
    l, o = min_max(l, o)
    n, q = min_max(n, q)
    p, s = min_max(p, s)
    r, u = min_max(r, u)
    t, w = min_max(t, w)

    b, c = min_max(b, c)
    d, e = min_max(d, e)
    f, g = min_max(f, g)
    h, i = min_max(h, i)
    j, k = min_max(j, k)
    l, m = min_max(l, m)
    n, o = min_max(n, o)
    p, q = min_max(p, q)
    r, s = min_max(r, s)
    t, u = min_max(t, u)
    v, w = min_max(v, w)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w
end

#    sort 24 values with 127 minmaxs in 15 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T, q::T, r::T, s::T, t::T,
              u::T, v::T, w::T, x::T) where {T}

    a, q = min_max(a, q)
    b, r = min_max(b, r)
    c, s = min_max(c, s)
    d, t = min_max(d, t)
    e, u = min_max(e, u)
    f, v = min_max(f, v)
    g, w = min_max(g, w)
    h, x = min_max(h, x)

    a, i = min_max(a, i)
    b, j = min_max(b, j)
    c, k = min_max(c, k)
    d, l = min_max(d, l)
    e, m = min_max(e, m)
    f, n = min_max(f, n)
    g, o = min_max(g, o)
    h, p = min_max(h, p)

    i, q = min_max(i, q)
    j, r = min_max(j, r)
    k, s = min_max(k, s)
    l, t = min_max(l, t)
    m, u = min_max(m, u)
    n, v = min_max(n, v)
    o, w = min_max(o, w)
    p, x = min_max(p, x)
    a, e = min_max(a, e)
    b, f = min_max(b, f)
    c, g = min_max(c, g)
    d, h = min_max(d, h)

    i, m = min_max(i, m)
    j, n = min_max(j, n)
    k, o = min_max(k, o)
    l, p = min_max(l, p)
    q, u = min_max(q, u)
    r, v = min_max(r, v)
    s, w = min_max(s, w)
    t, x = min_max(t, x)
    a, c = min_max(a, c)
    b, d = min_max(b, d)

    e, q = min_max(e, q)
    f, r = min_max(f, r)
    g, s = min_max(g, s)
    h, t = min_max(h, t)
    u, w = min_max(u, w)
    v, x = min_max(v, x)
    a, b = min_max(a, b)

    e, i = min_max(e, i)
    f, j = min_max(f, j)
    g, k = min_max(g, k)
    h, l = min_max(h, l)
    m, q = min_max(m, q)
    n, r = min_max(n, r)
    o, s = min_max(o, s)
    p, t = min_max(p, t)
    w, x = min_max(w, x)

    e, g = min_max(e, g)
    f, h = min_max(f, h)
    i, k = min_max(i, k)
    j, l = min_max(j, l)
    m, o = min_max(m, o)
    n, p = min_max(n, p)
    q, s = min_max(q, s)
    r, t = min_max(r, t)

    c, q = min_max(c, q)
    d, r = min_max(d, r)
    g, u = min_max(g, u)
    h, v = min_max(h, v)

    c, i = min_max(c, i)
    d, j = min_max(d, j)
    g, m = min_max(g, m)
    h, n = min_max(h, n)
    k, q = min_max(k, q)
    l, r = min_max(l, r)
    o, u = min_max(o, u)
    p, v = min_max(p, v)

    c, e = min_max(c, e)
    d, f = min_max(d, f)
    g, i = min_max(g, i)
    h, j = min_max(h, j)
    k, m = min_max(k, m)
    l, n = min_max(l, n)
    o, q = min_max(o, q)
    p, r = min_max(p, r)
    s, u = min_max(s, u)
    t, v = min_max(t, v)

    c, d = min_max(c, d)
    e, f = min_max(e, f)
    g, h = min_max(g, h)
    i, j = min_max(i, j)
    k, l = min_max(k, l)
    m, n = min_max(m, n)
    o, p = min_max(o, p)
    q, r = min_max(q, r)
    s, t = min_max(s, t)
    u, v = min_max(u, v)

    b, q = min_max(b, q)
    d, s = min_max(d, s)
    f, u = min_max(f, u)
    h, w = min_max(h, w)

    b, i = min_max(b, i)
    d, k = min_max(d, k)
    f, m = min_max(f, m)
    h, o = min_max(h, o)
    j, q = min_max(j, q)
    l, s = min_max(l, s)
    n, u = min_max(n, u)
    p, w = min_max(p, w)

    b, e = min_max(b, e)
    d, g = min_max(d, g)
    f, i = min_max(f, i)
    h, k = min_max(h, k)
    j, m = min_max(j, m)
    l, o = min_max(l, o)
    n, q = min_max(n, q)
    p, s = min_max(p, s)
    r, u = min_max(r, u)
    t, w = min_max(t, w)

    b, c = min_max(b, c)
    d, e = min_max(d, e)
    f, g = min_max(f, g)
    h, i = min_max(h, i)
    j, k = min_max(j, k)
    l, m = min_max(l, m)
    n, o = min_max(n, o)
    p, q = min_max(p, q)
    r, s = min_max(r, s)
    t, u = min_max(t, u)
    v, w = min_max(v, w)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x
end

#    sort 25 values with 138 minmaxs in 15 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T, q::T, r::T, s::T, t::T,
              u::T, v::T, w::T, x::T, y::T) where {T}

    a, q = min_max(a, q)
    b, r = min_max(b, r)
    c, s = min_max(c, s)
    d, t = min_max(d, t)
    e, u = min_max(e, u)
    f, v = min_max(f, v)
    g, w = min_max(g, w)
    h, x = min_max(h, x)
    i, y = min_max(i, y)

    a, i = min_max(a, i)
    b, j = min_max(b, j)
    c, k = min_max(c, k)
    d, l = min_max(d, l)
    e, m = min_max(e, m)
    f, n = min_max(f, n)
    g, o = min_max(g, o)
    h, p = min_max(h, p)
    q, y = min_max(q, y)

    i, q = min_max(i, q)
    j, r = min_max(j, r)
    k, s = min_max(k, s)
    l, t = min_max(l, t)
    m, u = min_max(m, u)
    n, v = min_max(n, v)
    o, w = min_max(o, w)
    p, x = min_max(p, x)
    a, e = min_max(a, e)
    b, f = min_max(b, f)
    c, g = min_max(c, g)
    d, h = min_max(d, h)

    i, m = min_max(i, m)
    j, n = min_max(j, n)
    k, o = min_max(k, o)
    l, p = min_max(l, p)
    q, u = min_max(q, u)
    r, v = min_max(r, v)
    s, w = min_max(s, w)
    t, x = min_max(t, x)
    a, c = min_max(a, c)
    b, d = min_max(b, d)

    e, q = min_max(e, q)
    f, r = min_max(f, r)
    g, s = min_max(g, s)
    h, t = min_max(h, t)
    m, y = min_max(m, y)
    v, x = min_max(v, x)
    a, b = min_max(a, b)

    e, i = min_max(e, i)
    f, j = min_max(f, j)
    g, k = min_max(g, k)
    h, l = min_max(h, l)
    m, q = min_max(m, q)
    n, r = min_max(n, r)
    o, s = min_max(o, s)
    p, t = min_max(p, t)
    u, y = min_max(u, y)

    e, g = min_max(e, g)
    f, h = min_max(f, h)
    i, k = min_max(i, k)
    j, l = min_max(j, l)
    m, o = min_max(m, o)
    n, p = min_max(n, p)
    q, s = min_max(q, s)
    r, t = min_max(r, t)
    u, w = min_max(u, w)

    c, q = min_max(c, q)
    d, r = min_max(d, r)
    g, u = min_max(g, u)
    h, v = min_max(h, v)
    k, y = min_max(k, y)

    c, i = min_max(c, i)
    d, j = min_max(d, j)
    g, m = min_max(g, m)
    h, n = min_max(h, n)
    k, q = min_max(k, q)
    l, r = min_max(l, r)
    o, u = min_max(o, u)
    p, v = min_max(p, v)
    s, y = min_max(s, y)

    c, e = min_max(c, e)
    d, f = min_max(d, f)
    g, i = min_max(g, i)
    h, j = min_max(h, j)
    k, m = min_max(k, m)
    l, n = min_max(l, n)
    o, q = min_max(o, q)
    p, r = min_max(p, r)
    s, u = min_max(s, u)
    t, v = min_max(t, v)
    w, y = min_max(w, y)

    c, d = min_max(c, d)
    e, f = min_max(e, f)
    g, h = min_max(g, h)
    i, j = min_max(i, j)
    k, l = min_max(k, l)
    m, n = min_max(m, n)
    o, p = min_max(o, p)
    q, r = min_max(q, r)
    s, t = min_max(s, t)
    u, v = min_max(u, v)
    w, x = min_max(w, x)

    b, q = min_max(b, q)
    d, s = min_max(d, s)
    f, u = min_max(f, u)
    h, w = min_max(h, w)
    j, y = min_max(j, y)

    b, i = min_max(b, i)
    d, k = min_max(d, k)
    f, m = min_max(f, m)
    h, o = min_max(h, o)
    j, q = min_max(j, q)
    l, s = min_max(l, s)
    n, u = min_max(n, u)
    p, w = min_max(p, w)
    r, y = min_max(r, y)

    b, e = min_max(b, e)
    d, g = min_max(d, g)
    f, i = min_max(f, i)
    h, k = min_max(h, k)
    j, m = min_max(j, m)
    l, o = min_max(l, o)
    n, q = min_max(n, q)
    p, s = min_max(p, s)
    r, u = min_max(r, u)
    t, w = min_max(t, w)
    v, y = min_max(v, y)

    b, c = min_max(b, c)
    d, e = min_max(d, e)
    f, g = min_max(f, g)
    h, i = min_max(h, i)
    j, k = min_max(j, k)
    l, m = min_max(l, m)
    n, o = min_max(n, o)
    p, q = min_max(p, q)
    r, s = min_max(r, s)
    t, u = min_max(t, u)
    v, w = min_max(v, w)
    x, y = min_max(x, y)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y
end


text:
const COLOR_SCHEME = Dict(:intrinsic => :blue, :routine => :blue)

function print_routine(io::IO, x::IntrinsicRoutine)
    print(io, routine_name(x))
    fnames = fieldnames(typeof(x))
    if !isempty(fnames)
        print(io, "(")
        for k in 1:length(fnames)
            print(io, getfield(x, fnames[k]))
            if k != length(fnames)
                print(io, ", ")
            end
        end
        print(io, ")")
    end
end

Base.show(io::IO, x::IntrinsicRoutine) = print_routine(io, x)

function Base.show(io::IO, ::MIME"text/plain", x::IntrinsicRoutine)
    print_routine(io, x)
    printstyled(io, " (intrinsic operation)"; color = COLOR_SCHEME[:intrinsic])
end

function Base.show(io::IO, ::GenericRoutine{name}) where {name}
    print(io, name)
end

function Base.show(io::IO, ::MIME"text/plain", fn::GenericRoutine{name}) where {name}
    print(io, name)
    printstyled(
        io,
        " (generic routine with ",
        length(methods(fn).ms),
        " methods)";
        color = COLOR_SCHEME[:routine],
    )
end

function Base.show(io::IO, ::Type{T}) where {name,T<:GenericRoutine{name}}
    mod = T.name.module
    print(io, "typeof(")
    mod === Main || print(io, mod, ".")
    print(io, name, ")")
end

function Base.show(io::IO, op::Operation)
    print(io, routine_name(op.parent))
    print(io, "(")
    join(io, op.args, ", ")
    print(io, ")")
end

function Base.show(io::IO, ::MIME"text/plain", g::Gate)
    print(io, g.locations, " => ", g.operation)
end

Base.show(io::IO, g::Chain) = show(io, MIME"text/plain"(), g)

function Base.show(io::IO, ::MIME"text/plain", g::Chain)
    print(io, "Chain(")
    join(io, g.args, ", ")
    print(io, ")")
end

ext: jl


text:
module Tweaks

# Submodules
include("BiMaps.jl")
include("GraphUtil.jl")
include("Time.jl")
include("Plot.jl")

# Tweaks
include("BaseTweaks.jl"); using .BaseTweaks
include("PyCallTweaks.jl"); using .PyCallTweaks

end


text:
using SliceMap
using Test
using ForwardDiff, Tracker, Zygote, JuliennedArrays

Zygote.refresh()

@testset "columns" begin

    mat = rand(1:9, 3,10)
    fun(x) = 2 .+ x.^2
    res = mapslices(fun, mat, dims=1)

    @test res ≈ mapcols(fun, mat)
    @test res ≈ MapCols{3}(fun, mat)
    @test res ≈ MapCols(fun, mat)

    @test res ≈ tmapcols(fun, mat)
    @test res ≈ ThreadMapCols{3}(fun, mat)
    @test res ≈ ThreadMapCols(fun, mat)

    grad = ForwardDiff.gradient(m -> sum(sin, mapslices(fun, m, dims=1)), mat)

    @test grad ≈ Tracker.gradient(m -> sum(sin, mapcols(fun, m)), mat)[1]
    @test grad ≈ Tracker.gradient(m -> sum(sin, MapCols{3}(fun, m)), mat)[1]
    @test grad ≈ Tracker.gradient(m -> sum(sin, MapCols(fun, m)), mat)[1]

    @test grad ≈ Tracker.gradient(m -> sum(sin, tmapcols(fun, m)), mat)[1]
    @test grad ≈ Tracker.gradient(m -> sum(sin, ThreadMapCols{3}(fun, m)), mat)[1]
    @test grad ≈ Tracker.gradient(m -> sum(sin, ThreadMapCols(fun, m)), mat)[1]

    @test grad ≈ Zygote.gradient(m -> sum(sin, mapcols(fun, m)), mat)[1]
    @test grad ≈ Zygote.gradient(m -> sum(sin, MapCols{3}(fun, m)), mat)[1]

    @test grad ≈ Zygote.gradient(m -> sum(sin, tmapcols(fun, m)), mat)[1]
    @test grad ≈ Zygote.gradient(m -> sum(sin, ThreadMapCols{3}(fun, m)), mat)[1]

    jcols(f,m) = Align(map(f, Slices(m, True(), False())), True(), False())
    @test res ≈ jcols(fun, mat)
    @test grad ≈ Zygote.gradient(m -> sum(sin, jcols(fun, m)), mat)[1]

    # Case of length 1, was an error for _threadmap
    @test res[:,1:1] ≈ mapcols(fun, mat[:, 1:1])
    @test res[:,1:1] ≈ tmapcols(fun, mat[:, 1:1])
    @test res[:,1:1] ≈ ThreadMapCols(fun, mat[:, 1:1])

end
@testset "columns -> scalar" begin

    mat = rand(1:9, 3,10)
    fun(x) = sum(x) # different function!
    res = mapslices(fun, mat, dims=1)

    @test res ≈ mapcols(fun, mat)
    @test res ≈ MapCols{3}(fun, mat)

    grad = ForwardDiff.gradient(m -> sum(sin, mapslices(fun, m, dims=1)), mat)

    @test grad ≈ Tracker.gradient(m -> sum(sin, mapcols(fun, m)), mat)[1]
    @test grad ≈ Tracker.gradient(m -> sum(sin, MapCols{3}(fun, m)), mat)[1]

    @test grad ≈ Zygote.gradient(m -> sum(sin, mapcols(fun, m)), mat)[1]
    @test grad ≈ Zygote.gradient(m -> sum(sin, MapCols{3}(fun, m)), mat)[1]

end
@testset "columns -> matrix" begin

    mat = rand(1:9, 3,10)
    fun(x) = x .* x' # different function! vector -> matrix
    res = mapslices(vec∘fun, mat, dims=1)

    @test res ≈ mapcols(fun, mat)
    @test res ≈ MapCols{3}(fun, mat)

    grad = ForwardDiff.gradient(m -> sum(sin, mapslices(vec∘fun, m, dims=1)), mat)

    @test grad ≈ Tracker.gradient(m -> sum(sin, mapcols(fun, m)), mat)[1]
    @test grad ≈ Tracker.gradient(m -> sum(sin, MapCols{3}(fun, m)), mat)[1]

    @test grad ≈ Zygote.gradient(m -> sum(sin, mapcols(fun, m)), mat)[1]
    @test grad ≈ Zygote.gradient(m -> sum(sin, MapCols{3}(fun, m)), mat)[1]

end
@testset "columns w args" begin

    mat = randn(Float32, 3,10)
    fun(x, s) = 1 .+ x .* s
    res = mapslices(x -> vec(fun(x,5)), mat, dims=1)

    @test res ≈ mapcols(fun, mat, 5)
    @test res ≈ MapCols{3}(fun, mat, 5)

    grad = ForwardDiff.gradient(m -> sum(sin, mapslices(x -> vec(fun(x,5)), m, dims=1)), mat)

    @test grad ≈ Tracker.gradient(m -> sum(sin, mapcols(fun, m, 5)), mat)[1]
    @test grad ≈ Tracker.gradient(m -> sum(sin, MapCols{3}(fun, m, 5)), mat)[1]

    @test grad ≈ Zygote.gradient(m -> sum(sin, mapcols(fun, m, 5)), mat)[1]
    @test grad ≈ Zygote.gradient(m -> sum(sin, MapCols{3}(fun, m, 5)), mat)[1]

end
@testset "rows" begin

    mat = randn(4,5)
    fun(x) = 2 .+ x.^2 ./ sum(x)

    res = mapslices(fun, mat, dims=2)
    @test res ≈ maprows(fun, mat)

    grad = ForwardDiff.gradient(m -> sum(sin, mapslices(fun, m, dims=2)), mat)
    @test grad ≈ Tracker.gradient(m -> sum(sin, maprows(fun, m)), mat)[1]
    @test grad ≈ Zygote.gradient(m -> sum(sin, maprows(fun, m)), mat)[1]

    jrows(f,m) = Align(map(f, Slices(m, False(), True())), False(), True())
    @test res ≈ jrows(fun, mat)
    @test grad ≈ Zygote.gradient(m -> sum(sin, jrows(fun, m)), mat)[1]

end
@testset "slices of a 4-tensor" begin

    ten = randn(3,4,5,2)
    fun(x::AbstractVector) = sqrt(3) .+ x.^3 ./ (sum(x)^2)
    res = mapslices(fun, ten, dims=3)

    @test res ≈ slicemap(fun, ten, dims=3)

    grad = ForwardDiff.gradient(x -> sum(sin, slicemap(fun, x, dims=3)), ten)
    @test grad ≈ Zygote.gradient(x -> sum(sin, slicemap(fun, x, dims=3)), ten)[1]

    jthree(f,m) = Align(map(f,
        Slices(m, False(), False(), True(), False())
        ), False(), False(), True(), False())
    @test res ≈ jthree(fun, ten)
    @test grad ≈ Zygote.gradient(m -> sum(sin, jthree(fun, m)), ten)[1]

    j3(f,m) = Align(map(f, Slices(m, 3)), 3)
    @test res ≈ j3(fun, ten)
    @test grad ≈ Zygote.gradient(m -> sum(sin, j3(fun, m)), ten)[1]

end
@testset "gradient of the function" begin

    struct F W end
    (f::F)(x) = f.W * x # toy version of e.g. Flux.Dense
    w = rand(3,2)
    x = rand(2,5)
    gradx = ForwardDiff.gradient(x -> sum(mapslices(F(w), x, dims=1)), x)
    gradw = ForwardDiff.gradient(w -> sum(mapslices(F(w), x, dims=1)), w)

    wp = Tracker.param(w)
    xp = Tracker.param(x)
    Tracker.back!(sum(mapcols(F(wp), xp)))
    @test Tracker.grad(xp) ≈ gradx
    @test Tracker.grad(wp) == 0 .* gradw # bug or a feature?

    # fp = F(wp)
    # wp.grad .= 0; xp.grad .= 0;
    # Tracker.back!(sum(mapcols(fp, xp)))
    # @test Tracker.grad(xp) ≈ gradx
    # @test_broken Tracker.grad(wp) ≈ gradw # zero

    f = F(w)
    grad_mapcols = Zygote.gradient(() -> sum(mapcols(f, x)), Zygote.Params([w,x]))
    @test grad_mapcols[x] ≈ gradx
    @test grad_mapcols[w] == nothing # bug or a feature?

    grad_slicemap = Zygote.gradient(() -> sum(slicemap(f, x, dims=1)), Zygote.Params([w,x]))
    @test grad_slicemap[x] ≈ gradx
    @test grad_slicemap[w] ≈ gradw
    @test gradw ≈ Zygote.gradient(w -> sum(slicemap(F(w), x, dims=1)), w)[1]
    # Using F(w) with Params() gives wrong answers:
    # https://github.com/FluxML/Zygote.jl/issues/522#issuecomment-605935652

end
@testset "slicemap rev=true" begin

    rec(store) = x -> (push!(store, first(x)); x)
    A = [1,1,1] .* (1:5)'

    store = []
    slicemap(rec(store), A; dims=1)
    @test store == 1:5

    store = []
    slicemap(rec(store), A; dims=1, rev=true)
    @test store == 5:-1:1

    # gradient check as above
    ten = randn(3,4,5,2)
    fun(x::AbstractVector) = sqrt(3) .+ x.^3 ./ (sum(x)^2)
    res = mapslices(fun, ten, dims=3)
    @test res ≈ slicemap(fun, ten; dims=3, rev=true)
    @test res ≈ slicemap(fun, ten; dims=3, rev=false)

    grad = ForwardDiff.gradient(x -> sum(sin, mapslices(fun, x, dims=3)), ten)
    @test grad ≈ Zygote.gradient(x -> sum(sin, slicemap(fun, x, dims=3, rev=false)), ten)[1]
    @test grad ≈ Zygote.gradient(x -> sum(sin, slicemap(fun, x, dims=3, rev=true)), ten)[1]

end


text:
#import Base.Broadcast: Broadcasted
import Base.Broadcast: Broadcasted, BroadcastStyle, ArrayStyle

BroadcastStyle(::Type{<:CuArray}) = ArrayStyle{CuArray}()

function Base.similar(bc::Broadcasted{ArrayStyle{CuArray}}, ::Type{T}) where T
    similar(CuArray, T, length.(axes(bc)))
end

function Base.copyto!(dest::CuArray, bc::Broadcasted{Nothing})
    throw("Not implemented yet.")
end


text:
using PyPlot

## Rescale linearly the dynamic of a vector to fit within a range [a,b]

function rescale(f, a = 0, b = 1
  v = maximum(f) - minimum(f)
  g = copy(f - minimum(f))
  if v > 0
    g = g / v
  end
  return a + g*(b - a)
end

function compute_max(X,d)
    """ compute_max - compute maximum along dimension d

       Y,I = compute_max(X,d);

       Copyright (c) 2008 Gabriel Peyre
    """
    if ndims(X)>=2
        Y = maximum(X,d)
        I = mapslices(indmax,X,d)[:]
    elseif ndims(X)==1
        Y = maximum(X)
        I = indmaw(X)
    end
    return Y,I
end

function compute_max(X)
    if ndims(X)>=2
        Y = maximum(X,1)
        I = mapslices(indmax,X,1)[:]
    elseif ndims(X)==1
        Y = maximum(X)
        I = indmaw(X)
    end
    return Y,I
end


text:
function DrivePOMDP(discount_factor::Float64, Δt::Float64, sRange::Vector{Float64}, vRange::Vector{Float64},
     rRange::Vector{UInt16}, Aset::NamedTuple{(:min, :max, :comfort),Tuple{Float64, Float64, Float64}},
     Routes::Vector{Route}, SsInit::Sts, Stopline::Float64, Smin::Float64)

     println("########## POMDP Model Initialization ###########")

     Δs = sRange[2] - sRange[1]
     Δv = vRange[2] - vRange[1]

     RoutesInit!(Routes)

     for R in Routes
         RouteAVrefset!(R, Δs, Aset)
     end

     SSpace = Sts[]
     EgoSSpace = CarSt[]
     OtherSSpace = CarSt[]

     for v in vRange, s in sRange
         push!(EgoSSpace, CarSt(s, v, SsInit.Ego.r))
     end
     for v in vRange, s in sRange, r in rRange
         push!(OtherSSpace, CarSt(s, v, r))
     end
     for ego in EgoSSpace, other in OtherSSpace
         push!(SSpace, Sts(ego, other))
     end

     OSpace = CarOb[]
     PosVec = Point2D[]
     s0_shared = Routes[rRange[1]].Geos[1].Length
     StdDevs = (0.3, 0.3)
     rng = MersenneTwister(rand(UInt32))
     """
     avoid generating too many duplicate Positions in OSpace because of R_shared
     """
     for s in sRange[1:end]
         if s <= s0_shared
             Point = GetGlobalPosition(Routes[rRange[1]], s)
             x, y = Sample2D(rng, (Point.x, Point.y), StdDevs)
             x = round(x, 0.01) # should we round (x, y)? Obspace and Stspace must not have the same precision. But for searching easily, the answer is currently "yes"
             y = round(y, 0.01)
             push!(PosVec, Point2D{Float64}(x, y))
         end
     end
     for r in rRange, s in sRange[2:end]
         if s > s0_shared
             Point = GetGlobalPosition(Routes[r], s)
             x, y = Sample2D(rng, (Point.x, Point.y), StdDevs)
             x = round(x, 0.01)
             y = round(y, 0.01)
             push!(PosVec, Point2D{Float64}(x, y))
         end
     end
     for v in vRange, Pos in PosVec
         push!(OSpace, CarOb(v, Pos.x, Pos.y))
     end

     DP = DrivePOMDP(discount_factor, Δt, Δs, Δv, sRange, vRange, rRange, Aset, Routes, SSpace, OSpace, SsInit, Stopline, Smin)
     println("###### POMDP Model Initialization Finished ######")
     return DP
 end


text:
function get_prob(outfile = "prob.xlsx")
  url = "https://www.philadelphiafed.org/-/media/research-and-data/real-time-center/survey-of-professional-forecasters/historical-data/prob.xlsx?la=en"
  urldownload(url, parser = identity, save_raw = outfile)
  return nothing
end


text:
"""
    JopRemoveDC(s::JetAbstractSpace{T,N})

F(x) = x - zero frequency component of x along the fastest dimension
Uses Fourier transform, as the DC component differs from the mean as a function of interval length and zero padding.
```
"""
function JopRemoveDC(spc::JetAbstractSpace{T,N}) where {T,N}
    JopLn(dom = spc, rng = spc, df! = JopRemoveDC_df!, df′! = JopRemoveDC_df′!)
end

export JopRemoveDC

function JopRemoveDC_df!(d::AbstractArray{T,2}, m::AbstractArray{T,2}; kwargs...) where {T}
    n1 = size(m,1)
    n2 = div(length(m),n1)
    _m = reshape(m,n1,n2)
    _d = reshape(d,n1,n2)
    for k2 = 1:n2
        mean = sum(_m[:,k2]) / length(_m[:,k2])
        _d[:,k2] .= _m[:,k2] .- mean
    end
    d
end

function JopRemoveDC_df′!(m::AbstractArray{T,2}, d::AbstractArray{T,2}; kwargs...) where {T}
    n1 = size(m,1)
    n2 = div(length(m),n1)
    _m = reshape(m,n1,n2)
    _d = reshape(d,n1,n2)
    for k2 = 1:n2
        mean = sum(_d[:,k2]) / length(_d[:,k2])
        _m[:,k2] .= _d[:,k2] .- mean
    end
    m
end


text:
##############################################################################################
module TestFilterExps
##############################################################################################
##############################################################################################
# imports and exports
using DataAssimilationBenchmarks.FilterExps
using JLD2, Statistics

##############################################################################################
##############################################################################################
# run and analyze the ETKF for state estimation with the Lorenz-96 model

function run_filter_state_L96()
    try
        path = joinpath(@__DIR__, "../src/data/time_series/")
        L96_test_ts = "L96_time_series_seed_0000_dim_40_diff_0.000_F_08.0" *
                      "_tanl_0.05_nanl_05000_spin_1500_h_0.050.jld2"

        args = (path * L96_test_ts, "etkf", 0, 3500, 1.0, 40, 1.00, 21, 1.02)
        filter_state(args)
        true
    catch
        false
    end
end

function analyze_filter_state_L96()
    try
        # test if the filter RMSE for standard simulation falls below adequate threshold
        path = joinpath(@__DIR__, "../src/data/etkf/")
        rmse = load(path * "etkf_L96_state_seed_0000_diff_0.000_sysD_40_obsD_40" *
                    "_obsU_1.00_gamma_001.0_nanl_03500_tanl_0.05_h_0.05_nens_021" *
                    "_stateInfl_1.02.jld2")
        rmse = rmse["filt_rmse"]

        # note, we use a small burn-in to reach more regular cycles
        if mean(rmse[501:end]) < 0.2
            true
        else
            false
        end
    catch
        false
    end
end

##############################################################################################
# run and analyze the ETKF for joint state-parameter estimation with the Lorenz-96 model

function run_filter_param_L96()
    try
        path = joinpath(@__DIR__, "../src/data/time_series/")
        L96_test_ts = "L96_time_series_seed_0000_dim_40_diff_0.000_F_08.0" *
                      "_tanl_0.05_nanl_05000_spin_1500_h_0.050.jld2"

        args = (path * L96_test_ts, "etkf", 0, 3500, 1.0, 40, 1.0, 0.10, 0.0010, 21, 1.02, 1.0)
        filter_param(args)
        true
    catch
        false
    end
end

function analyze_filter_param_L96()
    try
        # test if the filter RMSE for standard simulation falls below adequate threshold
        path = joinpath(@__DIR__, "../src/data/etkf/")
        rmse = load(path * "etkf_L96_param_seed_0000_diff_0.000_sysD_41_stateD_40_obsD_40_" *
                    "obsU_1.00_gamma_001.0_paramE_0.10_paramW_0.0010_nanl_03500_tanl_0.05_" *
                    "h_0.05_nens_021_stateInfl_1.02_paramInfl_1.00.jld2")
        filt_rmse = rmse["filt_rmse"]
        para_rmse = rmse["param_rmse"]

        # note, we use a small burn-in to reach more regular cycles
        if (mean(filt_rmse[501:end]) < 0.2) && (mean(para_rmse[501:end]) < 0.01)
            true
        else
            false
        end
    catch
        false
    end
end

##############################################################################################
# run and analyzed the ETKF for state estimateion with the IEEE39bus model

# static version for test cases
function run_filter_state_IEEE39bus()
    try
        path = joinpath(@__DIR__, "../src/data/time_series/")
        IEEE39bus_test_ts = "IEEE39bus_time_series_seed_0000_diff_0.000_tanl" *
                            "_0.01_nanl_05000_spin_1500_h_0.010.jld2"
        args = (path * IEEE39bus_test_ts, "etkf", 0, 3500, 0.1, 20, 1.00, 21, 1.02)
        filter_state(args)
        true
    catch
        false
    end
end

function analyze_filter_state_IEEE39bus()
    try
        # test if the filter RMSE for standard simulation falls below adequate threshold
        path = joinpath(@__DIR__, "../src/data/etkf/")
        rmse = load(path * "etkf_IEEE39bus_state_seed_0000_diff_0.000_sysD_20_obsD_20_" *
                    "obsU_0.10_gamma_001.0_nanl_03500_tanl_0.01_h_0.01_nens_021_" *
                    "stateInfl_1.02.jld2")
        rmse = rmse["filt_rmse"]

        # note, we use a small burn-in to reach more regular cycles
        if mean(rmse[501:end]) < 0.02
            true
        else
            false
        end
    catch
        false
    end
end

##############################################################################################
# end module

end

max_stars_repo_name: sukhreens/DataAssimilationBenchmarks

max_issues_repo_name: sukhreens/DataAssimilationBenchmarks

max_forks_repo_name: sukhreens/DataAssimilationBenchmarks


text:
PIDNames = Dict(
    1 => "d",
    -1 => "d̄",
    2 => "u",
    -2 => "ū",
    3 => "s",
    -3 => "s̄",
    4 => "c",
    -4 => "c̄",
    5 => "b",
    -5 => "b̄",
    6 => "t",
    -6 => "t̄",
    11 => "e⁻",
    -11 => "e⁺",
    12 => "νₑ",
    -12 => "ν̄ₑ",
    13 => "μ⁻",
    -13 => "μ⁺",
    14 => "ν(μ)",
    -14 => "ν̄(μ)",
    15 => "τ⁻",
    -15 => "τ⁺",
    16 => "ν(τ)",
    -16 => "ν̄(τ)",
    17 => "τ'",
    -17 => "τ'",
    18 => "ν(τ')",
    -18 => "ν̄(τ')",
    21 => "g",
    22 => "γ",
    23 => "Z⁰",
    24 => "W⁺",
    -24 => "W⁻",
    25 => "H⁰",
    111 => "π⁰",
    113 => "ρ(770)⁰",
    115 => "a₂(1320)⁰",
    117 => "ρ₃(1690)⁰",
    119 => "a₄(1970)⁰",
    130 => "K(L)",
    211 => "π⁺",
    -211 => "π⁻",
    213 => "ρ(770)⁺",
    -213 => "ρ(770)⁻",
    215 => "a₂(1320)⁺",
    -215 => "a₂(1320)⁻",
    217 => "ρ₃(1690)⁺",
    -217 => "ρ₃(1690)⁻",
    219 => "a₄(1970)⁺",
    -219 => "a₄(1970)⁻",
    221 => "η",
    223 => "ω(782)",
    225 => "f₂(1270)",
    227 => "ω₃(1670)",
    229 => "f₄(2050)",
    310 => "K(S)",
    311 => "K⁰",
    -311 => "K̄⁰",
    313 => "K*(892)",
    -313 => "K̄*(892)",
    315 => "K(2)*(1430)",
    -315 => "K̄(2)*(1430)",
    317 => "K(3)*(1780)",
    -317 => "K̄(3)*(1780)",
    319 => "K(4)*(2045)",
    -319 => "K̄(4)*(2045)",
    321 => "K⁺",
    -321 => "K⁻",
    323 => "K*(892)",
    -323 => "K̄*(892)",
    325 => "K(2)*(1430)",
    -325 => "K̄(2)*(1430)",
    327 => "K(3)*(1780)",
    -327 => "K̄(3)*(1780)",
    329 => "K(4)*(2045)",
    -329 => "K̄(4)*(2045)",
    331 => "eta'(958)",
    333 => "ϕ(1020)",
    335 => "f(2)'(1525)",
    337 => "ϕ₃(1850)",
    411 => "D⁺",
    -411 => "D⁻",
    413 => "D*(2010)",
    -413 => "D̄*(2010)",
    415 => "D(2)*(2460)",
    -415 => "D̄(2)*(2460)",
    421 => "D⁰",
    -421 => "D̄⁰",
    423 => "D*(2007)",
    -423 => "D̄*(2007)",
    425 => "D(2)*(2460)",
    -425 => "D̄(2)*(2460)",
    431 => "Dₛ⁺",
    -431 => "Dₛ⁻",
    433 => "D(s)*",
    -433 => "D̄(s)*",
    435 => "D(s2)*(2573)",
    -435 => "D̄(s2)*(2573)",
    441 => "eta(c)(1S)",
    443 => "J/ψ(1S)",
    445 => "χ(c2)(1P)",
    511 => "B⁰",
    -511 => "B̄⁰",
    513 => "B*",
    -513 => "B*",
    515 => "B(2)*(5747)",
    -515 => "B̄(2)*(5747)",
    521 => "B⁺",
    -521 => "B⁻",
    523 => "B*",
    -523 => "B̄*",
    525 => "B(2)*(5747)",
    -525 => "B̄(2)*(5747)",
    531 => "Bₛ⁰",
    -531 => "B̄ₛ⁰",
    533 => "B(s)*",
    -533 => "B̄(s)*",
    535 => "B(s2)*(5840)",
    -535 => "B̄(s2)*(5840)",
    541 => "B(c)",
    -541 => "B̄(c)",
    553 => "Υ(1S)",
    555 => "χ(b2)(1P)",
    1103 => "(dd)₁",
    -1103 => "(dd)₁",
    1112 => "Δ(1620)⁻",
    -1112 => "Δ̄(1620)⁺",
    1114 => "Δ(1232)⁻",
    -1114 => "Δ̄(1232)⁺",
    1116 => "Δ(1905)⁻",
    -1116 => "Δ̄(1905)⁺",
    1118 => "Δ(1950)⁻",
    -1118 => "Δ̄(1950)⁺",
    1212 => "Δ(1620)⁰",
    -1212 => "Δ̄(1620)⁰",
    1214 => "N(1520)⁰",
    -1214 => "N̄(1520)⁰",
    1216 => "Δ(1905)⁰",
    -1216 => "Δ̄(1905)⁰",
    1218 => "N(2190)⁰",
    -1218 => "N̄(2190)⁰",
    2101 => "(ud)₀",
    -2101 => "(ud)₀",
    2103 => "(ud)₁",
    -2103 => "(ud)₁",
    2112 => "n",
    -2112 => "n̄",
    2114 => "Δ(1232)⁰",
    -2114 => "Δ̄(1232)⁰",
    2116 => "N(1675)⁰",
    -2116 => "N̄(1675)⁰",
    2118 => "Δ(1950)⁰",
    -2118 => "Δ̄(1950)⁰",
    2122 => "Δ(1620)⁺",
    -2122 => "Δ̄(1620)⁻",
    2124 => "N(1520)⁺",
    -2124 => "N̄(1520)⁻",
    2126 => "Δ(1905)⁺",
    -2126 => "Δ̄(1905)⁻",
    2128 => "N(2190)⁺",
    -2128 => "N̄(2190)⁻",
    2203 => "(uu)₁",
    -2203 => "(uu)₁",
    2212 => "p",
    -2212 => "p̄",
    2214 => "Δ(1232)⁺",
    -2214 => "Δ̄(1232)⁻",
    2216 => "N(1675)⁺",
    -2216 => "N̄(1675)⁻",
    2218 => "Δ(1950)⁺",
    -2218 => "Δ̄(1950)⁻",
    2222 => "Δ(1620)⁺⁺",
    -2222 => "Δ̄(1620)⁻⁻",
    2224 => "Δ(1232)⁺⁺",
    -2224 => "Δ̄(1232)⁻⁻",
    2226 => "Δ(1905)⁺⁺",
    -2226 => "Δ̄(1905)⁻⁻",
    2228 => "Δ(1950)⁺⁺",
    -2228 => "Δ̄(1950)⁻⁻",
    3101 => "(sd)₀",
    -3101 => "(sd)₀",
    3103 => "(sd)₁",
    -3103 => "(sd)₁",
    3112 => "Σ⁻",
    -3112 => "Σ̄⁺",
    3114 => "Σ(1385)⁻",
    -3114 => "Σ̄(1385)⁺",
    3116 => "Σ(1775)⁻",
    -3116 => "Σ̄(1775)⁺",
    3118 => "Σ(2030)⁻",
    -3118 => "Σ̄(2030)⁺",
    3122 => "Λ",
    -3122 => "Λ̄",
    3124 => "Λ(1520)",
    -3124 => "Λ̄(1520)",
    3126 => "Λ(1820)",
    -3126 => "Λ̄(1820)",
    3128 => "Λ(2100)",
    -3128 => "Λ̄(2100)",
    3201 => "(su)₀",
    -3201 => "(su)₀",
    3203 => "(su)₁",
    -3203 => "(su)₁",
    3212 => "Σ⁰",
    -3212 => "Σ̄⁰",
    3214 => "Σ(1385)⁰",
    -3214 => "Σ̄(1385)⁰",
    3216 => "Σ(1775)⁰",
    -3216 => "Σ̄(1775)⁰",
    3218 => "Σ(2030)⁰",
    -3218 => "Σ̄(2030)⁰",
    3222 => "Σ⁺",
    -3222 => "Σ̄⁻",
    3224 => "Σ(1385)⁺",
    -3224 => "Σ̄(1385)⁻",
    3226 => "Σ(1775)⁺",
    -3226 => "Σ̄(1775)⁻",
    3228 => "Σ(2030)⁺",
    -3228 => "Σ̄(2030)⁻",
    3303 => "(ss)₁",
    -3303 => "(ss)₁",
    3312 => "Ξ⁻",
    -3312 => "Ξ̄⁺",
    3314 => "Ξ(1530)⁻",
    -3314 => "Ξ̄(1530)⁺",
    3322 => "Ξ⁰",
    -3322 => "Ξ̄⁰",
    3324 => "Ξ(1530)⁰",
    -3324 => "Ξ̄(1530)⁰",
    3334 => "Ω⁻",
    -3334 => "Ω̄⁺",
    4101 => "(cd)₀",
    -4101 => "(cd)₀",
    4103 => "(cd)₁",
    -4103 => "(cd)₁",
    4112 => "Σ(c)(2455)",
    -4112 => "Σ̄(c)(2455)",
    4114 => "Σ(c)(2520)",
    -4114 => "Σ̄(c)(2520)",
    4122 => "Λ(c)",
    -4122 => "Λ̄(c)",
    4132 => "Ξ(c)",
    -4132 => "Ξ̄(c)",
    4201 => "(cu)₀",
    -4201 => "(cu)₀",
    4203 => "(cu)₁",
    -4203 => "(cu)₁",
    4212 => "Σ(c)(2455)",
    -4212 => "Σ̄(c)(2455)",
    4214 => "Σ(c)(2520)",
    -4214 => "Σ̄(c)(2520)",
    4222 => "Σ(c)(2455)",
    -4222 => "Σ̄(c)(2455)",
    4224 => "Σ(c)(2520)",
    -4224 => "Σ̄(c)(2520)",
    4232 => "Ξ(c)",
    -4232 => "Ξ̄(c)",
    4301 => "(cs)₀",
    -4301 => "(cs)₀",
    4303 => "(cs)₁",
    -4303 => "(cs)₁",
    4312 => "Ξ(c)'",
    -4312 => "Ξ̄(c)'",
    4314 => "Ξ(c)(2645)",
    -4314 => "Ξ̄(c)(2645)",
    4322 => "Ξ(c)'",
    -4322 => "Ξ̄(c)'",
    4324 => "Ξ(c)(2645)",
    -4324 => "Ξ̄(c)(2645)",
    4332 => "Ω(c)",
    -4332 => "Ω̄(c)",
    4334 => "Ω(c)(2770)",
    -4334 => "Ω̄(c)(2770)",
    4403 => "(cc)₁",
    -4403 => "(cc)₁",
    5101 => "(bd)₀",
    -5101 => "(bd)₀",
    5103 => "(bd)₁",
    -5103 => "(bd)₁",
    5112 => "Σ(b)",
    -5112 => "Σ̄(b)",
    5114 => "Σ(b)*",
    -5114 => "Σ̄(b)*",
    5122 => "Λ(b)",
    -5122 => "Λ̄(b)",
    5132 => "Ξ(b)",
    -5132 => "Ξ̄(b)",
    5201 => "(bu)₀",
    -5201 => "(bu)₀",
    5203 => "(bu)₁",
    -5203 => "(bu)₁",
    5222 => "Σ(b)",
    -5222 => "Σ̄(b)",
    5224 => "Σ(b)*",
    -5224 => "Σ̄(b)*",
    5232 => "Ξ(b)",
    -5232 => "Ξ̄(b)",
    5301 => "(bs)₀",
    -5301 => "(bs)₀",
    5303 => "(bs)₁",
    -5303 => "(bs)₁",
    5332 => "Ω(b)",
    -5332 => "Ω̄(b)",
    5401 => "(bc)₀",
    -5401 => "(bc)₀",
    5403 => "(bc)₁",
    -5403 => "(bc)₁",
    5503 => "(bb)₁",
    -5503 => "(bb)₁",
    10111 => "a₀(1450)⁰",
    10113 => "b₁(1235)⁰",
    10115 => "π₂(1670)⁰",
    10211 => "a₀(1450)⁺",
    -10211 => "a₀(1450)⁻",
    10213 => "b₁(1235)⁺",
    -10213 => "b₁(1235)⁻",
    10215 => "π₂(1670)⁺",
    -10215 => "π₂(1670)⁻",
    10221 => "f₀(1370)",
    10223 => "h₁(1170)",
    10225 => "η₂(1645)",
    10311 => "K(0)*(1430)",
    -10311 => "K̄(0)*(1430)",
    10313 => "K₁(1270)⁰",
    -10313 => "K̄₁(1270)⁰",
    10315 => "K₂(1770)⁰",
    -10315 => "K̄₂(1770)⁰",
    10321 => "K(0)*(1430)",
    -10321 => "K̄(0)*(1430)",
    10323 => "K₁(1270)⁺",
    -10323 => "K₁(1270)⁻",
    10325 => "K₂(1770)⁺",
    -10325 => "K₂(1770)⁻",
    10331 => "f₀(1710)",
    10333 => "h₁(1415)",
    10411 => "D(0)*(2300)",
    -10411 => "D̄(0)*(2300)",
    10421 => "D(0)*(2300)",
    -10421 => "D̄(0)*(2300)",
    10423 => "D₁(2420)⁰",
    -10423 => "D̄₁(2420)⁰",
    10431 => "D(s0)*(2317)",
    -10431 => "D̄(s0)*(2317)",
    10433 => "Dₛ₁(2536)⁺",
    -10433 => "Dₛ₁(2536)⁻",
    10441 => "χ(c0)(1P)",
    10443 => "h(c)(1P)",
    10551 => "χ(b0)(1P)",
    10553 => "h(b)(1P)",
    11112 => "Δ(1900)⁻",
    -11112 => "Δ̄(1900)⁺",
    11114 => "Δ(1700)⁻",
    -11114 => "Δ̄(1700)⁺",
    11116 => "Δ(1930)⁻",
    -11116 => "Δ̄(1930)⁺",
    11212 => "Δ(1900)⁰",
    -11212 => "Δ̄(1900)⁰",
    11216 => "Δ(1930)⁰",
    -11216 => "Δ̄(1930)⁰",
    12112 => "N(1440)⁰",
    -12112 => "N̄(1440)⁰",
    12114 => "Δ(1700)⁰",
    -12114 => "Δ̄(1700)⁰",
    12116 => "N(1680)⁰",
    -12116 => "N̄(1680)⁰",
    12122 => "Δ(1900)⁺",
    -12122 => "Δ̄(1900)⁻",
    12126 => "Δ(1930)⁺",
    -12126 => "Δ̄(1930)⁻",
    12212 => "N(1440)⁺",
    -12212 => "N̄(1440)⁻",
    12214 => "Δ(1700)⁺",
    -12214 => "Δ̄(1700)⁻",
    12216 => "N(1680)⁺",
    -12216 => "N̄(1680)⁻",
    12222 => "Δ(1900)⁺⁺",
    -12222 => "Δ̄(1900)⁻⁻",
    12224 => "Δ(1700)⁺⁺",
    -12224 => "Δ̄(1700)⁻⁻",
    12226 => "Δ(1930)⁺⁺",
    -12226 => "Δ̄(1930)⁻⁻",
    13112 => "Σ(1660)⁻",
    -13112 => "Σ̄(1660)⁺",
    13114 => "Σ(1670)⁻",
    -13114 => "Σ̄(1670)⁺",
    13116 => "Σ(1915)⁻",
    -13116 => "Σ̄(1915)⁺",
    13122 => "Λ(1404)",
    -13122 => "Λ̄(1404)",
    13124 => "Λ(1690)",
    -13124 => "Λ̄(1690)",
    13126 => "Λ(1830)",
    -13126 => "Λ̄(1830)",
    13212 => "Σ(1660)⁰",
    -13212 => "Σ̄(1660)⁰",
    13214 => "Σ(1670)⁰",
    -13214 => "Σ̄(1670)⁰",
    13216 => "Σ(1915)⁰",
    -13216 => "Σ̄(1915)⁰",
    13222 => "Σ(1660)⁺",
    -13222 => "Σ̄(1660)⁻",
    13224 => "Σ(1670)⁺",
    -13224 => "Σ̄(1670)⁻",
    13226 => "Σ(1915)⁺",
    -13226 => "Σ̄(1915)⁻",
    13314 => "Ξ(1820)⁻",
    -13314 => "Ξ̄(1820)⁺",
    13324 => "Ξ(1820)⁰",
    -13324 => "Ξ̄(1820)⁰",
    14122 => "Λ(c)(2595)",
    -14122 => "Λ(c)(2595)",
    20113 => "a₁(1260)⁰",
    20213 => "a₁(1260)⁺",
    -20213 => "a₁(1260)⁻",
    20223 => "f₁(1285)",
    20313 => "K₁(1400)⁰",
    -20313 => "K̄₁(1400)⁰",
    20315 => "K₂(1820)⁰",
    -20315 => "K̄₂(1820)⁰",
    20323 => "K₁(1400)⁺",
    -20323 => "K₁(1400)⁻",
    20325 => "K₂(1820)⁺",
    -20325 => "K₂(1820)⁻",
    20333 => "f₁(1420)",
    20433 => "Dₛ₁(2460)⁺",
    -20433 => "Dₛ₁(2460)⁻",
    20443 => "χ(c1)(1P)",
    20553 => "χ(b1)(1P)",
    20555 => "Υ₂(1D)",
    21112 => "Δ(1910)⁻",
    -21112 => "Δ̄(1910)⁺",
    21114 => "Δ(1920)⁻",
    -21114 => "Δ̄(1920)⁺",
    21212 => "Δ(1910)⁰",
    -21212 => "Δ̄(1910)⁰",
    21214 => "N(1700)⁰",
    -21214 => "N̄(1700)⁰",
    22112 => "Δ(1910)⁰",
    -22112 => "Δ̄(1910)⁰",
    22114 => "Δ(1920)⁰",
    -22114 => "Δ̄(1920)⁰",
    22122 => "Δ(1910)⁺",
    -22122 => "Δ̄(1910)⁻",
    22124 => "N(1700)⁺",
    -22124 => "N̄(1700)⁻",
    22212 => "N(1535)⁺",
    -22212 => "N̄(1535)⁻",
    22214 => "Δ(1920)⁺",
    -22214 => "Δ̄(1920)⁻",
    22222 => "Δ(1910)⁺⁺",
    -22222 => "Δ̄(1910)⁻⁻",
    22224 => "Δ(1920)⁺⁺",
    -22224 => "Δ̄(1920)⁻⁻",
    23112 => "Σ(1750)⁻",
    -23112 => "Σ̄(1750)⁺",
    23114 => "Σ(1940)⁻",
    -23114 => "Σ̄(1940)⁺",
    23122 => "Λ(1600)",
    -23122 => "Λ̄(1600)",
    23124 => "Λ(1890)",
    -23124 => "Λ̄(1890)",
    23126 => "Λ(2110)",
    -23126 => "Λ̄(2110)",
    23212 => "Σ(1750)⁰",
    -23212 => "Σ̄(1750)⁰",
    23214 => "Σ(1940)⁰",
    -23214 => "Σ̄(1940)⁰",
    23222 => "Σ(1750)⁺",
    -23222 => "Σ̄(1750)⁻",
    23224 => "Σ(1940)⁺",
    -23224 => "Σ̄(1940)⁻",
    30113 => "ρ(1700)⁰",
    30213 => "ρ(1700)⁺",
    -30213 => "ρ(1700)⁻",
    30223 => "ω(1650)",
    30313 => "K*(1680)",
    -30313 => "K*(1680)",
    30323 => "K*(1680)",
    -30323 => "K*(1680)",
    30443 => "ψ(3770)",
    31114 => "Δ(1600)⁻",
    -31114 => "Δ̄(1600)⁺",
    31214 => "N(1720)⁰",
    -31214 => "N̄(1720)⁰",
    32112 => "N(1650)⁰",
    -32112 => "N̄(1650)⁰",
    32114 => "Δ(1600)⁰",
    -32114 => "Δ̄(1600)⁰",
    32124 => "N(1720)⁺",
    -32124 => "N̄(1720)⁻",
    32212 => "N(1650)⁺",
    -32212 => "N̄(1650)⁻",
    32214 => "Δ(1600)⁺",
    -32214 => "Δ̄(1600)⁻",
    32224 => "Δ(1600)⁺⁺",
    -32224 => "Δ̄(1600)⁻⁻",
    33122 => "Λ(1670)",
    -33122 => "Λ̄(1670)",
    42112 => "N(1710)⁰",
    -42112 => "N̄(1710)⁰",
    42212 => "N(1710)⁺",
    -42212 => "N̄(1710)⁻",
    43122 => "Λ(1800)",
    -43122 => "Λ̄(1800)",
    53122 => "Λ(1810)",
    -53122 => "Λ̄(1810)",
    100111 => "π(1300)⁰",
    100113 => "ρ(1450)⁰",
    100211 => "π(1300)⁺",
    -100211 => "π(1300)⁻",
    100213 => "ρ(1450)⁺",
    -100213 => "ρ(1450)⁻",
    100221 => "η(1295)",
    100313 => "K*(1410)",
    -100313 => "K̄*(1410)",
    100321 => "K(1460)⁺",
    -100321 => "K(1460)⁻",
    100323 => "K*(1410)",
    -100323 => "K̄*(1410)",
    100331 => "η(1475)",
    100333 => "ϕ(1680)",
    100441 => "eta(c)(2S)",
    100443 => "ψ(2S)",
    100445 => "χ(c2)(3930)",
    100553 => "Υ(2S)",
    100555 => "χ(b2)(2P)",
    103316 => "Ξ(1950)⁻",
    -103316 => "Ξ̄(1950)⁺",
    103326 => "Ξ(1950)⁰",
    -103326 => "Ξ̄(1950)⁰",
    104122 => "Λ(c)(2625)",
    -104122 => "Λ̄(c)(2625)",
    104312 => "Ξ(c)(2815)",
    -104312 => "Ξ̄(c)(2815)",
    104314 => "Ξ(c)(2790)",
    -104314 => "Ξ̄(c)(2790)",
    104322 => "Ξ(c)(2815)",
    -104322 => "Ξ̄(c)(2815)",
    104324 => "Ξ(c)(2790)",
    -104324 => "Ξ̄(c)(2790)",
    110551 => "χ(b0)(2P)",
    120553 => "χ(b1)(2P)",
    200553 => "Υ(3S)",
    203312 => "Ξ(1690)⁻",
    -203312 => "Ξ̄(1690)⁺",
    203316 => "Ξ(2030)⁻",
    -203316 => "Ξ̄(2030)⁺",
    203322 => "Ξ(1690)⁰",
    -203322 => "Ξ̄(1690)⁰",
    203326 => "Ξ(2030)⁰",
    -203326 => "Ξ̄(2030)⁰",
    203338 => "Ω(2250)⁻",
    -203338 => "Ω̄(2250)⁺",
    204126 => "Λ(c)(2880)",
    -204126 => "Λ̄(c)(2880)",
    300553 => "Υ(4S)",
    1000223 => "ω(1420)",
    9000111 => "a₀(980)⁰",
    9000113 => "π₁(1400)⁰",
    9000115 => "a₂(1700)⁰",
    9000211 => "a₀(980)⁺",
    -9000211 => "a₀(980)⁻",
    9000213 => "π₁(1400)⁺",
    -9000213 => "π₁(1400)⁻",
    9000215 => "a₂(1700)⁺",
    -9000215 => "a₂(1700)⁻",
    9000221 => "f₀(500)",
    9000311 => "K(0)*(700)",
    -9000311 => "K(0)*(700)",
    9000321 => "K(0)*(700)",
    -9000321 => "K(0)*(700)",
    9000323 => "K₁(1650)⁺",
    -9000323 => "K₁(1650)⁻",
    9000325 => "K₂(1580)⁺",
    -9000325 => "K₂(1580)⁻",
    9000329 => "K₄(2500)⁺",
    -9000329 => "K₄(2500)⁻",
    9000443 => "ψ(4040)",
    9000553 => "Υ(10860)",
    9010111 => "π(1800)⁰",
    9010113 => "π₁(1600)⁰",
    9010211 => "π(1800)⁺",
    -9010211 => "π(1800)⁻",
    9010213 => "π₁(1600)⁺",
    -9010213 => "π₁(1600)⁻",
    9010221 => "f₀(980)",
    9010315 => "K(2)*(1980)",
    -9010315 => "K̄(2)*(1980)",
    9010321 => "K(1830)⁺",
    -9010321 => "K(1830)⁻",
    9010325 => "K(2)*(1980)",
    -9010325 => "K̄(2)*(1980)",
    9010327 => "K₃(2320)⁺",
    -9010327 => "K₃(2320)⁻",
    9010443 => "ψ(4160)",
    9010553 => "Υ(11020)",
    9020113 => "a₁(1640)⁰",
    9020213 => "a₁(1640)⁺",
    -9020213 => "a₁(1640)⁻",
    9020221 => "η(1405)",
    9020321 => "K(0)*(1950)",
    -9020321 => "K̄(0)*(1950)",
    9020325 => "K₂(2250)⁺",
    -9020325 => "K₂(2250)⁻",
    9020443 => "ψ(4415)",
    9030221 => "f₀(1500)",
    9050225 => "f₂(1950)",
    9060225 => "f₂(2010)",
    9080225 => "f₂(2300)",
    9090225 => "f₂(2340)"
)


text:
# # Create a 3D LaMEM setup for La Palma
#
# ## Aim
# In this tutorial, your will learn how to use real data to create a geodynamic model setup with LaMEM. We will use the data of La Palma, which is a volcanic island that started erupting in mid september 2021.
# LaMEM is a cartesian geodynamic model, which implies that you will have to transfer the data from `GeoData` to `CartData`.
#
#
# ## 1. Load data
# We will use two types of data to create the model
#  1) Topography
#  2) Earthquake locations

# We start with loading the required packages
using GeophysicalModelGenerator, JLD2
using GMT

# In case you managed to install GMT on your machine, you can automatically download the topography with
Topo = ImportTopo(lon = [-18.7, -17.1], lat=[28.0, 29.2], file="@earth_relief_03s.grd")

# In case you did not manage that, we have prepared a JLD2 file [here](https://seafile.rlp.net/f/03286a46a9f040a69b0f/?dl=1).
# Download it, and doublecheck that you are in the same directory as the data file with:
pwd()

# Load the data:
Topo=load("Topo_LaPalma.jld2","Topo")

# Next, lets load the seismicity. We have prepared a file with earthquake locations up to early November (from [https://www.ign.es/web/ign/portal/vlc-catalogo](https://www.ign.es/web/ign/portal/vlc-catalogo)).
# Download that [here](https://seafile.rlp.net/f/03286a46a9f040a69b0f/?dl=1)
data_all_EQ = load("EQ_Data.jld2","data_all_EQ")

# Write the data to paraview with
Write_Paraview(data_all_EQ,"data_all_EQ",PointsData=true)
Write_Paraview(Topo,"Topo")

# As earthquakes are point-wise data, you have to specify this.
# ![LaPalma_EQTopo_GeoData](../assets/img/TopoEQs_LaPalma_GeoData.png)
# Note that this data is not in "easy" coordinates (see coordinate axis in the plot, where `z` is *not* pointing upwards).

# ## 2. Convert data
# In order to create model setups, it is helpful to first transfer the data to Cartesian.
# This requires us to first determine a *projection point*, that is fixed. Often, it is helpful to use the center of the topography for this. In the present example, we will center the model around La Palma itself:
proj = ProjectionPoint(Lon=-17.84, Lat=28.56)

# Once this is done you can convert the topographic data to the cartesian reference frame
EQ_cart   = Convert2CartData(data_all_EQ, proj);
Topo_cart = Convert2CartData(Topo, proj)

# It is important to realize that the cartesian coordinates of the topographic grid is no longer strictly orthogonal after this conversion. You don't notice that in the current example, as the model domain is rather small.
# In other cases, however, this is quite substantial (e.g., India-Asia collision zone).
# LaMEM needs an orthogonal grid of topography, which we can create with:
Topo_LaMEM = CartData(XYZGrid(-70:.2:70,-60:.2:70,0));

# In a next step, the routine `ProjectCartData` projects a `GeoData` structure to a `CartData` struct
Topo_LaMEM = ProjectCartData(Topo_LaMEM, Topo, proj)

# Let's have a look at the data:
Write_Paraview(EQ_cart,"EQ_cart",PointsData=true)
Write_Paraview(Topo_LaMEM,"Topo_LaMEM")

# ## 3. Create LaMEM setup
#
# In a next step, we need to read the LaMEM input file of the simulation. In particular, this will read the lateral dimensions of the grid and
# the number of control volumes (elements), you want to apply in every direction.
# The LaMEM input file can be downloaded [here](../scripts/LaPalma.dat).
# Make sure you are in the same directory as the *.dat file & execute the following command
Grid    = ReadLaMEM_InputFile("LaPalma.dat")

# The `LaMEM_grid` structure contains the number of elements in every direction and the number of markers in every cell.
# It also contains `Grid.X`, `Grid.Y`, `Grid.Z`, which are the coordinates of each of the markers in the 3 directions.
#
# In a next step we need to give each of these points a `Phase` number (which is an integer, that indicates the type of the rock that point has), as well as the temperature (in Celcius).
Phases  = ones(Int32,size(Grid.X))*2;
Temp    = ones(Float64,size(Grid.X));

# In this example we set the temperature based on the depth, assuming a constant geotherm. Depth is given in kilometers
Geotherm =  30;
Temp     =  -Grid.Z.*Geotherm;

# Since we are in La Palma, we assume that temperatures are not less than 20 degrees. Moreover, in the mantle, we have the mantle adiabat (1350 C)
Temp[Temp.<20]    .=  20;
Temp[Temp.>1350]  .=  1350;

# Because of lacking data, we have pretty much no idea where the Moho is in La Palma.
# Somewhat arbitrary, we assume it to be at 40 km depth, and set the rocks below that to "mantle":
ind = findall( Grid.Z .< -40);
Phases[ind] .= 3;

# Everything above the free surface is assumed to be "air"
ind     =   AboveSurface(Grid, Topo_cart);
Phases[ind] .= 0;

# And all "air" points that are below sea-level becomes "water"
ind = findall( (Phases.==0) .& (Grid.Z .< 0));
Phases[ind] .= 1;

# You can interpret the seismicity in different ways. If you assume that there are fully molten magma chambers, there should be no earthquakes within the magma chamber (as stresses are liklely to be very small).
# If, however, this is a partially molten mush, extraction of melt of that mush will change the fluid pressure and may thus release build-up stresses.
# We will assume the second option in our model setup.
#
# Looking at the seismicity, there is a swarm at around 35 km depth
AddSphere!(Phases,Temp,Grid, cen=(0,0,-35), radius=5, phase=ConstantPhase(5), T=ConstantTemp(1200));

# A shallower one exists as well
AddEllipsoid!(Phases,Temp,Grid, cen=(-1,0,-11), axes=(3,3,8), StrikeAngle=225, DipAngle=45, phase=ConstantPhase(5), T=ConstantTemp(1200));

# And there might be a mid-crustal one
AddEllipsoid!(Phases,Temp,Grid, cen=(-0,0,-23), axes=(8,8,2), StrikeAngle=0, DipAngle=0, phase=ConstantPhase(5), T=ConstantTemp(1200));

# We can generate a 3D model setup, which must include the `Phases` and `Temp` arrays
Model3D = CartData(Grid, (Phases=Phases,Temp=Temp))
Write_Paraview(Model3D,"Model3D")

# The model setup looks like this
# ![LaMEM_ModelSetup_LaPalma](../assets/img/LaMEM_ModelSetup_LaPalma.png)

# We can create a LaMEM marker file from the `Model3D` setup and the (cartesian) topography
Save_LaMEMTopography(Topo_cart, "Topography.txt")
Save_LaMEMMarkersParallel(Model3D)

# Next, you can use this to run the LaMEM model and visualize the model results in Paraview as well.
# If you are interested in doing this, have a look at the LaMEM [wiki](https://bitbucket.org/bkaus/lamem/wiki/Home) pages.

#src Note: The markdown page is generated using:
#src Literate.markdown("LaPalma_example.jl",outputdir="../man/", execute=true)


text:
"""
≪(μ,ν)

# Absolute continuity

The following are equivalent:
1. μ ≪ ν
2. μ is absolutely continuous wrt ν
3. There exists a function f such that μ = ∫f dν
"""
function ≪ end

"""
≃(μ,ν)

# Equivalence of Measure

Measures μ and ν on the same space X are equivalent, written `μ ≃ ν`, if `μ ≪ ν` and `ν ≪ μ`.

Note that this is often written `~` in the literature, but this is overloaded in probabilistic programming, so we use alternate notation.
"""
function ≃(μ,ν)
    return (μ≪ν && ν≪μ)
end

export representative

function representative(μ)
    # Check if we're done
    isprimitive(μ) && return μ

    ν = basemeasure(μ)

    # Make sure not to leave the equivalence class
    (ν ≪ μ) || return μ

    # Fall back on a recusive call
    return representative(ν)
end

≪(μ, ν) = representative(μ) ≪ representative(ν)


text:
println("test adaptive mel-cepstrum analysis")

function test_acep_base(f::Function, order, args...; kargs...)
    srand(98765)
    dummy_input = rand(64)

    c = zeros(order + 1)
    @assert applicable(f, c, first(dummy_input), args...)
    for x in dummy_input
        f(c, x, args...; kargs...)
        @test all(isfinite(c))
    end
end

println("-- test_acep")
for order in [20, 22, 24]
    for pd = 4:5
        println("  where order = $order, pd = $pd")
        test_acep_base(acep!, order; pd=pd)
    end
end

let
    println("-- test_acep exceptions")
    c = zeros(21)
    @test_throws ArgumentError acep!(c, 0.0, pd=3)
    @test_throws ArgumentError acep!(c, 0.0, pd=6)
end

println("-- test_agcep")
for order in [20, 22, 24]
    for stage in 1:10
        println("  where order = $order, stage = $stage")
        test_acep_base(agcep!, order, stage)
    end
end

let
    println("-- test_agcep exceptions")
    c = zeros(21)
    @test_throws ArgumentError agcep!(c, 0.0, 0)
    @test_throws ArgumentError agcep!(c, 0.0, -1)
end

println("-- test_amcep")
for order in [20, 22, 24]
    for α in [0.35, 0.41, 0.5]
        for pd = 4:5
            println("  where order = $order, α = $α, pd = $pd")
            test_acep_base(amcep!, order, α, pd=pd)
        end
    end
end

let
    println("-- test_amcep exceptions")
    c = zeros(21)
    @test_throws ArgumentError amcep!(c, 0.0, 0.35, pd=3)
    @test_throws ArgumentError amcep!(c, 0., 0.35, pd=6)
end