text:
# Linear algebra functions for dense matrices in column major format

scale!{T<:BlasFloat}(X::Array{T}, s::Number) = BLAS.scal!(length(X), convert(T,s), X, 1)
scale!{T<:BlasComplex}(X::Array{T}, s::Real) = BLAS.scal!(length(X), oftype(real(zero(T)),s), X, 1)

#Test whether a matrix is positive-definite
isposdef!{T<:BlasFloat}(A::Matrix{T}, UL::Char) = LAPACK.potrf!(UL, A)[2] == 0
isposdef!(A::Matrix) = ishermitian(A) && isposdef!(A, 'U')

isposdef{T<:BlasFloat}(A::Matrix{T}, UL::Char) = isposdef!(copy(A), UL)
isposdef{T<:BlasFloat}(A::Matrix{T}) = isposdef!(copy(A))
isposdef{T<:Number}(A::Matrix{T}, UL::Char) = isposdef!(float64(A), UL)
isposdef{T<:Number}(A::Matrix{T}) = isposdef!(float64(A))
isposdef(x::Number) = imag(x)==0 && real(x) > 0

function norm{T<:BlasFloat, TI<:Integer}(x::StridedVector{T}, rx::Union(UnitRange{TI},Range{TI}))
    (minimum(rx) < 1 || maximum(rx) > length(x)) && throw(BoundsError())
    BLAS.nrm2(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx))
end

vecnorm1{T<:BlasReal}(x::Union(Array{T},StridedVector{T})) = BLAS.asum(x)
vecnorm2{T<:BlasFloat}(x::Union(Array{T},StridedVector{T})) = BLAS.nrm2(x)

function triu!{T}(M::Matrix{T}, k::Integer)
    m, n = size(M)
    idx = 1
    for j = 0:n-1
        ii = min(max(0, j+1-k), m)
        for i = (idx+ii):(idx+m-1)
            M[i] = zero(T)
        end
        idx += m
    end
    M
end

triu(M::Matrix, k::Integer) = triu!(copy(M), k)

function tril!{T}(M::Matrix{T}, k::Integer)
    m, n = size(M)
    idx = 1
    for j = 0:n-1
        ii = min(max(0, j-k), m)
        for i = idx:(idx+ii-1)
            M[i] = zero(T)
        end
        idx += m
    end
    M
end

tril(M::Matrix, k::Integer) = tril!(copy(M), k)

function gradient(F::Vector, h::Vector)
    n = length(F)
    T = typeof(one(eltype(F))/one(eltype(h)))
    g = Array(T,n)
    if n == 1
        g[1] = zero(T)
    elseif n > 1
        g[1] = (F[2] - F[1]) / (h[2] - h[1])
        g[n] = (F[n] - F[n-1]) / (h[end] - h[end-1])
        if n > 2
            h = h[3:n] - h[1:n-2]
            g[2:n-1] = (F[3:n] - F[1:n-2]) ./ h
        end
    end
    g
end

function diagind(m::Integer, n::Integer, k::Integer=0)
    if 0 < k < n
        return range(k*m+1, m+1, min(m, n-k))
    elseif 0 <= -k <= m
        return range(1-k, m+1, min(m+k,n))
    end
    throw(BoundsError())
end

diagind(A::AbstractMatrix, k::Integer=0) = diagind(size(A,1), size(A,2), k)

diag(A::Matrix, k::Integer=0) = A[diagind(A,k)]

function diagm{T}(v::AbstractVector{T}, k::Integer=0)
    n = length(v) + abs(k)
    A = zeros(T,n,n)
    A[diagind(A,k)] = v
    A
end

diagm(x::Number) = (X = Array(typeof(x),1,1); X[1,1] = x; X)

function trace{T}(A::Matrix{T})
    n = chksquare(A)
    t = zero(T)
    for i=1:n
        t += A[i,i]
    end
    t
end

function kron{T,S}(a::Matrix{T}, b::Matrix{S})
    R = Array(promote_type(T,S), size(a,1)*size(b,1), size(a,2)*size(b,2))
    m = 1
    for j = 1:size(a,2), l = 1:size(b,2), i = 1:size(a,1)
        aij = a[i,j]
        for k = 1:size(b,1)
            R[m] = aij*b[k,l]
            m += 1
        end
    end
    R
end

kron(a::Number, b::Union(Number, Vector, Matrix)) = a * b
kron(a::Union(Vector, Matrix), b::Number) = a * b
kron(a::Vector, b::Vector)=vec(kron(reshape(a,length(a),1),reshape(b,length(b),1)))
kron(a::Matrix, b::Vector)=kron(a,reshape(b,length(b),1))
kron(a::Vector, b::Matrix)=kron(reshape(a,length(a),1),b)

^(A::Matrix, p::Integer) = p < 0 ? inv(A^-p) : Base.power_by_squaring(A,p)

function ^(A::Matrix, p::Number)
    isinteger(p) && return A^integer(real(p))

    chksquare(A)
    v, X = eig(A)
    any(v.<0) && (v = complex(v))
    Xinv = ishermitian(A) ? X' : inv(X)
    scale(X, v.^p)*Xinv
end

function rref{T}(A::Matrix{T})
    nr, nc = size(A)
    U = copy!(similar(A, T <: Complex ? Complex128 : Float64), A)
    e = eps(norm(U,Inf))
    i = j = 1
    while i <= nr && j <= nc
        (m, mi) = findmax(abs(U[i:nr,j]))
        mi = mi+i - 1
        if m <= e
            U[i:nr,j] = 0
            j += 1
        else
            for k=j:nc
                U[i, k], U[mi, k] = U[mi, k], U[i, k]
            end
            d = U[i,j]
            for k = j:nc
                U[i,k] /= d
            end
            for k = 1:nr
                if k != i
                    d = U[k,j]
                    for l = j:nc
                        U[k,l] -= d*U[i,l]
                    end
                end
            end
            i += 1
            j += 1
        end
    end
    U
end

rref(x::Number) = one(x)

# Matrix exponential
expm{T<:BlasFloat}(A::StridedMatrix{T}) = expm!(copy(A))
expm{T<:Integer}(A::StridedMatrix{T}) = expm!(float(A))
expm(x::Number) = exp(x)

## Destructive matrix exponential using algorithm from Higham, 2008,
## "Functions of Matrices: Theory and Computation", SIAM
function expm!{T<:BlasFloat}(A::StridedMatrix{T})
    n = chksquare(A)
    n<2 && return exp(A)
    ilo, ihi, scale = LAPACK.gebal!('B', A)    # modifies A
    nA   = norm(A, 1)
    I    = eye(T,n)
    ## For sufficiently small nA, use lower order Padé-Approximations
    if (nA <= 2.1)
        if nA > 0.95
            C = T[17643225600.,8821612800.,2075673600.,302702400.,
                     30270240.,   2162160.,    110880.,     3960.,
                           90.,         1.]
        elseif nA > 0.25
            C = T[17297280.,8648640.,1995840.,277200.,
                     25200.,   1512.,     56.,     1.]
        elseif nA > 0.015
            C = T[30240.,15120.,3360.,
                    420.,   30.,   1.]
        else
            C = T[120.,60.,12.,1.]
        end
        A2 = A * A
        P  = copy(I)
        U  = C[2] * P
        V  = C[1] * P
        for k in 1:(div(size(C, 1), 2) - 1)
            k2 = 2 * k
            P *= A2
            U += C[k2 + 2] * P
            V += C[k2 + 1] * P
        end
        U = A * U
        X = V + U
        LAPACK.gesv!(V-U, X)
    else
        s  = log2(nA/5.4)               # power of 2 later reversed by squaring
        if s > 0
            si = iceil(s)
            A /= oftype(T,2^si)
        end
        CC = T[64764752532480000.,32382376266240000.,7771770303897600.,
                1187353796428800.,  129060195264000.,  10559470521600.,
                    670442572800.,      33522128640.,      1323241920.,
                        40840800.,           960960.,           16380.,
                             182.,                1.]
        A2 = A * A
        A4 = A2 * A2
        A6 = A2 * A4
        U  = A * (A6 * (CC[14]*A6 + CC[12]*A4 + CC[10]*A2) +
                  CC[8]*A6 + CC[6]*A4 + CC[4]*A2 + CC[2]*I)
        V  = A6 * (CC[13]*A6 + CC[11]*A4 + CC[9]*A2) +
                   CC[7]*A6 + CC[5]*A4 + CC[3]*A2 + CC[1]*I

        X = V + U
        LAPACK.gesv!(V-U, X)

        if s > 0            # squaring to reverse dividing by power of 2
            for t=1:si X *= X end
        end
    end

	# Undo the balancing
    for j = ilo:ihi
        scj = scale[j]
        for i = 1:n
            X[j,i] *= scj
        end
        for i = 1:n
            X[i,j] /= scj
        end
    end

    if ilo > 1       # apply lower permutations in reverse order
        for j in (ilo-1):-1:1 rcswap!(j, int(scale[j]), X) end
    end
    if ihi < n       # apply upper permutations in forward order
        for j in (ihi+1):n    rcswap!(j, int(scale[j]), X) end
    end
    X
end

## Swap rows i and j and columns i and j in X
function rcswap!{T<:Number}(i::Integer, j::Integer, X::StridedMatrix{T})
    for k = 1:size(X,1)
        X[k,i], X[k,j] = X[k,j], X[k,i]
    end
    for k = 1:size(X,2)
        X[i,k], X[j,k] = X[j,k], X[i,k]
    end
end

function sqrtm{T<:Real}(A::StridedMatrix{T})
    issym(A) && return sqrtm(Symmetric(A))
    n = chksquare(A)
    SchurF = schurfact(complex(A))
    R = full(sqrtm(Triangular(SchurF[:T])))
    retmat = SchurF[:vectors]*R*SchurF[:vectors]'
    all(imag(retmat) .== 0) ? real(retmat) : retmat
end
function sqrtm{T<:Complex}(A::StridedMatrix{T})
    ishermitian(A) && return sqrtm(Hermitian(A))
    n = chksquare(A)
    SchurF = schurfact(A)
    R = full(sqrtm(Triangular(SchurF[:T])))
    SchurF[:vectors]*R*SchurF[:vectors]'
end
sqrtm(a::Number) = (b = sqrt(complex(a)); imag(b) == 0 ? real(b) : b)
sqrtm(a::Complex) = sqrt(a)

function det(A::Matrix)
    (istriu(A) || istril(A)) && return det(Triangular(A, :U, false))
    return det(lufact(A))
end
det(x::Number) = x

logdet(A::Matrix) = logdet(lufact(A))

function inv(A::Matrix)
    if istriu(A) return inv(Triangular(A, :U, false)) end
    if istril(A) return inv(Triangular(A, :L, false)) end
    return inv(lufact(A))
end

function factorize{T}(A::Matrix{T})
    m, n = size(A)
    if m == n
        if m == 1 return A[1] end
        utri    = true
        utri1   = true
        herm    = true
        sym     = true
        for j = 1:n-1, i = j+1:m
            if utri1
                if A[i,j] != 0
                    utri1 = i == j + 1
                    utri = false
                end
            end
            if sym
                sym &= A[i,j] == A[j,i]
            end
            if herm
                herm &= A[i,j] == conj(A[j,i])
            end
            if !(utri1|herm|sym) break end
        end
        ltri = true
        ltri1 = true
        for j = 3:n, i = 1:j-2
            ltri1 &= A[i,j] == 0
            if !ltri1 break end
        end
        if ltri1
            for i = 1:n-1
                if A[i,i+1] != 0
                    ltri &= false
                    break
                end
            end
            if ltri
                if utri
                    return Diagonal(A)
                end
                if utri1
                    return Bidiagonal(diag(A), diag(A, -1), false)
                end
                return Triangular(A, :L)
            end
            if utri
                return Bidiagonal(diag(A), diag(A, 1), true)
            end
            if utri1
                if (herm & (T <: Complex)) | sym
                    try
                        return ldltfact!(SymTridiagonal(diag(A), diag(A, -1)))
                    end
                end
                return lufact(Tridiagonal(diag(A, -1), diag(A), diag(A, 1)))
            end
        end
        if utri
            return Triangular(A, :U)
        end
        if herm
            try
                return cholfact(A)
            end
            return factorize(Hermitian(A))
        end
        if sym
            return factorize(Symmetric(A))
        end
        return lufact(A)
    end
    qrfact(A,pivot=T<:BlasFloat)
end

(\)(a::Vector, B::StridedVecOrMat) = (\)(reshape(a, length(a), 1), B)
function (\)(A::StridedMatrix, B::StridedVecOrMat)
    m, n = size(A)
    if m == n
        if istril(A)
            return istriu(A) ? \(Diagonal(A),B) : \(Triangular(A, :L),B)
        end
        istriu(A) && return \(Triangular(A, :U),B)
        return \(lufact(A),B)
    end
    return qrfact(A,pivot=eltype(A)<:BlasFloat)\B
end

## Moore-Penrose inverse
function pinv{T}(A::StridedMatrix{T})
    SVD         = svdfact(A, thin=true)
    S           = eltype(SVD[:S])
    m, n        = size(A)
    (m == 0 || n == 0) && return Array(S, n, m)
    Sinv        = zeros(S, length(SVD[:S]))
    index       = SVD[:S] .> eps(real(float(one(T))))*max(m,n)*maximum(SVD[:S])
    Sinv[index] = one(S) ./ SVD[:S][index]
    return SVD[:Vt]'scale(Sinv, SVD[:U]')
end
pinv(a::StridedVector) = pinv(reshape(a, length(a), 1))
pinv(x::Number) = one(x)/x

## Basis for null space
function null{T}(A::StridedMatrix{T})
    m, n = size(A)
    (m == 0 || n == 0) && return eye(T, n)
    SVD = svdfact(A, thin=false)
    indstart = sum(SVD[:S] .> max(m,n)*maximum(SVD[:S])*eps(eltype(SVD[:S]))) + 1
    return SVD[:V][:,indstart:end]
end
null(a::StridedVector) = null(reshape(a, length(a), 1))

function cond(A::StridedMatrix, p::Real=2)
    if p == 2
        v = svdvals(A)
        maxv = maximum(v)
        return maxv == 0.0 ? inf(typeof(real(A[1,1]))) : maxv / minimum(v)
    elseif p == 1 || p == Inf
        chksquare(A)
        return cond(lufact(A), p)
    end
    throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2 or Inf"))
end

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia


text:
tol = 0.5

# test the Three DOF Vehicle Model
include("ThreeDOF.jl")

ext: jl


text:
export contains_expression

# SUBEXPRESSION SEARCH

# if haystack is a logical symbol, the kneedle must be equal to be a subexpression
"""
    contains_expression(haystack::AbstractExpression, kneedle::AbstractExpression)

Searches for a subexpression `kneedle` recursively within the `haystack` expression. Returns true if a match can be found
and false otherwise.
"""
contains_expression(haystack::LogicalSymbol, kneedle::AbstractExpression) = isequal(haystack, kneedle)

# TODO: come up with quicker subexpression algorithm for LogicalExpressions
# if kneedle is an expression, use a recursive equality check
function contains_expression(haystack::LogicalExpression, kneedle::LogicalExpression)
    if isequal(haystack, kneedle)
        return true
    end
    return any(contains_expression.(arguments(haystack), Iterators.repeat([kneedle], length(arguments(haystack)))))
end

# if kneedle is a symbol, use the variables set to quickly check if its a subexpression
contains_expression(haystack::LogicalExpression, kneedle::LogicalSymbol) = kneedle ∈ variables(haystack)

# dfs for expression search; returns path to first found expression, otherwise empty vector
locate_expression(haystack::LogicalSymbol, kneedle::AbstractExpression) = []
function locate_expression(haystack::LogicalExpression, kneedle::AbstractExpression)
    i = 1
    for arg ∈ arguments(haystack)
        if isequal(arg, kneedle)
            return [i]
        end
        located = locate_expression(arg, kneedle)
        if length(located) > 0
            return [i, located...]
        end
        i += 1
    end
    return []
end

function select_by_path(expr::AbstractExpression, path::Vector{Int})
    if length(path) == 0
        return expr
    end
    return select_by_path(arguments(expr)[path[1]], path[2:end])
end


text:
#############################################################################
#############################################################################
#
# This file contains units tests for polynomial factorization
#
#############################################################################
#############################################################################

"""
Test factorization of polynomials.
"""
function factor_test_poly(;N::Int = 25, seed::Int = 0, primes::Vector{Int} = [5,7,17,19])
    Random.seed!(seed)
    for prime in primes
        print("\ndoing prime = $prime \t")
        for _ in 1:N
            print(".")
            p = rand(Polynomial)
            factorization = factor(p, prime)
            pr = mod(expand_factorization(factorization),prime)
            @assert mod(p-pr,prime) == 0
        end
    end

    println("\nfactor_test_poly - PASSED")
end


text:
immutable Config
  bots::Int
  field_radius::Int
  move::Int
  start_hp::Int
  cannon::Int
  radar::Int
  see::Int
  max_count::Int
  loop_time::Int
  function Config(bots::Integer=3, fieldRadius::Integer=14, move=2, startHp=10, cannon=1,
                 radar=3, see=2, maxCount=200, loopTime=300)
    new(bots, fieldRadius, move, startHp, cannon, radar, see, maxCount, loopTime)
  end
end

function Config(dict::Dict)
  params = [dict["bots"], dict["fieldRadius"], dict["move"], dict["startHp"], dict["cannon"], dict["radar"],
                dict["see"], dict["maxCount"], dict["loopTime"]]
  params = map(Int, params)
  return Config(params...)
end

max_stars_repo_name: Condimentum/spacetyckiting

max_issues_repo_name: Condimentum/spacetyckiting

max_forks_repo_name: Condimentum/spacetyckiting


text:
function GenericTrace(x::AbstractArray, y::AbstractArray;
                      kind="scatter", kwargs...)
    GenericTrace(kind; x=x, y=y, kwargs... )
end

"""
$(SIGNATURES)
Build a trace of kind `kind`, using the columns of `df` where possible. In
particular for all keyword arguments, if the value of the keyword argument is a
Symbol and matches one of the column names of `df`, replace the value of the
keyword argument with the column of `df`
"""
function GenericTrace(df::AbstractDataFrame; kind="scatter", kwargs...)
    d = Dict{Symbol,Any}(kwargs)

    for (k, v) in d
        if isa(v, Symbol) && haskey(df, v)
            d[k] = df[v]
        end
    end
    GenericTrace(kind; d...)
end

"""
$(SIGNATURES)
Pass the provided values of `x` and `y` as keyword arguments for constructing
the trace from `df`. See other method for more information
"""
function GenericTrace(df::AbstractDataFrame, x::Symbol, y::Symbol; kwargs...)
    GenericTrace(df; x=x, y=y, kwargs...)
end

"""
$(SIGNATURES)
Pass the provided value `y` as keyword argument for constructing the trace from
`df`. See other method for more information
"""
function GenericTrace(df::AbstractDataFrame, y::Symbol; kwargs...)
    GenericTrace(df; y=y, kwargs...)
end

"""
$(SIGNATURES)
Build a plot of with one trace of type `kind`and set `x` to x and `y` to y. All
keyword arguments are passed directly as keyword arguments to the constructed
trace.

**NOTE**: If `y` is a matrix, one trace is constructed for each column of `y`

**NOTE**: If `x` and `y` are both matrices, they must have the same number of
columns (say `N`). Then `N` traces are constructed, where the `i`th column of
`x` is paired with the `i`th column of `y`.
"""
function Plot{T<:_Scalar}(x::AbstractVector{T}, y::AbstractVector, l::Layout=Layout();
              kind="scatter", style::Style=DEFAULT_STYLE[1], kwargs...)
    Plot(GenericTrace(x, y; kind=kind, kwargs...), l, style=style)
end

function Plot{T<:_Scalar}(x::AbstractVector{T}, y::AbstractMatrix, l::Layout=Layout();
              style::Style=DEFAULT_STYLE[1], kwargs...)
    traces = GenericTrace[GenericTrace(x, view(y, :, i); kwargs...)
                          for i in 1:size(y, 2)]
    Plot(traces, l, style=style)
end

function Plot{T<:AbstractVector}(x::AbstractVector{T}, y::AbstractMatrix, l::Layout=Layout();
              style::Style=DEFAULT_STYLE[1], kwargs...)
    size(x, 1) == size(y, 2) || error("x and y must have same number of cols")

    traces = GenericTrace[GenericTrace(x[i], view(y, :, i); kwargs...)
                          for i in 1:size(y,2)]
    Plot(traces, l; style=style)
end

function Plot(x::AbstractMatrix, y::AbstractMatrix, l::Layout=Layout();
              style::Style=DEFAULT_STYLE[1], kwargs...)
    if size(x, 2) == 1
        # use method above
        Plot(view(x, :, 1), y, l; style=style, kwargs...)
    end

    size(x, 2) == size(y, 2) || error("x and y must have same number of cols")

    traces = GenericTrace[GenericTrace(view(x, :, i), view(y, :, i); kwargs...)
                          for i in 1:size(y,2)]
    Plot(traces, l; style=style)
end

"""
$(SIGNATURES)
Build a scatter plot and set  `y` to y. All keyword arguments are passed directly
as keyword arguments to the constructed scatter.
"""
# AbstractArray{T,N}
function Plot{T<:_Scalar}(y::AbstractArray{T}, l::Layout=Layout(); kwargs...)
    # call methods above to get many traces if y is >1d
    Plot(1:size(y, 1), y, l; kwargs...)
end

"""
$(SIGNATURES)
Construct a plot of `f` from `x0` to `x1`, using the layout `l`. All
keyword arguments are applied to the constructed trace.
"""
function Plot(f::Function, x0::Number, x1::Number, l::Layout=Layout();
              style::Style=DEFAULT_STYLE[1],
              kwargs...)
    x = linspace(x0, x1, 50)
    y = [f(_) for _ in x]
    Plot(GenericTrace(x, y; name=Symbol(f), kwargs...), l, style=style)
end

"""
$(SIGNATURES)
For each function in `f` in `fs`, construct a scatter trace that plots `f` from
`x0` to `x1`, using the layout `l`. All keyword arguments are applied to all
constructed traces.
"""
function Plot(fs::AbstractVector{Function}, x0::Number, x1::Number,
              l::Layout=Layout();
              style::Style=DEFAULT_STYLE[1],
              kwargs...)
    x = linspace(x0, x1, 50)
    traces = GenericTrace[GenericTrace(x, map(f, x); name=Symbol(f), kwargs...)
                          for f in fs]
    Plot(traces, l; style=style)
end

"""
$(SIGNATURES)
Construct a plot using the columns of `df` if possible. For each keyword
argument, if the value of the argument is a Symbol and the `df` has a column
whose name matches the value, replace the value with the column of the `df`.

If `group` is passed and is a Symbol that is one of the column names of `df`,
then call `by(df, group)` and construct one trace per SubDataFrame, passing
all other keyword arguments. This means all keyword arguments are passed
applied to all traces
"""
function Plot(df::AbstractDataFrame, l::Layout=Layout(); group=nothing,
              style::Style=DEFAULT_STYLE[1], kwargs...)
    if group != nothing
        # the user passed a group argument, we actually have to do something...
        if isa(group, Symbol) && haskey(df, group)
            _traces = by(df, group) do dfg
                GenericTrace(dfg; name=dfg[1, group], kwargs...)
            end
            traces = GenericTrace[t for t in _traces[:x1]]
            # for some reason I need to bypass calling `plot` due to a type
            # inference error??
            return Plot(traces, l, style=style)
        else
            warn("Unknown group $(group), skipping")
        end
    end
    Plot(GenericTrace(df; kwargs...), l, style=style)
end

"""
$(SIGNATURES)
Construct a plot from `df`, passing the provided values of x and y as keyword
arguments. See docstring for other method for more information.
"""
function Plot(d::AbstractDataFrame, x::Symbol, y::Symbol, l::Layout=Layout();
              style::Style=DEFAULT_STYLE[1], kwargs...)
    Plot(d, l; x=x, y=y, style=style, kwargs...)
end

"""
$(SIGNATURES)
Construct a plot from `df`, passing the provided value y as a keyword argument.
See docstring for other method for more information.
"""
function Plot(d::AbstractDataFrame, y::Symbol, l::Layout=Layout();
              style::Style=DEFAULT_STYLE[1], kwargs...)
    Plot(d, l; y=y, style=style, kwargs...)
end

"""
$(SIGNATURES)
Creates a "stem" or "lollipop" trace. It is implemented using plotly.js's
`scatter` type, using the error bars to draw the stem.

## Keyword Arguments:
* All properties accepted by `scatter` except `error_y`, which is used to draw
    the stems
* stem_color - sets the color of the stems
* stem_thickness - sets the thickness of the stems
"""
function stem(;y=nothing, stem_color="grey", stem_thickness=1, kwargs...)
    line_up = -min(y, 0)
    line_down = max(y, 0)
    trace = scatter(; y=y, text=y, marker_size=10, mode="markers", hoverinfo="text", kwargs...)
    trace.fields[:error_y] = Dict(
        :type => "data",
        :symmetric => false,
        :array => line_up,
        :arrayminus => line_down,
        :visible => true,
        :color => stem_color,
        :width => 0,
        :thickness => stem_thickness)
    trace
end


text:
struct RenderInfo
    shader_dependencies::ShaderDependencies
    pipeline::Ref{Created{Pipeline,GraphicsPipelineCreateInfo}}
    draw_args::NTuple{4,Int}
    push_data::Any
end

function RenderInfo(
        deps::ShaderDependencies,
        draw_args;
        pipeline = Ref{Created{Pipeline,GraphicsPipelineCreateInfo}}(),
        push_data = nothing,
    )
    RenderInfo(deps, draw_args, pipeline, push_data)
end

struct RenderSet
    pipeline_cache::GraphicsPipelineCache
    render_infos::Dictionary{Any,RenderInfo}
end

RenderSet(device) = RenderSet(GraphicsPipelineCache(device))
RenderSet(cache::GraphicsPipelineCache) = RenderSet(cache, Dictionary{Any,RenderInfo}())

@forward RenderSet.render_infos Base.insert!, Dictionaries.set!, Base.setindex!

function prepare!(f, set::RenderSet)
    batch = Dictionary{GraphicsPipelineCreateInfo,Ref{Pipeline}}()
    for (object, info) in pairs(set.render_infos)
        prepare_pipeline!(f, batch, set, object, info)
        prepare_resources!(set, object, info)
    end
    if !isempty(batch)
        pipelines = get_graphics_pipelines!(set.pipeline_cache, collect(keys(batch)))
        for ((info, ref), pipeline) in zip(pairs(batch), pipelines)
            ref[] = Created(pipeline, info)
        end
    end
end

function prepare_pipeline!(f, batch, set::RenderSet, object, info)
    if !isdefined(info.pipeline, 1)
        insert!(batch, info.pipeline, f(object))
    end
end

function BindRequirements(info::RenderInfo)
    pipeline_info = render_info.pipeline.info
    push_ranges = pipeline_info.push_constant_ranges == C_NULL ? nothing : pipeline_info.push_constant_ranges
    BindRequirements(info.shader_dependencies, PipelineBindingState(pipeline, push_ranges), pipeline_info.pipeline_layout, info.push_data)
end


text:
using Images, Random
import Base.Threads: @threads, @spawn

function escapetime(z; maxiter=80)
    c = z
    for n = 1:maxiter
        if abs(z) > 2
            return n-1
        end
        z = z^2 + c
    end
    return maxiter
end

function mandel(; width=80, height=20, maxiter=80)
    out = zeros(Int, height, width)
    real = range(-2.0, 0.5, length=width)
    imag = range(-1.0, 1.0, length=height)
    for x in 1:width
        for y in 1:height
            z = real[x] + imag[y]*im
            out[y,x] = escapetime(z, maxiter=maxiter)
        end
    end
    return out
end

@time m = mandel(width=1200,height=900,maxiter=400)
Gray.(m./80)

## CSV parsing

using CSV, BenchmarkTools

run(`python3 -m timeit -s "import pandas" -p "pandas.read_csv('mixed.csv')"`)

@btime CSV.read("mixed.csv", threaded=false)

@btime CSV.read("mixed.csv", threaded=true)


text:
struct QuickMDP{ID,S,A,D<:NamedTuple} <: MDP{S,A}
    data::D
end

"""
    QuickMDP(gen::Function, [id]; kwargs...)

Construct a generative MDP model with the function `gen` and keyword arguments.

`gen` should take three arguments: a state, an action, and a random number generator. It should return a `NamedTuple` with keys `sp` for the next state and `r` for the reward.

Keywords can be static objects or functions. See the QuickPOMDPs.jl documentation for more information.
"""
QuickMDP(gen::Function, id=uuid4(); kwargs...) = QuickMDP(id; gen=gen, kwargs...)

"""
    QuickMDP([id]; kwargs...)

Construct an MDP model with keyword arguments. Keywords can be static objects or functions. See the QuickPOMDPs.jl documentation for more information.
"""
function QuickMDP(id=uuid4(); kwargs...)
    kwd = Dict{Symbol, Any}(kwargs)

    quick_defaults!(kwd)

    S = infer_statetype(kwd)
    A = infer_actiontype(kwd)
    d = namedtuple(keys(kwd)...)(values(kwd)...)
    qm = QuickMDP{id, S, A, typeof(d)}(d)
    return qm
end

id(::QuickMDP{ID}) where ID = ID

struct QuickPOMDP{ID,S,A,O,D<:NamedTuple} <: POMDP{S,A,O}
    data::D
end

"""
    QuickPOMDP(gen::Function, [id]; kwargs...)

Construct a generative POMDP model with the function `gen` and keyword arguments.

`gen` should take three arguments: a state, an action, and a random number generator. It should return a `NamedTuple` with keys `sp` for the next state, `o` for the observation, and `r` for the reward.

Keywords can be static objects or functions. See the QuickPOMDPs.jl documentation for more information.
"""
QuickPOMDP(gen::Function, id=uuid4(); kwargs...) = QuickPOMDP(id; gen=gen, kwargs...)

"""
    QuickPOMDP([id]; kwargs...)

Construct an POMDP model with keyword arguments. Keywords can be static objects or functions. See the QuickPOMDPs.jl documentation for more information.
"""
function QuickPOMDP(id=uuid4(); kwargs...)
    kwd = Dict{Symbol, Any}(kwargs)

    quick_defaults!(kwd)

    S = infer_statetype(kwd)
    A = infer_actiontype(kwd)
    O = infer_obstype(kwd)
    d = namedtuple(keys(kwd)...)(values(kwd)...)
    qm = QuickPOMDP{id, S, A, O, typeof(d)}(d)
    return qm
end

id(::QuickPOMDP{ID}) where ID = ID

const QuickModel = Union{QuickMDP, QuickPOMDP}

function quick_defaults!(kwd::Dict)
    kwd[:discount] = get(kwd, :discount, 1.0)
    kwd[:isterminal] = get(kwd, :isterminal, false)

    # memoize initialstate_distribution since it should be constant (so we can use it below for initialstate)
    if haskey(kwd, :initialstate_distribution) && kwd[:initialstate_distribution] isa Function
        kwd[:initialstate_distribution] = kwd[:initialstate_distribution]()
    end

    if !haskey(kwd, :initialstate)
        if haskey(kwd, :initialstate_distribution)
            kwd[:initialstate] = rng -> rand(rng, kwd[:initialstate_distribution])
        end
    end

    # default for initialobs must be in the method below because the method table might change

    if !haskey(kwd, :stateindex)
        if haskey(kwd, :states)
            states = _call(Val(:states), kwd[:states], ())
            if hasmethod(length, typeof((states,))) && length(states) < Inf
                kwd[:stateindex] = Dict(s=>i for (i,s) in enumerate(states))
            end
        end
    end

    if !haskey(kwd, :actionindex)
        if haskey(kwd, :actions)
            actions = _call(Val(:actions), kwd[:actions], ())
            if hasmethod(length, typeof((actions,))) && length(actions) < Inf
                kwd[:actionindex] = Dict(s=>i for (i,s) in enumerate(actions))
            end
        end
    end

    if !haskey(kwd, :obsindex)
        if haskey(kwd, :observations)
            observations = _call(Val(:observations), kwd[:observations], ())
            if hasmethod(length, typeof((observations,))) && length(observations) < Inf
                kwd[:obsindex] = Dict(s=>i for (i,s) in enumerate(observations))
            end
        end
    end
end

function infer_statetype(kwd)
    if haskey(kwd, :statetype)
        st = _call(Val(:statetype), kwd[:statetype], (), NamedTuple())
    elseif haskey(kwd, :states)
        st = eltype(_call(Val(:states), kwd[:states], (), NamedTuple()))
    elseif haskey(kwd, :initialstate)
        st = typeof(_call(Val(:initialstate), kwd[:initialstate], (MersenneTwister(0),), NamedTuple()))
    else
        st = Any
    end
    if st == Any
        @warn("Unable to infer state type for a Quick(PO)MDP; using Any. This may have significant performance consequences. Use the statetype keyword argument to specify a concrete state type.")
    end
    return st
end

function infer_actiontype(kwd)
    if haskey(kwd, :actiontype)
        at = _call(Val(:actiontype), kwd[:actiontype], (), NamedTuple())
    elseif haskey(kwd, :actions)
        at = eltype(_call(Val(:actions), kwd[:actions], (), NamedTuple()))
    else
        at = Any
    end
    if at == Any
        @warn("Unable to infer action type for a Quick(PO)MDP; using Any. This may have significant performance consequences. Use the actiontype keyword argument to specify a concrete action type.")
    end
    return at
end

function infer_obstype(kwd)
    if haskey(kwd, :obstype)
        ot = _call(Val(:obstype), kwd[:obstype], (), NamedTuple())
    elseif haskey(kwd, :observations)
        ot = eltype(_call(Val(:observations), kwd[:observations], (), NamedTuple()))
    elseif haskey(kwd, :initialobs) && haskey(kwd, :initialstate)
        s0 = _call(Val(:initialstate), kwd[:initialstate], (MersenneTwister(0),), NamedTuple())
        ot = typeof(_call(Val(:initialobs), kwd[:initialobs], (s0, MersenneTwister(0),), NamedTuple()))
    else
        ot = Any
    end
    if ot == Any
        @warn("Unable to infer observation type for a QuickPOMDP; using Any. This may have significant performance consequences. Use the obstype keyword argument to specify a concrete observation type.")
    end
    return ot
end

function _call(namev::Val{name}, m::QuickModel, args, kwargs=NamedTuple()) where name
    _call(namev,
          get(m.data, name) do
              throw(MissingQuickArgument(m, name))
          end,
          args,
          kwargs)
end

_call(::Val, f::Function, args, kwargs=NamedTuple()) = f(args...; kwargs...)
_call(v::Val, object, args, kwargs=NamedTuple()) = object
_call(v::Val, d::Dict, args, kwargs=NamedTuple()) = d[args]

macro forward_to_data(f)
    quote
        $f(m::QuickModel, args...; kwargs...) = _call(Val(nameof($f)), m, args, kwargs)
    end
end

function POMDPs.transition(m::QuickModel, s, a)
    if haskey(m.data, :transition)
        return m.data.transition(s, a)
    else
        throw(MissingQuickArgument(m, :transition, types=[Function], also=[:gen]))
    end
end

function POMDPs.observation(m::QuickPOMDP, args...)
    if haskey(m.data, :observation)
        return m.data.observation(args...)
    else
        throw(MissingQuickArgument(m, :observation, types=[Function], also=[:gen]))
    end
end

@forward_to_data POMDPs.initialstate_distribution
@forward_to_data POMDPs.reward

function POMDPs.gen(m::QuickModel, s, a, rng)
    if haskey(m.data, :gen)
        return m.data.gen(s, a, rng)
    else
        return NamedTuple()
    end
end

POMDPs.initialstate(m::QuickModel, rng::AbstractRNG) = _call(Val(:initialstate), m, (rng,))

function POMDPs.initialobs(m::QuickPOMDP, s, rng::AbstractRNG)
    if haskey(m.data, :initialobs)
        return _call(Val(:initialobs), m, (s, rng))
    elseif haskey(m.data, :observation) && hasmethod(m.data.observation, typeof((s,)))
        return rand(rng, m.data.observation(s))
    else
        throw(MissingQuickArgument(m, :initialobs, types=[obstype(m), Function], also=[:observation]))
    end
end

@forward_to_data POMDPs.states
@forward_to_data POMDPs.actions
@forward_to_data POMDPs.observations
@forward_to_data POMDPs.discount

@forward_to_data POMDPs.stateindex
@forward_to_data POMDPs.actionindex
@forward_to_data POMDPs.obsindex
@forward_to_data POMDPs.isterminal

@forward_to_data POMDPModelTools.obs_weight
@forward_to_data POMDPModelTools.render


text:
"""
    BetaBinomial(n,α,β)

A *Beta-binomial distribution* is the compound distribution of the [`Binomial`](@ref) distribution where the probability of success `p` is distributed according to the [`Beta`](@ref). It has three parameters: `n`, the number of trials and two shape parameters `α`, `β`

```math
P(X = k) = {n \\choose k} B(k + \\alpha, n - k + \\beta) / B(\\alpha, \\beta),  \\quad \\text{ for } k = 0,1,2, \\ldots, n.
```

```julia
BetaBinomial(n, α, β)      # BetaBinomial distribution with n trials and shape parameters α, β

params(d)       # Get the parameters, i.e. (n, α, β)
ntrials(d)      # Get the number of trials, i.e. n
```

External links:

* [Beta-binomial distribution on Wikipedia](https://en.wikipedia.org/wiki/Beta-binomial_distribution)
"""
struct BetaBinomial{T<:Real} <: DiscreteUnivariateDistribution
    n::Int
    α::T
    β::T

    BetaBinomial{T}(n::Integer, α::T, β::T) where {T <: Real} = new{T}(n, α, β)
end

function BetaBinomial(n::Integer, α::T, β::T; check_args=true) where {T <: Real}
    check_args && @check_args(BetaBinomial, n >= zero(n) && α >= zero(α) && β >= zero(β))
    return BetaBinomial{T}(n, α, β)
end

BetaBinomial(n::Integer, α::Real, β::Real) = BetaBinomial(n, promote(α, β)...)
BetaBinomial(n::Integer, α::Integer, β::Integer) = BetaBinomial(n, float(α), float(β))

@distr_support BetaBinomial 0 d.n

#### Conversions
function convert(::Type{BetaBinomial{T}}, n::Int, α::S, β::S) where {T <: Real, S <: Real}
    BetaBinomial(n, T(α), T(β))
end
function convert(::Type{BetaBinomial{T}}, d::BetaBinomial{S}) where {T <: Real, S <: Real}
    BetaBinomial(d.n, T(d.α), T(d.β), check_args=false)
end

#### Parameters

ntrials(d::BetaBinomial) = d.n

params(d::BetaBinomial) = (d.n, d.α, d.β)
partype(::BetaBinomial{T}) where {T} = T

#### Properties

mean(d::BetaBinomial) = (d.n * d.α) / (d.α + d.β)
function var(d::BetaBinomial)
    n, α, β = d.n, d.α, d.β
    numerator = n * α * β * (α + β + n)
    denominator = (α + β)^2 * (α + β + 1)
    return numerator / denominator
end

function skewness(d::BetaBinomial)
    n, α, β = d.n, d.α, d.β
    t1 = (α + β + 2n) * (β - α) / (α + β + 2)
    t2 = sqrt((1 + α +β) / (n * α * β * (n + α + β)))
    return t1 * t2
end

function kurtosis(d::BetaBinomial)
    n, α, β = d.n, d.α, d.β
    alpha_beta_sum = α + β
    alpha_beta_product = α * β
    numerator = ((alpha_beta_sum)^2) * (1 + alpha_beta_sum)
    denominator = (n * alpha_beta_product) * (alpha_beta_sum + 2) * (alpha_beta_sum + 3) * (alpha_beta_sum + n)
    left = numerator / denominator
    right = (alpha_beta_sum) * (alpha_beta_sum - 1 + 6n) + 3*alpha_beta_product * (n - 2) + 6n^2
    right -= (3*alpha_beta_product * n * (6 - n)) / alpha_beta_sum
    right -= (18*alpha_beta_product * n^2) / (alpha_beta_sum)^2
    return (left * right) - 3
end

function logpdf(d::BetaBinomial, k::Real)
    n, α, β = d.n, d.α, d.β
    _insupport = insupport(d, k)
    _k = _insupport ? round(Int, k) : 0
    logbinom = - log1p(n) - logbeta(_k + 1, n - _k + 1)
    lognum   = logbeta(_k + α, n - _k + β)
    logdenom = logbeta(α, β)
    result = logbinom + lognum - logdenom
    return _insupport ? result : oftype(result, -Inf)
end

entropy(d::BetaBinomial) = entropy(Categorical(pdf.(Ref(d),support(d))))
median(d::BetaBinomial) = median(Categorical(pdf.(Ref(d),support(d)))) - 1
mode(d::BetaBinomial) = argmax(pdf.(Ref(d),support(d))) - 1
modes(d::BetaBinomial) = modes(Categorical(pdf.(Ref(d),support(d)))) .- 1

quantile(d::BetaBinomial, p::Float64) = quantile(Categorical(pdf.(Ref(d), support(d))), p) - 1

#### Sampling

rand(rng::AbstractRNG, d::BetaBinomial) =
    rand(rng, Binomial(d.n, rand(rng, Beta(d.α, d.β))))


text:
using Raylib

function main()
    screenWidth = 800
    screenHeight = 450

    Raylib.InitWindow(screenWidth, screenHeight, "raylib [core] example - drop files")

    Raylib.SetTargetFPS(60)

    droppedFiles = Vector{String}()

    while !Raylib.WindowShouldClose()

        if Raylib.IsFileDropped()
            droppedFiles = Raylib.GetDroppedFiles()
        end

        Raylib.BeginDrawing()
        Raylib.ClearBackground(Raylib.RAYWHITE)
        if iszero(length(droppedFiles))
            Raylib.DrawText("Drop your files to this window!", 100, 40, 20, Raylib.DARKGRAY)
        else
            Raylib.DrawText("Dropped files:", 100, 40, 20, Raylib.DARKGRAY)

            for i = 1:length(droppedFiles)
                if i % 2 == 1
                    Raylib.DrawRectangle(0, 85 + 40*(i-1), screenWidth, 40, Raylib.Fade(Raylib.LIGHTGRAY, 0.5))
                else
                    Raylib.DrawRectangle(0, 85 + 40*(i-1), screenWidth, 40, Raylib.Fade(Raylib.LIGHTGRAY, 0.3))
                end

                Raylib.DrawText(droppedFiles[i], 120, 100 + 40*(i-1), 10, Raylib.GRAY)

            end
            Raylib.DrawText("Dropped new files...", 100, 110 + 40length(droppedFiles), 20, Raylib.DARKGRAY)
        end

        Raylib.EndDrawing()
    end

    Raylib.ClearDroppedFiles()
    Raylib.CloseWindow()
end


text:
# Autogenerated wrapper script for TauDEM_jll for x86_64-apple-darwin
export aread8, areadinf, d8flowdir, d8flowpathextremeup, d8hdisttostrm, dinfavalanche, dinfconclimaccum, dinfdecayaccum, dinfdistdown, dinfdistup, dinfflowdir, dinfrevaccum, dinftranslimaccum, dinfupdependence, dropanalysis, gagewatershed, gridnet, lengtharea, moveoutletstostrm, peukerdouglas, pitremove, slopearea, slopearearatio, slopeavedown, streamnet, threshold

using GDAL_jll
using MPICH_jll
using MicrosoftMPI_jll
JLLWrappers.@generate_wrapper_header("TauDEM")
JLLWrappers.@declare_executable_product(aread8)
JLLWrappers.@declare_executable_product(areadinf)
JLLWrappers.@declare_executable_product(d8flowdir)
JLLWrappers.@declare_executable_product(d8flowpathextremeup)
JLLWrappers.@declare_executable_product(d8hdisttostrm)
JLLWrappers.@declare_executable_product(dinfavalanche)
JLLWrappers.@declare_executable_product(dinfconclimaccum)
JLLWrappers.@declare_executable_product(dinfdecayaccum)
JLLWrappers.@declare_executable_product(dinfdistdown)
JLLWrappers.@declare_executable_product(dinfdistup)
JLLWrappers.@declare_executable_product(dinfflowdir)
JLLWrappers.@declare_executable_product(dinfrevaccum)
JLLWrappers.@declare_executable_product(dinftranslimaccum)
JLLWrappers.@declare_executable_product(dinfupdependence)
JLLWrappers.@declare_executable_product(dropanalysis)
JLLWrappers.@declare_executable_product(gagewatershed)
JLLWrappers.@declare_executable_product(gridnet)
JLLWrappers.@declare_executable_product(lengtharea)
JLLWrappers.@declare_executable_product(moveoutletstostrm)
JLLWrappers.@declare_executable_product(peukerdouglas)
JLLWrappers.@declare_executable_product(pitremove)
JLLWrappers.@declare_executable_product(slopearea)
JLLWrappers.@declare_executable_product(slopearearatio)
JLLWrappers.@declare_executable_product(slopeavedown)
JLLWrappers.@declare_executable_product(streamnet)
JLLWrappers.@declare_executable_product(threshold)
function __init__()
    JLLWrappers.@generate_init_header(GDAL_jll, MPICH_jll, MicrosoftMPI_jll)
    JLLWrappers.@init_executable_product(
        aread8,
        "bin/aread8",
    )

    JLLWrappers.@init_executable_product(
        areadinf,
        "bin/areadinf",
    )

    JLLWrappers.@init_executable_product(
        d8flowdir,
        "bin/d8flowdir",
    )

    JLLWrappers.@init_executable_product(
        d8flowpathextremeup,
        "bin/d8flowpathextremeup",
    )

    JLLWrappers.@init_executable_product(
        d8hdisttostrm,
        "bin/d8hdisttostrm",
    )

    JLLWrappers.@init_executable_product(
        dinfavalanche,
        "bin/dinfavalanche",
    )

    JLLWrappers.@init_executable_product(
        dinfconclimaccum,
        "bin/dinfconclimaccum",
    )

    JLLWrappers.@init_executable_product(
        dinfdecayaccum,
        "bin/dinfdecayaccum",
    )

    JLLWrappers.@init_executable_product(
        dinfdistdown,
        "bin/dinfdistdown",
    )

    JLLWrappers.@init_executable_product(
        dinfdistup,
        "bin/dinfdistup",
    )

    JLLWrappers.@init_executable_product(
        dinfflowdir,
        "bin/dinfflowdir",
    )

    JLLWrappers.@init_executable_product(
        dinfrevaccum,
        "bin/dinfrevaccum",
    )

    JLLWrappers.@init_executable_product(
        dinftranslimaccum,
        "bin/dinftranslimaccum",
    )

    JLLWrappers.@init_executable_product(
        dinfupdependence,
        "bin/dinfupdependence",
    )

    JLLWrappers.@init_executable_product(
        dropanalysis,
        "bin/dropanalysis",
    )

    JLLWrappers.@init_executable_product(
        gagewatershed,
        "bin/gagewatershed",
    )

    JLLWrappers.@init_executable_product(
        gridnet,
        "bin/gridnet",
    )

    JLLWrappers.@init_executable_product(
        lengtharea,
        "bin/lengtharea",
    )

    JLLWrappers.@init_executable_product(
        moveoutletstostrm,
        "bin/moveoutletstostrm",
    )

    JLLWrappers.@init_executable_product(
        peukerdouglas,
        "bin/peukerdouglas",
    )

    JLLWrappers.@init_executable_product(
        pitremove,
        "bin/pitremove",
    )

    JLLWrappers.@init_executable_product(
        slopearea,
        "bin/slopearea",
    )

    JLLWrappers.@init_executable_product(
        slopearearatio,
        "bin/slopearearatio",
    )

    JLLWrappers.@init_executable_product(
        slopeavedown,
        "bin/slopeavedown",
    )

    JLLWrappers.@init_executable_product(
        streamnet,
        "bin/streamnet",
    )

    JLLWrappers.@init_executable_product(
        threshold,
        "bin/threshold",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the conversion from Euler angles to DCM.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export angle_to_dcm

"""
    angle_to_dcm(θ₁::Number, θ₂::Number, θ₃::Number, rot_seq::Symbol = :ZYX)
    angle_to_dcm(Θ::EulerAngles)

Convert the Euler angles `θ₁`, `θ₂`, and `θ₃` [rad] with the rotation sequence
`rot_seq` to a direction cosine matrix.

The input values of the origin Euler angles can also be passed inside the
structure `Θ` (see [`EulerAngles`](@ref)).

The rotation sequence is defined by a `:Symbol`. The possible values are:
`:XYX`, `XYZ`, `:XZX`, `:XZY`, `:YXY`, `:YXZ`, `:YZX`, `:YZY`, `:ZXY`, `:ZXZ`,
`:ZYX`, and `:ZYZ`. If no value is specified, then it defaults to `:ZYX`.

# Remarks

This function assigns `dcm = A3 * A2 * A1` in which `Ai` is the DCM related with
the *i*-th rotation, `i Є [1,2,3]`.

# Example

```jldoctest
julia> dcm = angle_to_dcm(pi / 2, pi / 3, pi / 4, :ZYX)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  3.06162e-17  0.5       -0.866025
 -0.707107     0.612372   0.353553
  0.707107     0.612372   0.353553
```
"""
function angle_to_dcm(θ₁::Number, θ₂::Number, θ₃::Number, rot_seq::Symbol = :ZYX)
    # Compute the sines and cosines.
    s₁, c₁ = sincos(θ₁)
    s₂, c₂ = sincos(θ₂)
    s₃, c₃ = sincos(θ₃)

    if rot_seq == :ZYX
        return DCM(
                 c₂ * c₁,                c₂ * s₁,             -s₂ ,
            s₃ * s₂ * c₁ - c₃ * s₁, s₃ * s₂ * s₁ + c₃ * c₁, s₃ * c₂,
            c₃ * s₂ * c₁ + s₃ * s₁, c₃ * s₂ * s₁ - s₃ * c₁, c₃ * c₂
        )'
    elseif rot_seq == :XYX
        return DCM(
              c₂,               s₁ * s₂,               -c₁ * s₂,
            s₂ * s₃, -s₁ * c₂ * s₃ + c₁ * c₃, c₁ * c₂ * s₃ + s₁ * c₃,
            s₂ * c₃, -s₁ * c₃ * c₂ - c₁ * s₃, c₁ * c₃ * c₂ - s₁ * s₃
        )'
    elseif rot_seq == :XYZ
        return DCM(
             c₂ * c₃,  s₁ * s₂ * c₃ + c₁ * s₃, -c₁ * s₂ * c₃ + s₁ * s₃,
            -c₂ * s₃, -s₁ * s₂ * s₃ + c₁ * c₃,  c₁ * s₂ * s₃ + s₁ * c₃,
                s₂,             -s₁ * c₂,                 c₁ * c₂
        )'
    elseif rot_seq == :XZX
        return DCM(
               c₂,               c₁ * s₂,                 s₁ * s₂,
            -s₂ * c₃,  c₁ * c₃ * c₂ - s₁ * s₃,  s₁ * c₃ * c₂ + c₁ * s₃,
             s₂ * s₃, -c₁ * c₂ * s₃ - s₁ * c₃, -s₁ * c₂ * s₃ + c₁ * c₃
        )'

    elseif rot_seq == :XZY
        return DCM(
            c₃ * c₂, c₁ * c₃ * s₂ + s₁ * s₃, s₁ * c₃ * s₂ - c₁ * s₃,
              -s₂,             c₁ * c₂,                s₁ * c₂,
            s₃ * c₂, c₁ * s₂ * s₃ - s₁ * c₃, s₁ * s₂ * s₃ + c₁ * c₃
        )'
    elseif rot_seq == :YXY
        return DCM(
            -s₁ * c₂ * s₃ + c₁ * c₃,  s₂ * s₃, -c₁ * c₂ * s₃ - s₁ * c₃,
                       s₁ * s₂,         c₂,              c₁ * s₂,
             s₁ * c₃ * c₂ + c₁ * s₃, -s₂ * c₃,  c₁ * c₃ * c₂ - s₁ * s₃
        )'
    elseif rot_seq == :YXZ
        return DCM(
             c₁ * c₃ + s₂ * s₁ * s₃, c₂ * s₃, -s₁ * c₃ + s₂ * c₁ * s₃,
            -c₁ * s₃ + s₂ * s₁ * c₃, c₂ * c₃,  s₁ * s₃ + s₂ * c₁ * c₃,
                  s₁ * c₂,             -s₂,         c₂ * c₁
        )'
    elseif rot_seq == :YZX
        return DCM(
                       c₁ * c₂,         s₂,              -s₁ * c₂,
            -c₃ * c₁ * s₂ + s₃ * s₁,  c₂ * c₃,  c₃ * s₁ * s₂ + s₃ * c₁,
             s₃ * c₁ * s₂ + c₃ * s₁, -s₃ * c₂, -s₃ * s₁ * s₂ + c₃ * c₁
        )'
    elseif rot_seq == :YZY
        return DCM(
            c₁ * c₃ * c₂ - s₁ * s₃, s₂ * c₃, -s₁ * c₃ * c₂ - c₁ * s₃,
                -c₁ * s₂,             c₂,               s₁ * s₂,
            c₁ * c₂ * s₃ + s₁ * c₃, s₂ * s₃, -s₁ * c₂ * s₃ + c₁ * c₃
        )'
    elseif rot_seq == :ZXY
        return DCM(
            c₃ * c₁ - s₂ * s₃ * s₁, c₃ * s₁ + s₂ * s₃ * c₁, -s₃ * c₂,
                -c₂ * s₁,                c₂ * c₁,              s₂,
            s₃ * c₁ + s₂ * c₃ * s₁, s₃ * s₁ - s₂ * c₃ * c₁,  c₂ * c₃
        )'
    elseif rot_seq == :ZXZ
        return DCM(
            -s₁ * c₂ * s₃ + c₁ * c₃, c₁ * c₂ * s₃ + s₁ * c₃, s₂ * s₃,
            -s₁ * c₃ * c₂ - c₁ * s₃, c₁ * c₃ * c₂ - s₁ * s₃, s₂ * c₃,
                       s₁ * s₂,               -c₁ * s₂,         c₂
        )'
    elseif rot_seq == :ZYZ
        return DCM(
             c₁ * c₃ * c₂ - s₁ * s₃,  s₁ * c₃ * c₂ + c₁ * s₃, -s₂ * c₃,
            -c₁ * c₂ * s₃ - s₁ * c₃, -s₁ * c₂ * s₃ + c₁ * c₃,  s₂ * s₃,
                       c₁ * s₂,                 s₁ * s₂,          c₂
        )'
    else
        throw(ArgumentError("The rotation sequence :$rot_seq is not valid."))
    end
end

angle_to_dcm(Θ::EulerAngles) = angle_to_dcm(Θ.a1, Θ.a2, Θ.a3, Θ.rot_seq)


text:
using Distributed, K8sClusterManagers
using DataFrames
using CSV

workers = parse(Int,ARGS[1])
cpu_fraction = parse(Int,ARGS[2])
addprocs(K8sClusterManager(workers,cpu=cpu_fraction); exeflags="--project")

#addprocs(K8sClusterManager(10, cpu=1, memory="300Mi", pending_timeout=300); exeflags="--project")
#nprocs() == 1 && addprocs(; exeflags = "--project")

@everywhere include("twoblocks.jl")
@everywhere using Main.TwoBlocksPipeline

dataset = getiris()
X = dataset[:,1:4]
Y = dataset[:,5] |> collect
X = dataset[:,2:end-1]

#dataset=CSV.read("environment_data.csv",DataFrame)
#X = select(dataset, Not([:state,:time]))
#Y = dataset.state

results=TwoBlocksPipeline.twoblockspipelinesearch(X,Y)
println(results)
println("best pipeline is: ", results[1,1],", with mean-sd ",results[1,2], " ± ",results[1,3])


text:
function make_initialize_paremeters(ℼ::Float64 = 1.02)

    function initialize_paremeters()

        eval("Inflation Target")
        eval(Meta.parse("const global 𝜋 = "*string(ℼ)) )
        # params
        @eval begin

            "Utility function parameter"
            const global σ = 2;
            "Discount Factor"
            const global β = 0.96;

            "Share of debt to be repaid each period"
            const global δᶠ = 0.2;
            "Share of share of debt that must be stored as deposit/liquidity"
            const global δˡ = .8;
            const global pₑᵤ = 0.08;
            const global pᵤₑ = 3/7;
            "Status probability matrix"
            const global λ = [1-pₑᵤ pₑᵤ;pᵤₑ 1-pᵤₑ] ;
            "Labor Supply"
            const global Lₛ = (λ^2000)[1]
            "Probability of reentering market if defaulted"
            const global ζ = 0.5;
            "unemployment benefit (real)"
            const global ϐ = 0.25 ;
            "Chi. Paremeter representing increased cost of consumption if defaulted (from 0 - absolute cost to 1 - no cost and decreased cost if 1+)"
            const global Χ = 1  ;
            "Banks - depositors monopolistic advatage parameter"
            const global ρ = .95 ;
            "Liquidity in advance punishment parameter"
            const global ξ = MathConstants.φ - 1

        end
        # const global AreBanksRetarded = false

    end

    return (initialize_paremeters)

end

initialize_paremeters = make_initialize_paremeters()
initialize_paremeters()


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>100.0, :alg=>TMJets(abs_tol=1e-5, orderT=5, orderQ=2, max_steps=50_000), :disjointness_method=>BoxEnclosure(), :intersect_source_invariant=>false),
    "X0"   => [(1, Hyperrectangle(low=[-65.0, -0.2], high=[-60.0, 0.2]))],
    "info" => Dict(
        "name" => "spiking_neuron",
        "dim" => 2,
        "linear" => false,
        "hybrid" => true
    )
])


text:
"""
Axisymmetric Solov'ev equilibra in (R/R₀,Z/R₀,ϕ) coordinates.
Based on Cerfon & Freidberg, Physics of Plasmas 17, 032502, 2010,
      and Freidberg, Ideal Magnetohydrodynamics, 2014.
"""
module Solovev

    using RecipesBase

    using SymEngine: N, symbols, diff, expand, subs

    import ..ElectromagneticFields
    import ..ElectromagneticFields: code
    import ..SolovevAbstract: AbstractSolovevEquilibrium, X, Y, Z, R, r, θ, ϕ, r²

    export  SolovevEquilibrium, SolovevXpointEquilibrium

    function ψ₀(x::AbstractVector{T}, a) where {T <: Number}
        x[1]^4 / 8 + a * (x[1]^2 * log(x[1]) / 2 - x[1]^4 / 8 )
    end

    function ψ₁(x::AbstractVector{T}) where {T <: Number}
        one(T)
    end

    function ψ₂(x::AbstractVector{T}) where {T <: Number}
        x[1]^2
    end

    function ψ₃(x::AbstractVector{T}) where {T <: Number}
        x[2]^2 - x[1]^2 * log(x[1])
    end

    function ψ₄(x::AbstractVector{T}) where {T <: Number}
        x[1]^4 - 4 * x[1]^2 * x[2]^2
    end

    function ψ₅(x::AbstractVector{T}) where {T <: Number}
        2 * x[2]^4 - 9 * x[2]^2 * x[1]^2 + 3 * x[1]^4 * log(x[1]) - 12 * x[1]^2 * x[2]^2 * log(x[1])
    end

    function ψ₆(x::AbstractVector{T}) where {T <: Number}
        x[1]^6 - 12 * x[1]^4 * x[2]^2 + 8 * x[1]^2 * x[2]^4
    end

    function ψ₇(x::AbstractVector{T}) where {T <: Number}
        8 * x[2]^6 - 140 * x[2]^4 * x[1]^2 + 75 * x[2]^2 * x[1]^4 - 15 * x[1]^6 * log(x[1]) +
                    180 * x[1]^4 * x[2]^2 * log(x[1]) - 120 * x[1]^2 * x[2]^4 * log(x[1])
    end

    function ψ₈(x::AbstractVector{T}) where {T <: Number}
        x[2]
    end

    function ψ₉(x::AbstractVector{T}) where {T <: Number}
        x[2] * x[1]^2
    end

    function ψ₁₀(x::AbstractVector{T}) where {T <: Number}
        x[2]^3 - 3 * x[2] * x[1]^2 * log(x[1])
    end

    function ψ₁₁(x::AbstractVector{T}) where {T <: Number}
        3 * x[2] * x[1]^4 - 4 * x[2]^3 * x[1]^2
    end

    function ψ₁₂(x::AbstractVector{T}) where {T <: Number}
        8 * x[2]^5 - 45 * x[2] * x[1]^4 - 80 * x[2]^3 * x[1]^2 * log(x[1]) + 60 * x[2] * x[1]^4 * log(x[1])
    end

    @doc raw"""
    Axisymmetric Solov'ev equilibra in (R/R₀,Z/R₀,ϕ) coordinates.
    Based on Cerfon & Freidberg, Physics of Plasmas 17, 032502, 2010,
        and Freidberg, Ideal Magnetohydrodynamics, 2014.

    The covariant components of the vector potential are given by
    ```math
    A (x, y, \phi) = \left( \frac{B_0 R_0}{2} \, \frac{y}{x} , \, - \frac{B_0 R_0}{2} \, \ln x , \, \psi(x,y) \right)^T ,
    ```

    with $x = R/R_0$ and $y = Z/R_0$.
    The normalised poloidal flux $\psi$ is given by
    ```math
    \psi (x,y) = \psi_0 + \sum \limits_{i=1}^{7} c_i \psi_i (x,y) ,
    ```

    with
    ```math
    \begin{aligned}
    \psi_{0} &= \frac{x^4}{8} + \alpha \left( \frac{1}{2} x^2 \, \ln x - \frac{x^4}{8} \right) , \\
    \psi_{1} &= 1 , \\
    \psi_{2} &= x^2 , \\
    \psi_{3} &= y^2 - x^2 \, \ln x , \\
    \psi_{4} &= x^4 - 4 x^2 y^2 , \\
    \psi_{5} &= 2 y^4 9 y^2 x^2 + 3 x^4 \, \ln x - 12 x^2 y^2 \, \ln x , \\
    \psi_{6} &= x^6 - 12 x^4 y^2 + 8 x^2 y^4 , \\
    \psi_{7} &= 8 y^6 - 140 y^4 x^2 + 75 y^2 x^4 - 15 x^6 \, \ln x + 180 x^4 y^2 \, \ln x - 120 x^2 y^4 \, \ln x .
    \end{aligned}
    ```

    This formula describes exact solutions of the Grad-Shafranov equation with up-down symmetry.
    The constants $c_i$ are determined from boundary constraints on $\psi$, that are derived from
    the following analytic model for a smooth, elongated "D" shaped cross section:
    ```math
    \begin{aligned}
    x &= 1 + \epsilon \, \cos (\tau + \delta_0 \, \sin \tau) , \\
    y &= \epsilon \kappa \, \sin (\tau) ,
    \end{aligned}
    ```
    where $0 \leq \tau < 2 \pi$, $\epsilon = a / R_0$ is the inverse aspect ratio, $\kappa$ the elongation,
    and $\sin \delta_0 = \delta$ is the triangularity.

    Defining three test points, namely
    - the high point $(1 - \delta \epsilon, \kappa \epsilon)$,
    - the inner equatorial point $(1 - \epsilon, 0)$,
    - and the outer equatorial point $(1 + \epsilon, 0)$,

    the following geometric constraints can be posed on the solution:
    ```math
    \begin{aligned}
    \psi (1 + \epsilon, 0) &= 0 , \\
    \psi (1 - \epsilon, 0) &= 0 , \\
    \psi (1 - \delta \epsilon, \kappa \epsilon) &= 0 , \\
    \psi_{x} (1 - \delta \epsilon, \kappa \epsilon) &= 0 , \\
    \psi_{yy} (1 + \epsilon, 0) &= - N_1 \psi_{x} (1 + \epsilon, 0) , \\
    \psi_{yy} (1 - \epsilon, 0) &= - N_2 \psi_{x} (1 - \epsilon, 0) , \\
    \psi_{xx} (1 - \delta \epsilon, \kappa \epsilon) &= - N_3 \psi_y (1 - \delta \epsilon, \kappa \epsilon) .
    \end{aligned}
    ```

    The first three equations define the three test points, the fourth equations enforces the high
    point to be a maximum, and the last three equations define the curvature at the test points.

    The coefficients $N_j$ can be found from the analytic model cross section as
    ```math
    \begin{aligned}
    N_1 &= \left[ \frac{d^2 x}{dy^2} \right]_{\tau = 0}     = - \frac{(1 + \delta_0)^2}{\epsilon \kappa^2} , \\
    N_2 &= \left[ \frac{d^2 x}{dy^2} \right]_{\tau = \pi}   = \hphantom{-} \frac{(1 - \delta_0)^2}{\epsilon \kappa^2} , \\
    N_3 &= \left[ \frac{d^2 x}{dy^2} \right]_{\tau = \pi/2} = - \frac{\kappa}{\epsilon \, \cos^2 \delta_0} .
    \end{aligned}
    ```

    For a given value of the constant $a$ above conditions reduce to a set of seven linear
    inhomogeneous algebraic equations for the unknown $c_i$, which can easily be solved.

    Parameters:
    * `R₀`: position of magnetic axis
    * `B₀`: B-field at magnetic axis
    * `ϵ`:  inverse aspect ratio
    * `κ`:  elongation
    * `δ`:  triangularity
    * `α`:  free constant, determined to match a given beta value
    """
    struct SolovevEquilibrium{T <: Number} <: AbstractSolovevEquilibrium
        name::String
        R₀::T
        B₀::T
        ϵ::T
        κ::T
        δ::T
        α::T
        c::Vector{T}

        function SolovevEquilibrium{T}(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, α::T, c::Vector{T}) where T <: Number
            new("Solovev Equilibrium", R₀, B₀, ϵ, κ, δ, α, c)
        end
    end

    function SolovevEquilibrium(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, α::T) where T <: Number

        n = 7
        A = zeros(n,n+1)

        x₁, x₂, x₃ = symbols("x₁, x₂, x₃")
        x = [x₁, x₂, x₃]
        ψ = [ψ₁(x), ψ₂(x), ψ₃(x), ψ₄(x), ψ₅(x), ψ₆(x), ψ₇(x), -ψ₀(x,α)]

        for i in axes(A,2)
            A[1,i] = N(subs(ψ[i], x[1]=>1+ϵ, x[2]=>0))
            A[2,i] = N(subs(ψ[i], x[1]=>1-ϵ, x[2]=>0))
            A[3,i] = N(subs(ψ[i], x[1]=>1-δ*ϵ, x[2]=>κ*ϵ))
            A[4,i] = N(subs(diff(ψ[i], x[1]), x[1]=>1-δ*ϵ, x[2]=>κ*ϵ))
            A[5,i] = N(subs(diff(ψ[i], x[2], 2), x[1]=>1+ϵ, x[2]=>0) - (1 + asin(δ))^2 / (ϵ * κ^2) * subs(diff(ψ[i], x[1]), x[1]=>1+ϵ, x[2]=>0))
            A[6,i] = N(subs(diff(ψ[i], x[2], 2), x[1]=>1-ϵ, x[2]=>0) + (1 - asin(δ))^2 / (ϵ * κ^2) * subs(diff(ψ[i], x[1]), x[1]=>1-ϵ, x[2]=>0))
            A[7,i] = N(subs(diff(ψ[i], x[1], 2), x[1]=>1-δ*ϵ, x[2]=>κ*ϵ) - κ / (ϵ * (1 - δ^2)) * subs(diff(ψ[i], x[2]), x[1]=>1-δ*ϵ, x[2]=>κ*ϵ))
        end

        c = A[1:n,1:n] \ A[1:n,n+1]

        SolovevEquilibrium{T}(R₀, B₀, ϵ, κ, δ, α, c)
    end

    function init(R₀, B₀, ϵ, κ, δ, α)
        SolovevEquilibrium(R₀, B₀, ϵ, κ, δ, α)
    end

    macro code(R₀, B₀, ϵ, κ, δ, α)
        code(init(R₀, B₀, ϵ, κ, δ, α); escape=true)
    end

    SolovevEquilibriumITER() = SolovevEquilibrium(6.2, 5.3, 0.32,  1.7, 0.33, -0.155)
    # SolovevEquilibriumTFTR() = SolovevEquilibrium(2.5, 5.6, 0.345, 1.0, 0.0,  )
    # SolovevEquilibriumJET()  = SolovevEquilibrium(3.0, 3.6, 0.333, 1.7, 0.25, )
    SolovevEquilibriumNSTX() = SolovevEquilibrium(0.85, 0.30, 0.78, 2.00, 0.35, 1.0)
    # SolovevEquilibriumMAST() = SolovevEquilibrium(0.85, 0.52, 0.77, 2.45, 0.50,  )
    SolovevEquilibriumFRC()  = SolovevEquilibrium(0.0, 0.0, 0.99, 10., 0.7, 0.0)
    # SolovevEquilibriumFRC2() = SolovevEquilibrium(0.0, 0.0, 1.00, 10., 1.0, 0.0)

    function Base.show(io::IO, equ::SolovevEquilibrium)
        print(io, "SolovevEquilibrium Equilibrium with\n")
        print(io, "  R₀ = ", equ.R₀, "\n")
        print(io, "  B₀ = ", equ.B₀, "\n")
        print(io, "  ϵ  = ", equ.ϵ,  "\n")
        print(io, "  κ  = ", equ.κ,  "\n")
        print(io, "  δ  = ", equ.δ,  "\n")
        print(io, "  α  = ", equ.α)
    end

    function ElectromagneticFields.A₃(x::AbstractArray{T,1}, equ::SolovevEquilibrium) where {T <: Number}
        ( ψ₀(x, equ.α) + equ.c[1] * ψ₁(x)
                       + equ.c[2] * ψ₂(x)
                       + equ.c[3] * ψ₃(x)
                       + equ.c[4] * ψ₄(x)
                       + equ.c[5] * ψ₅(x)
                       + equ.c[6] * ψ₆(x)
                       + equ.c[7] * ψ₇(x) )
    end

    @doc raw"""
    Axisymmetric Solov'ev equilibra with X-point in (R/R₀,Z/R₀,phi) coordinates.
    Based on Cerfon & Freidberg, Physics of Plasmas 17, 032502, 2010,
        and Freidberg, Ideal Magnetohydrodynamics, 2014.

    The covariant components of the vector potential are given by
    ```math
    A (x, y, \phi) = \left( \frac{B_0 R_0}{2} \, \frac{y}{x} , \, - \frac{B_0 R_0}{2} \, \ln x , \, \psi(x,y) \right)^T ,
    ```

    with $x = R/R_0$ and $y = Z/R_0$.
    The normalised poloidal flux $\psi$ is given by
    ```math
    \psi (x,y) = \psi_0 + \sum \limits_{i=1}^{12} c_i \psi_i (x,y) ,
    ```

    with
    ```math
    \begin{aligned}
    \psi_{0}  &= \frac{x^4}{8} + \alpha \left( \frac{1}{2} x^2 \, \ln x - \frac{x^4}{8} \right) , \\
    \psi_{1}  &= 1 , \\
    \psi_{2}  &= x^2 , \\
    \psi_{3}  &= y^2 - x^2 \, \ln x , \\
    \psi_{4}  &= x^4 - 4 x^2 y^2 , \\
    \psi_{5}  &= 2 y^4 9 y^2 x^2 + 3 x^4 \, \ln x - 12 x^2 y^2 \, \ln x , \\
    \psi_{6}  &= x^6 - 12 x^4 y^2 + 8 x^2 y^4 , \\
    \psi_{7}  &= 8 y^6 - 140 y^4 x^2 + 75 y^2 x^4 - 15 x^6 \, \ln x + 180 x^4 y^2 \, \ln x - 120 x^2 y^4 \, \ln x , \\
    \psi_{8}  &= y , \\
    \psi_{9}  &= y x^2 , \\
    \psi_{10} &= y^3 - 3 y x^2 \, \ln x , \\
    \psi_{11} &= 3 y x^4 - 4 y^3 x^2 , \\
    \psi_{12} &= 8 y^5 - 45 y x^4 - 80 y^3 x^2 \, \ln x + 60 y x^4 \, \ln x .
    \end{aligned}
    ```

    This formula describes exact solutions of the Grad-Shafranov equation with up-down asymmetry.
    The constants $c_i$ are determined from boundary constraints on $\psi$, that are derived from
    the following analytic model for a smooth, elongated "D" shaped cross section:
    ```math
    \begin{aligned}
    x &= 1 + \epsilon \, \cos (\tau + \arcsin \delta \, \sin \tau) , \\
    y &= \epsilon \kappa \, \sin (\tau) ,
    \end{aligned}
    ```
    where $0 \leq \tau < 2 \pi$, $\epsilon = a / R_0$ is the inverse aspect ratio, $\kappa$ the elongation,
    and $\sin \delta_0 = \delta$ is the triangularity.

    Defining four test points, namely
    - the high point $(1 - \delta \epsilon, \kappa \epsilon)$,
    - the inner equatorial point $(1 - \epsilon, 0)$,
    - and the outer equatorial point $(1 + \epsilon, 0)$,
    - the position of the X-point $(x_{\mathrm{sep}}, y_{\mathrm{sep}})$,

    the following geometric constraints can be posed on the solution:
    ```math
    \begin{aligned}
    \psi (1 + \epsilon, 0) &= 0 , \\
    \psi (1 - \epsilon, 0) &= 0 , \\
    \psi (1 - \delta \epsilon, \kappa \epsilon) &= 0 , \\
    \psi (x_{\mathrm{sep}}, y_{\mathrm{sep}}) &= 0 , \\
    \psi_{y} (1 + \epsilon, 0) &= 0 , \\
    \psi_{y} (1 - \epsilon, 0) &= 0 , \\
    \psi_{x} (1 - \delta \epsilon, \kappa \epsilon) &= 0 , \\
    \psi_{x} (x_{\mathrm{sep}}, y_{\mathrm{sep}}) &= 0 , \\
    \psi_{y} (x_{\mathrm{sep}}, y_{\mathrm{sep}}) &= 0 , \\
    \psi_{yy} (1 + \epsilon, 0) &= - N_1 \psi_{x} (1 + \epsilon, 0) , \\
    \psi_{yy} (1 - \epsilon, 0) &= - N_2 \psi_{x} (1 - \epsilon, 0) , \\
    \psi_{xx} (1 - \delta \epsilon, \kappa \epsilon) &= - N_3 \psi_y (1 - \delta \epsilon, \kappa \epsilon) .
    \end{aligned}
    ```

    The first four equations define the four test points, the fifth and sixth equations define the up-down
    symmetry, the seventh equations enforces the high point to be a maximum, the eighth and ninth eqaution
    set the $x$- and $y$-components of the magnetic field at the X-point to zero, and the last three
    equations define the curvature at the first three test points.

    The coefficients $N_j$ can be found from the analytic model cross section as
    ```math
    \begin{aligned}
    N_1 &= \left[ \frac{d^2 x}{dy^2} \right]_{\tau = 0}     = - \frac{(1 + \delta_0)^2}{\epsilon \kappa^2} , \\
    N_2 &= \left[ \frac{d^2 x}{dy^2} \right]_{\tau = \pi}   = \hphantom{-} \frac{(1 - \delta_0)^2}{\epsilon \kappa^2} , \\
    N_3 &= \left[ \frac{d^2 x}{dy^2} \right]_{\tau = \pi/2} = - \frac{\kappa}{\epsilon \, \cos^2 \delta_0} .
    \end{aligned}
    ```

    For a given value of the constant $a$ above conditions reduce to a set of seven linear
    inhomogeneous algebraic equations for the unknown $c_i$, which can easily be solved.

    Parameters:
    * `R₀`: position of magnetic axis
    * `B₀`: B-field at magnetic axis
    * `ϵ`:  inverse aspect ratio
    * `κ`:  elongation
    * `δ`:  triangularity
    * `α`:  free constant, determined to match a given beta value
    * `xsep`: x position of the X point
    * `ysep`: y position of the X point
    """
    struct SolovevXpointEquilibrium{T <: Number} <: AbstractSolovevEquilibrium
        name::String
        R₀::T
        B₀::T
        ϵ::T
        κ::T
        δ::T
        α::T
        xsep::T
        ysep::T
        c::Vector{T}

        function SolovevXpointEquilibrium{T}(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, α::T, xsep::T, ysep::T, c::Vector{T}) where T <: Number
            new("Solovev Equilibrium with X-point", R₀, B₀, ϵ, κ, δ, α, xsep, ysep, c)
        end
    end

    function SolovevXpointEquilibrium(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, α::T, xsep::T, ysep::T) where T <: Number

        n = 12
        A = zeros(n,n+1)

        x₁, x₂, x₃ = symbols("x₁, x₂, x₃")
        x = [x₁, x₂, x₃]
        ψ = [ψ₁(x), ψ₂(x), ψ₃(x), ψ₄(x), ψ₅(x), ψ₆(x), ψ₇(x), ψ₈(x), ψ₉(x), ψ₁₀(x), ψ₁₁(x), ψ₁₂(x), -ψ₀(x,α)]

        for i in axes(A,2)
            A[ 1,i] = N(subs(ψ[i], x[1]=>1+ϵ, x[2]=>0))
            A[ 2,i] = N(subs(ψ[i], x[1]=>1-ϵ, x[2]=>0))
            A[ 3,i] = N(subs(ψ[i], x[1]=>1-δ*ϵ, x[2]=>κ*ϵ))
            A[ 4,i] = N(subs(ψ[i], x[1]=>xsep, x[2]=>ysep))
            A[ 5,i] = N(subs(diff(ψ[i], x[2]), x[1]=>1+ϵ, x[2]=>0))
            A[ 6,i] = N(subs(diff(ψ[i], x[2]), x[1]=>1-ϵ, x[2]=>0))
            A[ 7,i] = N(subs(diff(ψ[i], x[1]), x[1]=>1-δ*ϵ, x[2]=>κ*ϵ))
            A[ 8,i] = N(subs(diff(ψ[i], x[1]), x[1]=>xsep, x[2]=>ysep))
            A[ 9,i] = N(subs(diff(ψ[i], x[2]), x[1]=>xsep, x[2]=>ysep))
            A[10,i] = N(subs(diff(ψ[i], x[2], 2), x[1]=>1+ϵ, x[2]=>0) - (1 + asin(δ))^2 / (ϵ * κ^2) * subs(diff(ψ[i], x[1]), x[1]=>1+ϵ, x[2]=>0))
            A[11,i] = N(subs(diff(ψ[i], x[2], 2), x[1]=>1-ϵ, x[2]=>0) + (1 - asin(δ))^2 / (ϵ * κ^2) * subs(diff(ψ[i], x[1]), x[1]=>1-ϵ, x[2]=>0))
            A[12,i] = N(subs(diff(ψ[i], x[1], 2), x[1]=>1-δ*ϵ, x[2]=>κ*ϵ) - κ / (ϵ * (1 - δ^2)) * subs(diff(ψ[i], x[2]), x[1]=>1-δ*ϵ, x[2]=>κ*ϵ))
        end

        c = A[1:n,1:n] \ A[1:n,n+1]

        SolovevXpointEquilibrium{T}(R₀, B₀, ϵ, κ, δ, α, xsep, ysep, c)
    end

    function SolovevDoubleXpointEquilibrium(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, α::T, xsep::T, ysep::T) where T <: Number

        n = 7
        A = zeros(n,n+1)

        x₁, x₂, x₃ = symbols("x₁, x₂, x₃")
        x = [x₁, x₂, x₃]
        ψ = [ψ₁(x), ψ₂(x), ψ₃(x), ψ₄(x), ψ₅(x), ψ₆(x), ψ₇(x), -ψ₀(x,α)]

        for i in axes(A,2)
            A[1,i] = N(subs(ψ[i], x[1]=>1+ϵ, x[2]=>0))
            A[2,i] = N(subs(ψ[i], x[1]=>1-ϵ, x[2]=>0))
            A[3,i] = N(subs(ψ[i], x[1]=>xsep, x[2]=>ysep))
            A[4,i] = N(subs(diff(ψ[i], x[1]), x[1]=>xsep, x[2]=>ysep))
            A[5,i] = N(subs(diff(ψ[i], x[2]), x[1]=>xsep, x[2]=>ysep))
            A[6,i] = N(subs(diff(ψ[i], x[2], 2), x[1]=>1+ϵ, x[2]=>0) - (1 + asin(δ))^2 / (ϵ * κ^2) * subs(diff(ψ[i], x[1]), x[1]=>1+ϵ, x[2]=>0))
            A[7,i] = N(subs(diff(ψ[i], x[2], 2), x[1]=>1-ϵ, x[2]=>0) + (1 - asin(δ))^2 / (ϵ * κ^2) * subs(diff(ψ[i], x[1]), x[1]=>1-ϵ, x[2]=>0))
        end

        c = vcat(A[1:n,1:n] \ A[1:n,n+1], zeros(5))

        SolovevXpointEquilibrium{T}(R₀, B₀, ϵ, κ, δ, α, xsep, ysep, c)
    end

    function init(R₀, B₀, ϵ, κ, δ, α, xsep, ysep, doublex=false)
        if doublex
            return SolovevDoubleXpointEquilibrium(R₀, B₀, ϵ, κ, δ, α, xsep, ysep)
        else
            return SolovevXpointEquilibrium(R₀, B₀, ϵ, κ, δ, α, xsep, ysep)
        end
    end

    macro code_xpoint(R₀, B₀, ϵ, κ, δ, α, xsep, ysep, doublex=false)
        code(SolovevXpointEquilibrium(R₀, B₀, ϵ, κ, δ, α, xsep, ysep, doublex); escape=true)
    end

    SolovevXpointEquilibriumITER() = SolovevXpointEquilibrium(6.2, 5.3, 0.32, 1.7, 0.33, -0.155, 0.88, -0.60)
    SolovevXpointEquilibriumNSTX() = SolovevXpointEquilibrium(0.85, 0.3, 0.78, 2.0, 0.35, -0.05, 0.70, -1.71)
    SolovevDoubleXpointEquilibriumNSTX() = SolovevDoubleXpointEquilibrium(0.85, 0.3, 0.78, 2.0, 0.35, 0.0, 0.70, -1.71)

    function Base.show(io::IO, equ::SolovevXpointEquilibrium)
        print(io, "Solovev Xpoint Equilibrium with\n")
        print(io, "  R₀ = ", equ.R₀, "\n")
        print(io, "  B₀ = ", equ.B₀, "\n")
        print(io, "  ϵ  = ", equ.ϵ,  "\n")
        print(io, "  κ  = ", equ.κ,  "\n")
        print(io, "  δ  = ", equ.δ,  "\n")
        print(io, "  α  = ", equ.α,  "\n")
        print(io, "  xsep  = ", equ.xsep, "\n")
        print(io, "  ysep  = ", equ.ysep)
    end

    function ElectromagneticFields.A₃(x::AbstractArray{T,1}, equ::SolovevXpointEquilibrium) where {T <: Number}
        ( ψ₀(x, equ.α) + equ.c[1]  * ψ₁(x)
                       + equ.c[2]  * ψ₂(x)
                       + equ.c[3]  * ψ₃(x)
                       + equ.c[4]  * ψ₄(x)
                       + equ.c[5]  * ψ₅(x)
                       + equ.c[6]  * ψ₆(x)
                       + equ.c[7]  * ψ₇(x)
                       + equ.c[8]  * ψ₈(x)
                       + equ.c[9]  * ψ₉(x)
                       + equ.c[10] * ψ₁₀(x)
                       + equ.c[11] * ψ₁₁(x)
                       + equ.c[12] * ψ₁₂(x) )
    end

    function ITER(; xpoint=false)
        if xpoint
            return SolovevXpointEquilibriumITER()
        else
            return SolovevEquilibriumITER()
        end
    end

    function NSTX(; xpoint=false)
        if xpoint
            return SolovevXpointEquilibriumNSTX()
        else
            return SolovevEquilibriumNSTX()
        end
    end

    function NSTXdoubleX()
        SolovevDoubleXpointEquilibriumNSTX()
    end

    function FRC()
        SolovevEquilibriumFRC()
    end

    macro code_iter(xpoint=false)
        code(ITER(xpoint=xpoint); escape=true)
    end

    macro code_iter_xpoint()
        code(ITER(xpoint=true); escape=true)
    end

    macro code_nstx(xpoint=false)
        code(NSTX(xpoint=xpoint); escape=true)
    end

    macro code_nstx_xpoint(doublex=false)
        if doublex
            equilibrium = NSTXdoubleX()
        else
            equilibrium = NSTX(xpoint=true)
        end
        code(equilibrium; escape=true)
    end

    macro code_nstx_double_xpoint()
        code(NSTXdoubleX(); escape=true)
    end

    macro code_frc()
        code(FRC(); escape=true)
    end

    @recipe function f(equ::SolovevEquilibrium;
                       nx = 100, ny = 120, nτ = 200, levels = 50, size = (300,400), aspect_ratio = :equal,
                       xlims = ( 0.50,  1.50),
                       ylims = (-0.75, +0.75))

        xgrid = LinRange(xlims[1], xlims[2], nx)
        zgrid = LinRange(ylims[1], ylims[2], ny)
        pot   = [ElectromagneticFields.A₃([xgrid[i], zgrid[j], 0.0], equ) / xgrid[i] for i in eachindex(xgrid), j in eachindex(zgrid)]

        τ = LinRange(0, 2π, nτ)
        boundary_X = 1 .+ equ.ϵ .* cos.(τ .+ asin(equ.δ) .* sin.(τ) )
        boundary_Y = equ.ϵ .* equ.κ .* sin.(τ)

        aspect_ratio := aspect_ratio
        size   := size
        xlims  := xlims
        ylims  := ylims
        levels := levels
        legend := :none

        @series begin
            seriestype := :contour
            (xgrid, zgrid, pot')
        end

        @series begin
            seriestype  := :path
            seriescolor := :red
            linewidth := 3
            (boundary_X, boundary_Y)
        end
    end

    @recipe function f(equ::SolovevXpointEquilibrium;
                       nx = 100, ny = 120, levels = 50, size = (300,400),
                       xlims = ( 0.50,  1.50),
                       ylims = (-0.75, +0.75))

        xgrid = LinRange(xlims[1], xlims[2], nx)
        zgrid = LinRange(ylims[1], ylims[2], ny)
        pot   = [ElectromagneticFields.A₃([xgrid[i], zgrid[j], 0.0], equ) / xgrid[i] for i in eachindex(xgrid), j in eachindex(zgrid)]

        seriestype := :contour
        aspect_ratio := :equal
        size   := size
        xlims  := xlims
        ylims  := ylims
        levels := levels
        legend := :none

        (xgrid, zgrid, pot')
    end

end

module SolovevFRC
    import ..Solovev: @code_frc, FRC
    export @code, init
    var"@code" = var"@code_frc"
    init = FRC
end

module SolovevITER
    import ..Solovev: @code_iter, ITER
    export @code, init
    var"@code" = var"@code_iter"
    init = ITER
end

module SolovevITERwXpoint
    import ..Solovev: @code_iter_xpoint, ITER
    export @code, init
    var"@code" = var"@code_iter_xpoint"
    init() = ITER(xpoint=true)
end

module SolovevNSTX
    import ..Solovev: @code_nstx, NSTX
    export @code, init
    var"@code" = var"@code_nstx"
    init = NSTX
end

module SolovevNSTXwXpoint
    import ..Solovev: @code_nstx_xpoint, NSTX
    export @code, init
    var"@code" = var"@code_nstx_xpoint"
    init() = NSTX(xpoint=true)
end

module SolovevNSTXwDoubleXpoint
    import ..Solovev: @code_nstx_double_xpoint, NSTXdoubleX
    export @code, init
    var"@code" = var"@code_nstx_double_xpoint"
    init = NSTXdoubleX
end


text:
"""
    Laplace(μ,β)

The *Laplace distribution* with location `μ` and scale `β` has probability density function

```math
f(x; \\mu, \\beta) = \\frac{1}{2 \\beta} \\exp \\left(- \\frac{|x - \\mu|}{\\beta} \\right)
```

```julia
Laplace()       # Laplace distribution with zero location and unit scale, i.e. Laplace(0, 1)
Laplace(μ)      # Laplace distribution with location μ and unit scale, i.e. Laplace(μ, 1)
Laplace(μ, β)   # Laplace distribution with location μ and scale β

params(d)       # Get the parameters, i.e., (μ, β)
location(d)     # Get the location parameter, i.e. μ
scale(d)        # Get the scale parameter, i.e. β
```

External links

* [Laplace distribution on Wikipedia](http://en.wikipedia.org/wiki/Laplace_distribution)

"""
struct Laplace{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    θ::T
    Laplace{T}(µ::T, θ::T) where {T} = new{T}(µ, θ)
end

function Laplace(μ::T, θ::T; check_args=true) where {T <: Real}
    check_args && @check_args(Laplace, θ > zero(θ))
    return Laplace{T}(μ, θ)
end

Laplace(μ::Real, θ::Real) = Laplace(promote(μ, θ)...)
Laplace(μ::Integer, θ::Integer) = Laplace(float(μ), float(θ))
Laplace(μ::T) where {T <: Real} = Laplace(μ, one(T))
Laplace() = Laplace(0.0, 1.0, check_args=false)

const Biexponential = Laplace

@distr_support Laplace -Inf Inf

#### Conversions
function convert(::Type{Laplace{T}}, μ::S, θ::S) where {T <: Real, S <: Real}
    Laplace(T(μ), T(θ))
end
function convert(::Type{Laplace{T}}, d::Laplace{S}) where {T <: Real, S <: Real}
    Laplace(T(d.μ), T(d.θ), check_args=false)
end

#### Parameters

location(d::Laplace) = d.μ
scale(d::Laplace) = d.θ
params(d::Laplace) = (d.μ, d.θ)
@inline partype(d::Laplace{T}) where {T<:Real} = T

#### Statistics

mean(d::Laplace) = d.μ
median(d::Laplace) = d.μ
mode(d::Laplace) = d.μ

var(d::Laplace) = 2d.θ^2
std(d::Laplace) = sqrt2 * d.θ
skewness(d::Laplace{T}) where {T<:Real} = zero(T)
kurtosis(d::Laplace{T}) where {T<:Real} = 3one(T)

entropy(d::Laplace) = log(2d.θ) + 1

#### Evaluations

zval(d::Laplace, x::Real) = (x - d.μ) / d.θ
xval(d::Laplace, z::Real) = d.μ + z * d.θ

pdf(d::Laplace, x::Real) = exp(-abs(zval(d, x))) / 2scale(d)
logpdf(d::Laplace, x::Real) = - (abs(zval(d, x)) + log(2scale(d)))

cdf(d::Laplace, x::Real) = (z = zval(d, x); z < 0 ? exp(z)/2 : 1 - exp(-z)/2)
ccdf(d::Laplace, x::Real) = (z = zval(d, x); z > 0 ? exp(-z)/2 : 1 - exp(z)/2)
logcdf(d::Laplace, x::Real) = (z = zval(d, x); z < 0 ? loghalf + z : loghalf + log2mexp(-z))
logccdf(d::Laplace, x::Real) = (z = zval(d, x); z > 0 ? loghalf - z : loghalf + log2mexp(z))

quantile(d::Laplace, p::Real) = p < 1/2 ? xval(d, log(2p)) : xval(d, -log(2(1 - p)))
cquantile(d::Laplace, p::Real) = p > 1/2 ? xval(d, log(2(1 - p))) : xval(d, -log(2p))
invlogcdf(d::Laplace, lp::Real) = lp < loghalf ? xval(d, logtwo + lp) : xval(d, -(logtwo + log1mexp(lp)))
invlogccdf(d::Laplace, lp::Real) = lp > loghalf ? xval(d, logtwo + log1mexp(lp)) : xval(d, -(logtwo + lp))

function gradlogpdf(d::Laplace, x::Real)
    μ, θ = params(d)
    x == μ && error("Gradient is undefined at the location point")
    g = 1 / θ
    x > μ ? -g : g
end
heslogpdf(d::Laplace{T}, x::Real) where {T<:Real} = zero(T)

function mgf(d::Laplace, t::Real)
    st = d.θ * t
    exp(t * d.μ) / ((1 - st) * (1 + st))
end
function cf(d::Laplace, t::Real)
    st = d.θ * t
    cis(t * d.μ) / (1+st*st)
end

#### Sampling

rand(rng::AbstractRNG, d::Laplace) =
    d.μ + d.θ*randexp(rng)*ifelse(rand(rng, Bool), 1, -1)

#### Fitting

function fit_mle(::Type{<:Laplace}, x::Array)
    xc = copy(x)
    a = median!(xc)
    Laplace(a, StatsBase.mad!(xc, center=a))
end


text:
#
# DEM simulatior for Julia lang. practice
#

using Printf
using Plots
ENV["GKSwstype"] = "nul"

const global num_active_particles = 20
const global num_boundary_particles = 100*100
const global num_particles = num_active_particles + num_boundary_particles
const global particle_type = 2 # active particle and boundary material

const global ms = 1.0  # mass of particle
const global rs = 0.01 # radius of particle
const global Is = fill(2.0/5.0*ms*rs*rs, num_particles, 3)  # moment of Inertia of particle
const global g = [0, 0, -9.81] # gravity
const global pID = fill(1, num_particles)

vb = [-1,+1];

# initialize particle state

xs = rand(Float64, num_particles, 3) .- 0.5 # position vector xyz:[-0.5,0.5]
n = num_active_particles
for i in 1:100, j in 1:100
   global n +=1
   xs[n,1] = -0.5 + (1.0/100) * (i-0.5)
   xs[n,2] = -0.5 + (1.0/100) * (j-0.5)
   xs[n,3] = -0.5
end
vs = zeros(Float64, num_particles, 3)       # velocity vector
ws = zeros(Float64, num_particles, 3)       # rotational velocity vector
ts = zeros(Float64, num_particles, 3)       # rotation vector

force  = zeros(Float64, num_particles, 3)
torque = zeros(Float64, num_particles, 3)

xsnew = xs
vsnew = vs
wsnew = ws
tsnew = ts

# time stepping (symplectic Euler scheme)
deltaT = 0.01 # s
endTime = 1.0  # s
timeSteps = Integer(endTime/deltaT)

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

function forces_update(force, state="init")
   if state == "init"
     # @printf("init \n")
     force  = zeros(Float64, num_particles, 3)
   else
     # @printf("update \n")
     force[1:num_active_particles,1] .= ms*g[1]
     force[1:num_active_particles,2] .= ms*g[2]
     force[1:num_active_particles,3] .= ms*g[3]
   end
end

function torque_update(torque, state="init")
   if state == "init"
     #@printf("init \n")
     torque  = zeros(Float64, num_particles, 3)
   else
     #@printf("update \n")
     torque[1:num_particles,1] .= 0.001
     torque[1:num_particles,2] .= 0.001
     torque[1:num_particles,3] .= 0.001
   end
end

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

anim = Animation()

forces_update(force, "init")
torque_update(torque, "init")

for i in 1:timeSteps

  # forces and torque  update
  forces_update(force,"runtime")
  torque_update(torque,"runtime")

  # velocities update
  global vsnew = vs + force/ms * deltaT

  for n=1:num_particles
    global wsnew[n,1] = ws[n,1] + torque[n,1]./Is[n,1]*deltaT
    global wsnew[n,2] = ws[n,2] + torque[n,2]./Is[n,2]*deltaT
    global wsnew[n,3] = ws[n,3] + torque[n,3]./Is[n,3]*deltaT
  end

  # position and rotation angle update
  global xsnew = xs + vsnew * deltaT
  global tsnew = ts + wsnew * deltaT

  #@printf("%d %f  \n", i, minimum(xsnew))

  # time increment
  global xs = xsnew
  global vs = vsnew
  global ws = wsnew
  global ts = tsnew

  # Plots
  plot(xs[:,1],xs[:,2], xs[:,3], xlim=vb, ylim=vb, zlim=vb, seriestype=:scatter, markersize=6, title="Particle", legend=:none)
  #display(plt3d)
  frame(anim)

end
gif(anim)


text:
# ---
# title: 514. Freedom Trail
# id: problem514
# author: AquaIndigo
# date: 2020-11-11
# difficulty: Hard
# categories: Divide and Conquer, Dynamic Programming, Depth-first Search
# link: <https://leetcode.com/problems/freedom-trail/description/>
# hidden: true
# ---
#
# In the video game Fallout 4, the quest "Road to Freedom" requires players to
# reach a metal dial called the "Freedom Trail Ring", and use the dial to spell
# a specific keyword in order to open the door.
#
# Given a string **ring** , which represents the code engraved on the outer ring
# and another string **key** , which represents the keyword needs to be spelled.
# You need to find the **minimum** number of steps in order to spell all the
# characters in the keyword.
#
# Initially, the first character of the **ring** is aligned at 12:00 direction.
# You need to spell all the characters in the string **key** one by one by
# rotating the ring clockwise or anticlockwise to make each character of the
# string **key** aligned at 12:00 direction and then by pressing the center
# button.
#
# At the stage of rotating the ring to spell the key character **key[i]** :
#
#   1. You can rotate the **ring** clockwise or anticlockwise **one place** , which counts as 1 step. The final purpose of the rotation is to align one of the string **ring 's** characters at the 12:00 direction, where this character must equal to the character **key[i]**.
#   2. If the character **key[i]** has been aligned at the 12:00 direction, you need to press the center button to spell, which also counts as 1 step. After the pressing, you could begin to spell the next character in the key (next stage), otherwise, you've finished all the spelling.
#
# **Example:**
#
# ![](https://assets.leetcode.com/uploads/2018/10/22/ring.jpg)
#
#
#     Input: ring = "godding", key = "gd"
#     Output: 4
#     Explanation:
#     For the first key character 'g', since it is already in place, we just need 1 step to spell this character.
#     For the second key character 'd', we need to rotate the ring "godding" anticlockwise by two steps to make it become "ddinggo".
#     Also, we need 1 more step for spelling.
#     So the final output is 4.
#
#
# **Note:**
#
#   1. Length of both ring and **key** will be in range 1 to 100.
#   2. There are only lowercase letters in both strings and might be some duplcate characters in both strings.
#   3. It's guaranteed that string **key** could always be spelled by rotating the string **ring**.
#
#
## @lc code=start
using LeetCode

function find_rotate_steps(ring::String, key::String)
    lenr, lenk = length(ring), length(key)
    dp1 = fill(typemax(Int), lenr)
    dp2 = dp1[:]
    mp = [Int[] for i in 1:26]
    for i in 1:lenr
        push!(mp[ring[i] - 'a' + 1], i)
    end
    for idx in mp[key[1] - 'a' + 1]
        dp1[idx] = min(idx - 1, lenr - idx)
    end
    for i in 2:lenk
        for idx1 in mp[key[i] - 'a' + 1]
            for idx2 in mp[key[i - 1] - 'a' + 1]
                dp2[idx1] = min(
                    dp2[idx1], dp1[idx2] + min(abs(idx2 - idx1), lenr - abs(idx2 - idx1))
                )
            end
        end
        dp1, dp2 = dp2, dp1
        dp2 .= typemax(Int)
    end
    return minimum(dp1) + lenk
end
## @lc code=end


text:
function compare_images(testfn::AbstractString, reffn::AbstractString; sigma = [1,1], eps = 0.02)

    result = VisualTestResult(testfn, nothing, reffn, nothing, PROCESSING_ERROR, 1.0, nothing)

    # load the test image... if we error, return immediately
    try
        result.testImage = Images.load(testfn)
    catch err
        result.err = err
        return result
    end

    # load the ref image... if we error, return immediately
    try
        result.referenceImage = Images.load(reffn)
    catch err
        result.err = err
        return result
    end

    # we loaded both images, now do the comparison
    try

        # run the comparison test... a difference will throw an error
        # NOTE: sigma is a 2-length vector with x/y values for the number of pixels
        #       to blur together when comparing images
        result.diff = Images.test_approx_eq_sigma_eps(result.testImage, result.referenceImage, sigma, eps)

        # we passed!
        info("Reference image $reffn matches.  Difference: $(result.diff)")

        result.status = if result.diff == 0
          EXACT_MATCH
        else
          CLOSE_MATCH
        end

    catch err

        warn("Got error: $err")
        result.err = err
        result.status = DOES_NOT_MATCH

        # HACK: this will fail if the Images error message changes
        if typeof(err) == ErrorException
          msg = split(err.msg)
          if length(msg) > 2 && msg[1] == "Arrays" && msg[2] == "differ."
            result.diff = parse(Float64, msg[end-2])
          end
        end

    end

    result
end

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

function test_images(testfn::AbstractString, reffn::AbstractString; popup=isinteractive(), newfn = reffn, kw...)
    result = compare_images(testfn, reffn; kw...)

    if !success(result)
        warn("Image did not match reference image $reffn. err: $(result.err)")

        if popup
            # open a popup and give us a chance to examine the images,
            # then ask to replace the reference
            warn("Should we make this the new reference image?")
            replace_refimg_dialog(testfn, reffn; reffn_new = newfn)
        end
    end

    result
end

function test_images(visualtest::VisualTest; popup=isinteractive(), kw...)
    visualtest.testFunction(visualtest.testFilename, visualtest.args...; visualtest.kw...)
    test_images(visualtest.testFilename, visualtest.referenceFilename; popup=popup, kw...)
end

function test_images(visualtests::AbstractVector{VisualTest}; popup=isinteractive(), kw...)
    VisualTestResult[test_images(vtest, popup=popup, kw...) for vtest in visualtests]
end


text:
# This file defines the ManpowerSimulation type. This type bundles all the
#   information of the simulation.

# The ManpowerSimulation type requires SimJulia, ResumableFunctions, and
#   Distributions.

# The ManpowerSimulation type requires all the other types.
requiredTypes = [ "prerequisite", "attrition", "state", "compoundState",
    "prerequisiteGroup", "recruitment", "retirement", "simulationReport" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

export ManpowerSimulation
type ManpowerSimulation
    # This is the name of the parameter configuration file. If this is an empty
    #   string, the simulation must be configured manually.
    catFileName::String
    parFileName::String

    # this is the name of the simulation database file.
    dbName::String

    # This is the name of this simulation instance.
    simName::String

    # This is a link to the database file holding the relevant tables for this
    #   simulation.
    simDB::SQLite.DB

    # These are the names of the personnel, history, and transition databases.
    personnelDBname::String
    historyDBname::String
    transitionDBname::String

    # This flag marks if the simulation has been properly initialised with the
    #   init function.
    isInitialised::Bool

    # This flag marks if the simulation has been properly initialised and has
    #   an empty results database.
    isVirgin::Bool

    # The ID key of the database.
    idKey::String

    # Target number of people in the simulation. If this is set to 0, there is
    #   no cap.
    personnelTarget::Int

    # The current number of active personnel members in the simulation.
    personnelSize::Int

    # The total number of personnel members in the simulation.
    resultSize::Int

    # The time between successive SQLite commits.
    commitFrequency::Float64

    # The additional attributes.
    initAttrList::Vector{PersonnelAttribute}
    otherAttrList::Vector{PersonnelAttribute}

    # The states.
    initStateList::Dict{State, Vector{Transition}}
    otherStateList::Dict{State, Vector{Transition}}
    stateList::Dict{String, State}

    # The compound states. The first is for the compound states from catalogue,
    #   so they can be properly processed at runtime and inserted into the
    #   the second list (which retains the component states of each compound
    #   state).
    compoundStates::Dict{String, State}
    compoundStateList::Dict{String, CompoundState}

    # The names of the transitions.
    transList::Dict{String, Int}

    # The recruitment schemes.
    recruitmentSchemes::Vector{Recruitment}

    # The default attrition scheme.
    defaultAttritionScheme::Attrition

    # Other attrition schemes.
    attritionSchemes::Dict{String, Attrition}

    # The time between two checks of the attrition process (purely technical)
    attritionTimeSkip::Float64

    # The retirement scheme.
    retirementScheme::Union{Void, Retirement}

    # A flag to track if the system is well defined.
    isWellDefined::Bool

    # SimJulia simulation.
    sim::Simulation

    # Processor time of the simulation.
    simTimeElapsed::Dates.Millisecond

    # Time of attrition execution processes.
    attrExecTimeElapsed::Dates.Millisecond

    # The priorities of the various simulation phases.
    phasePriorities::Dict{Symbol, Int}

    # The length of the simulation.
    simLength::Float64

    # The start date of the simulation.
    simStartDate::Date

    # The generated reports of the simulation.
    # simReports::Dict{Float64, SimulationReport}

    function ManpowerSimulation( ; dbName::String = "simDB",
        simName::String = "testSim" )

        newMPsim = new()
        newMPsim.catFileName = ""
        newMPsim.parFileName = ""
        newMPsim.simName = simName

        # This block creates the database file if necessary and opens a link to
        #   it.  XXX Persistent??
        tmpDBname = ( dbName == "" ) || endswith( dbName, ".sqlite" ) ?
            dbName : dbName * ".sqlite"

        if tmpDBname == ""
            println( "Results database kept in memory." )
        else
            println( "Database file \"$tmpDBname\" ",
                isfile( tmpDBname ) ? "exists" : "does not exist", "." )
        end  # if tmpDBname == ""

        newMPsim.dbName = tmpDBname
        newMPsim.simDB = SQLite.DB( tmpDBname )

        # This line ensures that foreign key logic works.
        SQLite.execute!( newMPsim.simDB, "PRAGMA foreign_keys = ON" )

        newMPsim.personnelDBname = "Personnel_$simName"
        newMPsim.historyDBname = "History_$simName"
        newMPsim.transitionDBname = "Transitions_$simName"

        # Check if databases are present and issue a warning if so.
        tmpTableList = SQLite.tables( newMPsim.simDB )[ :name ]

        if ( newMPsim.personnelDBname ∈ tmpTableList ) ||
                ( newMPsim.historyDBname ∈ tmpTableList )
            warn( "Results for a simulation called \"$(newMPsim.simName)\" already in database. These will be overwritten." )
        end  # if newMPsim.personnelDBname ∈ tmpTableList

        newMPsim.isInitialised = false
        newMPsim.isVirgin = false
        newMPsim.idKey = "id"
        newMPsim.personnelTarget = 0
        newMPsim.personnelSize = 0
        newMPsim.resultSize = 0
        newMPsim.commitFrequency = 1.0
        newMPsim.initAttrList = Vector{PersonnelAttribute}()
        newMPsim.otherAttrList = Vector{PersonnelAttribute}()
        newMPsim.initStateList = Dict{State, Vector{Transition}}()
        newMPsim.otherStateList = Dict{State, Vector{Transition}}()
        newMPsim.stateList = Dict{String, State}()
        newMPsim.compoundStates = Dict{String, State}()
        newMPsim.compoundStateList = Dict{String, CompoundState}()
        newMPsim.transList = Dict{String, Int}()
        newMPsim.recruitmentSchemes = Vector{Recruitment}()
        newMPsim.defaultAttritionScheme = Attrition( "default" )
        newMPsim.attritionSchemes = Dict{String, Attrition}()
        newMPsim.attritionTimeSkip = 1.0
        newMPsim.retirementScheme = Retirement()
        newMPsim.isWellDefined = true
        newMPsim.sim = Simulation()
        newMPsim.simTimeElapsed = Dates.Millisecond( 0 )
        newMPsim.attrExecTimeElapsed = Dates.Millisecond( 0 )
        newMPsim.phasePriorities = Dict( :attrCheck => typemin( Int ),
            :recruitment => typemin( Int ) + 1,
            # :transition => 30,
            :retirement => 40,
            :attrition => 50 )
        newMPsim.simLength = 1.0
        newMPsim.simStartDate = Date( now() )
        # newMPsim.simReports = Dict{Float64, SimulationReport}()
        return newMPsim

    end  # ManpowerSimulation( ; dbName, simName )

    function ManpowerSimulation( configFileName::String )

        newMPsim = ManpowerSimulation()
        tmpPath = Base.source_path()
        tmpPath = tmpPath isa Void ? "" : dirname( tmpPath )
        newMPsim.parFileName = joinpath( tmpPath, configFileName )
        initialiseFromExcel( newMPsim, configFileName )
        # initialise( newMPsim )
        return newMPsim

    end  # ManpowerSimulation( configFileName )

end  # type ManpowerSimulation