text:
""" The implementation is primarily based on a paper by Glover and Doyle, and
the code is best read with this docment at hand. All references to equations
in docstrings and comments are to the version of the paper given below:

  @article{glover1988state,
    title={State-space formulae for all stabilizing controllers that satisfy an
           H-infinity norm bound and relations to relations to risk sensitivity},
    author={Glover, Keith and Doyle, John C},
    journal={Systems & control letters},
    volume={11},
    number={3},
    pages={167--172},
    year={1988},
    publisher={Citeseer}
  }
"""

"""
    flag = hinfassumptions(P::ExtendedStateSpace; verbose=true)

Check the assumptions for using the γ-iteration synthesis in Theorem 1.
"""
function hinfassumptions(P::ExtendedStateSpace; verbose = true)

    A, B1, B2, C1, C2, D11, D12, D21, D22 = ssdata_e(P)

    # Check assumption A1
    if !_stabilizable(A, B2)
        verbose && @warn(
            "The system A is not stabilizable through B2, ",
            "violation of assumption A1."
        )
        return false
    end
    if !_detectable(A, C2)
        verbose && @warn(
            "The system A is not detectable through C2, ",
            "violation of assumption A1."
        )
        return false
    end

    # Check assumption A2
    if rank(D12) < size(D12, 2)
        verbose &&
            @warn("The matrix D12 does not have full rank, ", "violation of assumption A2. The full control signal must have nonzero penalty at infinite frequency.")
        return false
    end
    if rank(D21) < size(D21, 1)
        verbose &&
            @warn("The matrix D21 does not have full rank, ", "violation of assumption A2. The whole measurement vector y must be corrupted by noise at infinite frequency.")
        return false
    end

    # Check assumption A5
    if rank(A - B2 * pinv(D12) * C1) < size(A, 1)
        verbose && @warn(
            "The matrix (A - B2*D12^-*C1) does not have full",
            "rank, violation of assumption A5."
        )
        return false
    end
    # Check assumption A6
    if rank(A - B1 * pinv(D21) * C2) < size(A, 1)
        verbose && @warn(
            "The matrix (A - B1*D21Pinv*C2) does not ",
            "have full rank, violation of assumption A6."
        )
        return false
    end

    # All assumptions have passed, and we may proceed with the synthesis
    verbose && println("All assumtions are satisfied!")
    return true
end

"""
    flag = _stabilizable(A::AbstractMatrix, B::AbstractMatrix)

Applies the Hautus lemma to check if the pair is stabilizable
"""
function _stabilizable(A::AbstractMatrix, B::AbstractMatrix)
    eigValsA = eigvals(A)
    for ii in eachindex(eigValsA)
        if real(eigValsA[ii]) >= 0
            if rank([eigValsA[ii] * I - A B]) != size(A, 1)
                return false
            end
        end
    end
    return true
end

"""
    flag = _detectable(A::AbstractMatrix, C::AbstractMatrix)

Applies the Hautus lemma to check if the pair is detectable
"""
function _detectable(A::AbstractMatrix, C::AbstractMatrix)
    eigValsA = eigvals(A)
    for ii = 1:length(eigValsA)
        if real(eigValsA[ii]) >= 0
            if rank([eigValsA[ii] * I - A; C]) != size(A, 1)
                return false
            end
        end
    end
    return true
end

"""
    flag, K, γ, mats = hinfsynthesize(P::ExtendedStateSpace; maxIter = 20, interval = (2 / 3, 20), verbose = false, tolerance = 1.0e-10, γrel = 1.01)

Computes an H-infinity optimal controller K for an extended plant P such that
||F_l(P, K)||∞ < γ for the smallest possible γ given P. The routine is
known as the γ-iteration, and is based on the paper "State-space formulae for
all stabilizing controllers that satisfy an H∞-norm bound and relations to
risk sensitivity" by Glover and Doyle.

# Arguments:
- `maxIter`: Maximum number of γ iterations
- `interval`: The starting interval for the bisection.
- `verbose`: Print progress?
- `tolerance`: Stop when the interval is this small.
- `γrel`: If `γrel > 1`, the optimal γ will be found by γ iteration after which a controller will be designed for `γ = γopt * γrel`. It is often a good idea to design a slightly suboptimal controller, both for numerical reasons, but also since the optimal controller may contain very fast dynamics. If `γrel → ∞`, the computed controller will approach the 𝑯₂ optimal controller. Getting a mix between 𝑯∞ and 𝑯₂ properties is another reason to choose `γrel > 1`.
"""
function hinfsynthesize(
    P::ExtendedStateSpace;
    maxIter = 20,
    interval = (2 / 3, 20),
    verbose = false,
    tolerance = 1e-10,
    γrel = 1.01,
)

    # Transform the system into a suitable form
    P̄, Ltrans12, Rtrans12, Ltrans21, Rtrans21 = _transformp2pbar(P)

    # Run the γ iterations
    X∞Feasible, Y∞Feasible, F∞Feasible, H∞Feasible, γFeasible =
        _γiterations(P̄, maxIter, interval, verbose, tolerance)

    if !isempty(γFeasible)
        # Synthesize the controller and trnasform it back into the original coordinates

        if γrel > 1
            γFeasible *= γrel
            X∞Feasible, Y∞Feasible, F∞Feasible, H∞Feasible = _solvematrixequations(P̄, γFeasible)
        end

        K = _synthesizecontroller(
            P̄,
            X∞Feasible,
            Y∞Feasible,
            F∞Feasible,
            H∞Feasible,
            γFeasible,
            Ltrans12,
            Rtrans12,
            Ltrans21,
            Rtrans21,
        )

        # Return the controller, the optimal gain γ, and a true flag
        γ = γFeasible
        flag = true
    else
        # Return and empty controller, empty gain γ, and a false flag
        K = ss(0.0)
        γ = Inf
        flag = false
    end
    return flag, K, γFeasible, (X=X∞Feasible, Y=Y∞Feasible, F=F∞Feasible, H=H∞Feasible)
end

"""
    K = _synthesizecontroller(P::ExtendedStateSpace, Xinf, Yinf, F, H, γ, Ltrans12, Rtrans12, Ltrans21, Rtrans21)

Syntheize a controller by operating on the scaled state-space description of the
system (i.e., the state-space realization of `P̄`) using the solutions from the
γ-iterations. The controller is synthesized in the coordinates of `P̄`, and then
transformed back to the original coordinates by the linear transformations
Ltrans12, Rtrans12, Ltrans21 and Rtrans21.
"""
function _synthesizecontroller(
    P::ExtendedStateSpace,
    Xinf::AbstractMatrix,
    Yinf::AbstractMatrix,
    F::AbstractMatrix,
    H::AbstractMatrix,
    γ::Number,
    Ltrans12::AbstractMatrix,
    Rtrans12::AbstractMatrix,
    Ltrans21::AbstractMatrix,
    Rtrans21::AbstractMatrix,
)

    A   = P.A
    B1  = P.B1
    B2  = P.B2
    C1  = P.C1
    C2  = P.C2
    D11 = P.D11
    D12 = P.D12
    D21 = P.D21
    D22 = P.D22

    γ² = γ * γ

    B = [B1 B2]
    C = [C1; C2]

    # Dimensionality
    P1 = size(C1, 1)
    P2 = size(C2, 1)
    M1 = size(B1, 2)
    M2 = size(B2, 2)

    # Equation (11)
    # F11 = F[1:(M1-P2), :]
    F12 = F[(M1-P2+1):M1, :]
    F2 = F[(M1+1):(M1+M2), :]

    # Equation (12)
    # H11 = H[:, 1:(P1-M2)]
    H12 = H[:, (P1-M2+1):P1]
    H2 = H[:, (P1+1):(P1+P2)]

    # Definition of D in the assumptions section
    D1111 = D11[1:(P1-M2), 1:(M1-P2)]
    D1112 = D11[1:(P1-M2), (M1-P2+1):M1]
    D1121 = D11[(P1-M2+1):P1, 1:(M1-P2)]
    D1122 = D11[(P1-M2+1):P1, (M1-P2+1):M1]

    # Equation 19
    D11hat = ((-D1121 * D1111') / (γ² * I - D1111 * D1111')) * D1112 - D1122

    # Equation 20
    D12hatD12hat = I - (D1121 / (γ² * I - D1111' * D1111)) * D1121'
    _assertrealandpsd(D12hatD12hat; msg = " in equation (20)")
    D12hat = cholesky(Hermitian(D12hatD12hat)).L

    # Equation 21
    D21hatD21hat = I - (D1112' / (γ² * I - D1111 * D1111')) * D1112
    _assertrealandpsd(D21hatD21hat; msg = " in equation (21)")
    D21hat = cholesky(Hermitian(D21hatD21hat)).U

    # Equation 27
    Zinv = (I - Yinf * Xinf / γ²)

    # Equation 22
    B2hat = (B2 + H12) * D12hat

    # Equation 23 (using the inverse of 27)
    C2hat = -D21hat * (C2 + F12) / Zinv

    # Equation 24
    B1hat = -H2 + (B2hat / D12hat) * D11hat

    # Equation 25 (using the inverse of 27)
    C1hat = F2 / Zinv + (D11hat / D21hat) * C2hat

    # Equation 26
    Ahat = A + H * C + (B2hat / D12hat) * C1hat

    Ac = Ahat

    B1c = B1hat * Ltrans21
    B2c = B2hat

    C1c = Rtrans12 * C1hat
    C2c = C2hat

    Bc = [B1c B2c]
    Cc = [C1c; C2c]

    D11c = Rtrans12 * D11hat * Ltrans21
    # TODO implement loop shift for any system not satisfying A4
    # D12c = D12hat * Ltrans21
    # D21c = Rtrans12 * D21hat
    # D22c = zeros(size(D11hat))
    # Dc = [D11c D12c; D21c D22c]
    # return ss(Ac, Bc[:, 1:P2], Cc[1:M2, :], Dc[1:M2, 1:P2])

    return ss(Ac, Bc[:, 1:P2], Cc[1:M2, :], D11c)
end

"""
    _assertrealandpsd(A::AbstractMatrix, msg::AbstractString)

Check that a matrix is real and PSD - throw an error otherwise.
"""
function _assertrealandpsd(A::AbstractMatrix; msg = "")
    if any(real(eigvals(A)) .<= 0)
        error(string("The matrix", msg, " is not PSD."))
    end
    if any(imag(eigvals(A)) .!= 0)
        error(string("The matrix", msg, " is not real."))
    end
end

"""
    flag =  _checkfeasibility(Xinf, Yinf, γ, tolerance, iteration; verbose=true)

Check the feasibility of the computed solutions Xinf, Yinf and the algebraic
Riccatti equations, return true if the solution is valid, and false otherwise.
"""
function _checkfeasibility(
    Xinf::AbstractMatrix,
    Yinf::AbstractMatrix,
    γ::Number,
    tolerance::Number,
    iteration::Number;
    verbose = true,
)

    minXev = minimum(real.(eigvals(Xinf)))
    minYev = minimum(real.(eigvals(Yinf)))
    specrad = maximum(abs.(eigvals(Xinf * Yinf))) / (γ^2)

    if verbose
        if iteration == 1
            println("iteration, γ")
        end
        println(iteration, " ", γ)
    end

    if minXev < -tolerance
        # Failed test, eigenvalues of Xinf must be positive real
        return false
    end
    if minYev < -tolerance
        # Failed test, eigenvalues of Yinf must be positive real
        return false
    end
    if specrad > 1
        # Failed test, spectral radius of XY must be greater than γ squared
        return false
    end
    return true
end

"""
    solution = _solvehamiltonianare(H)

Solves a hamiltonian Alebraic Riccati equation using the Schur-decomposition,
for additional details, see

  @article{laub1979schur,
    title={A Schur method for solving algebraic Riccati equations},
    author={Laub, Alan},
    journal={IEEE Transactions on automatic control},
    volume={24},
    number={6},
    pages={913--921},
    year={1979},
    publisher={IEEE}
  }
"""
function _solvehamiltonianare(H)

    S = schur(H)
    S = ordschur(S, real.(S.values) .< 0)

    (m, n) = size(S.Z)
    U11 = S.Z[1:div(m, 2), 1:div(n, 2)]
    U21 = S.Z[div(m, 2)+1:m, 1:div(n, 2)]

    return U21 / U11
end

"""
    solution = _solvematrixequations(P::ExtendedStateSpace, γ::Number)

Solves the dual matrix equations in the γ-iterations (equations 7-12 in Doyle).
"""
function _solvematrixequations(P::ExtendedStateSpace, γ::Number)

    A = P.A
    B1 = P.B1
    B2 = P.B2
    C1 = P.C1
    C2 = P.C2
    D11 = P.D11
    D12 = P.D12
    D21 = P.D21
    D22 = P.D22

    P1 = size(C1, 1)
    P2 = size(C2, 1)
    M1 = size(B1, 2)
    M2 = size(B2, 2)

    γ² = γ * γ

    B = [B1 B2]
    C = [C1; C2]

    # Equation (7)
    D1dot = [D11 D12]
    R = [-γ²*I zeros(M1, M2); zeros(M2, M1) zeros(M2, M2)] + D1dot' * D1dot

    # Equation (8)
    Ddot1 = [D11; D21]
    Rbar = [-γ²*I zeros(P1, P2); zeros(P2, P1) zeros(P2, P2)] + Ddot1 * Ddot1'

    # Form hamiltonian for X and Y, equation (9) and (10)
    HX = [A zeros(size(A)); -C1'*C1 -A'] - ([B; -C1' * D1dot] / R) * [D1dot' * C1 B']
    HY = [A' zeros(size(A)); -B1*B1' -A] - ([C'; -B1 * Ddot1'] / Rbar) * [Ddot1 * B1' C]

    # Solve matrix equations
    Xinf = _solvehamiltonianare(HX)
    Yinf = _solvehamiltonianare(HY)

    # Equation (11)
    F = -R \ (D1dot' * C1 + B' * Xinf)

    # Equation (12)
    H = -(B1 * Ddot1' + Yinf * C') / Rbar

    return Xinf, Yinf, F, H
end

"""
    flag = _γiterations(A, B1, B2, C1, C2, D11, D12, D21, D22, maxIter, interval, verbose, tolerance)

Rune the complete set of γ-iterations over a specified search interval with a
set number of iterations. It is possible to break the algorithm if the number
of iterations become too large. This should perhaps be tken care of by
specifying an interval and tolerance for γ. In addition, the algorithm simply
terminates without a solution if the maximum possible γ on the defined
interval is infeasible. Here we could consider increasing the bounds somewhat
and warn the user if this occurrs.
"""
function _γiterations(
    P::ExtendedStateSpace,
    maxIter::Number,
    interval::Tuple,
    verbose::Bool,
    tolerance::Number,
)

    T = typeof(P.A)
    XinfFeasible, YinfFeasible, FinfFeasible, HinfFeasible, gammFeasible =
        T(undef,0,0), T(undef,0,0), T(undef,0,0), T(undef,0,0), T(undef,0,0)

    γ = maximum(interval)

    for iteration = 1:maxIter

        # Solve the matrix equations
        Xinf, Yinf, F, H = _solvematrixequations(P, γ)

        # Check Feasibility
        isFeasible =
            _checkfeasibility(Xinf, Yinf, γ, tolerance, iteration; verbose = verbose)

        if isFeasible
            XinfFeasible = Xinf
            YinfFeasible = Yinf
            FinfFeasible = F
            HinfFeasible = H
            gammFeasible = γ
            γ = γ - abs.(interval[2] - interval[1]) / (2^iteration)
        else
            γ = γ + abs.(interval[2] - interval[1]) / (2^iteration)
            if γ > maximum(interval)
                break
            end
        end
    end
    return XinfFeasible, YinfFeasible, FinfFeasible, HinfFeasible, gammFeasible
end

"""
    Pbar, Ltrans12, Rtrans12, Ltrans21, Rtrans21 = _transformP2Pbar(P::ExtendedStateSpace)

Transform the original system P to a new system Pbar, in which D12bar = [0; I]
and D21bar = [0 I] in order to satisfy the feasibility assumption A3 (see Doyle)
"""
function _transformp2pbar(P::ExtendedStateSpace)

    # Compute the transformation
    Ltrans12, Rtrans12 = _scalematrix(P.D12)
    Ltrans21, Rtrans21 = _scalematrix(P.D21)

    # Transform the system
    Abar = P.A
    B1bar = P.B1 * Rtrans21
    B2bar = P.B2 * Rtrans12
    C1bar = Ltrans12 * P.C1
    C2bar = Ltrans21 * P.C2
    D11bar = Ltrans12 * P.D11 * Rtrans21
    D12bar = Ltrans12 * P.D12 * Rtrans12
    D21bar = Ltrans21 * P.D21 * Rtrans21
    D22bar = Ltrans21 * P.D22 * Rtrans12
    Pbar = ss(Abar, B1bar, B2bar, C1bar, C2bar, D11bar, D12bar, D21bar, D22bar)

    return Pbar, Ltrans12, Rtrans12, Ltrans21, Rtrans21
end

"""
    Tl, Tr = _scalematrix(A::AbstractMatrix; method::String)

Find a left and right transform of A such that Tl*A*Tr = [I, 0], or
Tl*A*Tr = [I; 0], depending on the dimensionality of A.
"""
function _scalematrix(A::AbstractMatrix; method = "QR")
    # Check the rank condition
    if (minimum(size(A)) > 0)
        if rank(A) != minimum(size(A))
            error(ErrorException("Cannot scale the system, assumption A2 is violated"))
        end
    else
        error(
            ErrorException(
                "Cannot scale the system, minimum size of A must begreater than 0",
            ),
        )
    end

    # Perform scaling with the cosen method
    if method == "QR"
        return _coordinatetransformqr(A)
    elseif method == "SVD"
        return _coordinatetransformsvd(A)
    else
        error(
            ErrorException(
                string(
                    "The method",
                    method,
                    " is not supported, use 'QR' or 'SVD' instad.",
                ),
            ),
        )
    end
end

"""
    Tl, Tr =  _computeCoordinateTransformQR(A::AbstractMatrix)

Use the QR decomposition to find a transformaiton [Tl, Tr] such that
Tl*A*Tr becomes [0;I], [0 I] or I depending on the dimensionality of A.
"""
function _coordinatetransformqr(A::AbstractMatrix)
    m, n = size(A)
    if m == n
        # Square matrix with full rank
        Q, R = qr(A)
        LeftTransform = Q'
        RightTransform = inv(R)
    elseif m > n
        # Rectangular matrix with rank = n
        Q, R = qr(A)
        LeftTransform = [Q[:, (n+1):(m)]'; Q[:, 1:n]']
        RightTransform = inv(R[1:n, 1:n])
    elseif m < n
        # Rectangular matrix with rank = m
        Q, R = qr(A')
        LeftTransform = inv(R[1:m, 1:m])'
        RightTransform = [Q[:, m+1:n] Q[:, 1:m]]
    end
    return LeftTransform, RightTransform
end

"""
    Tl, Tr =  _computeCoordinateTransformSVD(A::AbstractMatrix)

Use the SVD to find a transformaiton [Tl, Tr] such that
Tl*A*Tr becomes [0;I], [0 I] or I depending on the dimensionality of A.
"""
function _coordinatetransformsvd(A::AbstractMatrix)
    m, n = size(A)
    if m == n
        # Square matrix with full rank
        U, S, V = svd(A)
        LeftTransform = inv(Diagonal(S)) * U'
        RightTransform = V
    elseif m > n
        # Rectangular matrix with rank = n
        U, S, V = svd(A; full = true)
        LeftTransform = [U[:, n+1:m]'; inv(Diagonal(S)) * U[:, 1:n]']
        RightTransform = V
    elseif m < n
        # Rectangular matrix with rank = m
        U, S, V = svd(A; full = true)
        LeftTransform = inv(Diagonal(S)) * U'
        RightTransform = [V[:, m+1:n] V[:, 1:m]]
    end
    return LeftTransform, RightTransform
end

"""
    P = hinfpartition(G, WS, WU, WT)

Transform a SISO or MIMO system G, with weighting functions WS, WU, WT into
and LFT with an isolated controller, and write the resulting system, P(s),
on a state-space form. Valid inputs for G are transfer functions (with dynamics,
can be both MIMO and SISO, both in tf and ss forms). Valid inputs for the
weighting functions are empty arrays, numbers (static gains), and `LTISystem`s.
"""
function hinfpartition(G, WS, WU, WT)
    # # Convert the systems into state-space form
    Ag, Bg, Cg, Dg = _input2ss(G)
    Asw, Bsw, Csw, Dsw = _input2ss(WS)
    Auw, Buw, Cuw, Duw = _input2ss(WU)
    Atw, Btw, Ctw, Dtw = _input2ss(WT)

    G = ss(Ag, Bg, Cg, Dg)
    WS = ss(Asw, Bsw, Csw, Dsw)
    WU = ss(Auw, Buw, Cuw, Duw)
    WT = ss(Atw, Btw, Ctw, Dtw)

    # Check that the system is realizable
    if size(Cg, 1) != size(Btw, 2) && size(Btw, 2) != 0
        if ControlSystems.issiso(WT)
            WT = ControlSystems.append(fill(WT, G.ny)...)
            return hinfpartition(G, WS, WU, WT)
        end
        println([size(Cg, 1), size(Btw, 2)])
        throw(DimensionMismatch(
            "You must have the same number of inputs to WT as there are outputs",
        ))
    end
    if size(Cg, 1) != size(Bsw, 2) && size(Bsw, 2) != 0
        if ControlSystems.issiso(WS)
            WS = ControlSystems.append(fill(WS, G.ny)...)
            return hinfpartition(G, WS, WU, WT)
        end
        println([size(Cg, 1), size(Bsw, 2)])
        throw(DimensionMismatch(
            "You must have the same number of inputs to WS as there are outputs",
        ))
    end
    if size(Bg, 2) != size(Buw, 2) && size(Buw, 2) != 0
        if ControlSystems.issiso(WU)
            WU = ControlSystems.append(fill(WU, G.nu)...)
            return hinfpartition(G, WS, WU, WT)
        end
        println([size(Bg, 2), size(Buw, 2)])
        throw(DimensionMismatch(
            "You must have the same number of inputs to WU as there are controls",
        ))
    end
    if (
        size(Ag, 1) == 0 ||
        size(Ag, 2) == 0 ||
        size(Bg, 1) == 0 ||
        size(Bg, 2) == 0 ||
        size(Cg, 1) == 0 ||
        size(Cg, 2) == 0 ||
        size(Dg, 1) == 0 ||
        size(Dg, 2) == 0
    )
        throw(DimensionMismatch(
            "Expansion of systems dimensionless A,B,C or D is not yet supported",
        ))
    end

    # Form A
    (mAg, nAg) = size(Ag)
    (mAsw, nAsw) = size(Asw)
    (mAuw, nAuw) = size(Auw)
    (mAtw, nAtw) = size(Atw)

    if size(Bsw, 1) != 0 && size(Bsw, 2) != 0
        BswCg = Bsw * Cg
    else
        BswCg = zeros(mAsw, nAg)
    end
    if size(Btw, 1) != 0 && size(Btw, 2) != 0
        BtwCg = Btw * Cg
    else
        BtwCg = zeros(mAtw, nAg)
    end

    @debug ([(mAg, nAg), (mAsw, nAsw), (mAuw, nAuw), (mAtw, nAtw)])
    A = [
        Ag zeros(mAg, nAsw) zeros(mAg, nAuw) zeros(mAg, nAtw)
        -BswCg Asw zeros(mAsw, nAuw) zeros(mAsw, nAtw)
        zeros(mAuw, nAg) zeros(mAuw, nAsw) Auw zeros(mAuw, nAtw)
        BtwCg zeros(mAtw, nAsw) zeros(mAtw, nAuw) Atw
    ]

    if size(Buw, 2) == 0
        Buw = zeros(0, size(Bg, 2))
    end

    (mBg, nBg) = size(Bg)
    (mBsw, nBsw) = size(Bsw)
    (mBuw, nBuw) = size(Buw)
    (mBtw, nBtw) = size(Btw)

    Bw = [zeros(mBg, nBsw); Bsw; zeros(mBuw, nBsw); zeros(mAtw, nBsw)]
    Bu = [Bg; zeros(mBsw, nBg); Buw; zeros(mAtw, nBg)]

    (mCg, nCg) = size(Cg)
    (mCsw, nCsw) = size(Csw)
    (mCuw, nCuw) = size(Cuw)
    (mCtw, nCtw) = size(Ctw)

    if size(Dsw, 1) != 0 && size(Dsw, 1) != 0
        DswCg = Dsw * Cg
    else
        DswCg = zeros(0, nAg)
    end
    if size(Dtw, 1) != 0 && size(Dtw, 1) != 0
        DtwCg = Dtw * Cg
    else
        DtwCg = zeros(0, nAg)
    end

    Cz = [
        -DswCg Csw zeros(mCsw, nAuw) zeros(mCsw, nAtw)
        zeros(mCuw, nAg) zeros(mCuw, nAsw) Cuw zeros(mCuw, nAtw)
        DtwCg zeros(mCtw, nAsw) zeros(mCtw, nAuw) Ctw
    ]
    Cy = [-Cg zeros(mCg, nAsw) zeros(mCg, nAuw) zeros(mCg, nAtw)]

    if size(Duw, 2) == 0
        Duw = zeros(0, size(Bg, 2))
    end

    (mDg, nDg) = size(Dg)
    (mDsw, nDsw) = size(Dsw)
    (mDuw, nDuw) = size(Duw)
    (mDtw, nDtw) = size(Dtw)

    Dzw = [Dsw; zeros(mDuw, nDsw); zeros(mDtw, nDsw)]
    Dzu = [zeros(mDsw, nDuw); Duw; zeros(mDtw, nDuw)]
    Dyw = Matrix{Float64}(I, mCg, nDuw)
    Dyu = -Dg

    P = ss(A, Bw, Bu, Cz, Cy, Dzw, Dzu, Dyw, Dyu)

end

function hinfpartition2(G, WS, WU, WT)
    # # Convert the systems into state-space form
    Ag, Bg, Cg, Dg = _input2ss(G)
    Asw, Bsw, Csw, Dsw = _input2ss(WS)
    Auw, Buw, Cuw, Duw = _input2ss(WU)
    Atw, Btw, Ctw, Dtw = _input2ss(WT)

    G = ss(Ag, Bg, Cg, Dg)
    WS = ss(Asw, Bsw, Csw, Dsw)
    WU = ss(Auw, Buw, Cuw, Duw)
    WT = ss(Atw, Btw, Ctw, Dtw)

    # Check that the system is realizable
    if size(Cg, 1) != size(Btw, 2) && size(Btw, 2) != 0
        if ControlSystems.issiso(WT)
            WT = ControlSystems.append(fill(WT, G.ny)...)
            return hinfpartition(G, WS, WU, WT)
        end
        println([size(Cg, 1), size(Btw, 2)])
        throw(DimensionMismatch(
            "You must have the same number of inputs to WT as there are outputs",
        ))
    end
    if size(Cg, 1) != size(Bsw, 2) && size(Bsw, 2) != 0
        if ControlSystems.issiso(WS)
            WS = ControlSystems.append(fill(WS, G.ny)...)
            return hinfpartition(G, WS, WU, WT)
        end
        println([size(Cg, 1), size(Bsw, 2)])
        throw(DimensionMismatch(
            "You must have the same number of inputs to WS as there are outputs",
        ))
    end
    if size(Bg, 2) != size(Buw, 2) && size(Buw, 2) != 0
        if ControlSystems.issiso(WU)
            WU = ControlSystems.append(fill(WU, G.nu)...)
            return hinfpartition(G, WS, WU, WT)
        end
        println([size(Bg, 2), size(Buw, 2)])
        throw(DimensionMismatch(
            "You must have the same number of inputs to WU as there are controls",
        ))
    end
    if (
        size(Ag, 1) == 0 ||
        size(Ag, 2) == 0 ||
        size(Bg, 1) == 0 ||
        size(Bg, 2) == 0 ||
        size(Cg, 1) == 0 ||
        size(Cg, 2) == 0 ||
        size(Dg, 1) == 0 ||
        size(Dg, 2) == 0
    )
        throw(DimensionMismatch(
            "Expansion of systems dimensionless A,B,C or D is not yet supported",
        ))
    end

    ny,nu = G.ny, G.nu
    Iny = I(ny)
    perm1 = [Iny; zeros(ny+nu,ny); Iny]
    perm2 = [zeros(ny,nu) ; I(nu) ; zeros(2*ny,nu)] + [-Iny ; zeros(nu,ny) ; Iny ; -Iny] * G
    W = ControlSystems.blockdiag(WS, WU, WT, ss(Iny))
    P = W * [perm1 perm2]
    nz = P.ny - ny
    nw = P.nu - nu
    zi = 1:nz
    yi = nz+1:P.ny
    wi = 1:nw
    ui = nw+1:P.nu
    ss(P.A, P.B[:, wi], P.B[:, ui], P.C[zi, :], P.C[yi, :],
        P.D[zi, wi], P.D[zi, ui], P.D[yi, wi], P.D[yi, ui], G.timeevol)
end

function ControlSystems.blockdiag(systems::AbstractStateSpace...)
    ST = promote_type(typeof.(systems)...)
    timeevol = common_timeevol(systems...)
    A = ControlSystems.blockdiag([s.A for s in systems]...)
    B = ControlSystems.blockdiag([_remove_empty_cols(s.B) for s in systems]...)
    C = ControlSystems.blockdiag([s.C for s in systems]...)
    D = ControlSystems.blockdiag([_remove_empty_cols(s.D) for s in systems]...)
    return ST(A, B, C, D, timeevol)
end

_remove_empty_cols(x) = size(x,2) == 0 ? zeros(size(x,1), 1) : x

ControlSystems.numeric_type(n::Number) = typeof(n)
ControlSystems.numeric_type(::Any) = Float64
ControlSystems.numeric_type(::AbstractVector{T}) where T = numeric_type(T)

"""
    convert_input_to_ss(H)

Helper function used for type conversion in hinfpartition()
"""
function _input2ss(H)
    T = ControlSystems.numeric_type(H) |> float
    if isa(H, LTISystem)
        if isa(H, TransferFunction)
            Hss = ss(H)
        else
            Hss = H
        end
        Ah, Bh, Ch, Dh = Hss.A, Hss.B, Hss.C, Hss.D
    elseif isa(H, Number)
        Ah, Bh, Ch, Dh = zeros(T, 0, 0), zeros(T, 0, 1), zeros(T, 1, 0), H * ones(T, 1, 1)
    else # e.g. H = nothing
        Ah, Bh, Ch, Dh = zeros(0, 0), zeros(0, 0), zeros(0, 0), zeros(0, 0)
    end
    return Ah, Bh, Ch, Dh
end

"""
    hinfsignals(P::ExtendedStateSpace, G::LTISystem, C::LTISystem)

Use the extended state-space model, a plant and the found controller to extract
the closed loop transfer functions operating solely on the state-space.

- `Pcl : w → z` : From input to the weighted functions
- `S   : w → e` : From input to error
- `CS  : w → u` : From input to control
- `T   : w → y` : From input to output
"""
function hinfsignals(P::ExtendedStateSpace, G::LTISystem, C::LTISystem)

    common_timeevol(P,G,C)
    A, B1, B2, C1, C2, D11, D12, D21, D22 = ssdata_e(P)
    Ag, Bg, Cg, Dg = ssdata(ss(G))
    Ac, Bc, Cc, Dc = ssdata(ss(C))

    # Precompute the inverse
    M = inv(I - Dc * D22)
    A11 = A + B2 * M * Dc * C2
    A12 = B2 * M * Cc
    A21 = Bc * C2 + Bc * D22 * M * Dc * C2
    A22 = Ac + Bc * D22 * M * Cc
    B11 = B1 + B2 * M * Dc * D21
    B21 = Bc * D21 + Bc * D22 * M * Dc * D21

    C_w2z_11 = C1 + D12 * M * Dc * C2
    C_w2z_12 = D12 * M * Cc
    D_w2z = D11 + D12 * M * Dc * D21

    Pw2z = ss([A11 A12; A21 A22], [B11; B21], [C_w2z_11 C_w2z_12], D_w2z, timeevol(P))

    C_w2e_11 = C2 + D22 * M * Dc * C2
    C_w2e_12 = D22 * M * Cc
    D_w2e = D21 + D22 * M * Dc * D21

    Pw2e = ss([A11 A12; A21 A22], [B11; B21], [C_w2e_11 C_w2e_12], D_w2e, timeevol(P))

    C_w2u_11 = M * Dc * C2
    C_w2u_12 = M * Cc
    D_w2u = M * Dc * D21

    Pw2u = ss([A11 A12; A21 A22], [B11; B21], [C_w2u_11 C_w2u_12], D_w2u, timeevol(P))

    Abar = [A11 A12; A21 A22]
    Bbar = [B11; B21]
    C_w2y_11 = M * Dc * C2
    C_w2y_12 = M * Cc
    D_w2y = M * Dc * D21
    Cbar1 = [Cg zeros(size(Cg, 1), size(Abar, 2) - size(Cg, 2))]
    Cbar2 = [Dg * C_w2y_11 Dg * C_w2y_12]
    Dbar = Dg * M * Dc * D21

    Pw2y = ss(Abar, Bbar, Cbar1 + Cbar2, Dbar, timeevol(P))

    Pcl = Pw2z
    S = Pw2e
    CS = Pw2u
    T = Pw2y

    return Pcl, S, CS, T
end

"""
    bilineard2c(Ad::AbstractArray, Bd::AbstractArray, Cd::AbstractArray, Dd::AbstractArray, Ts::Number; tolerance=1e-12)

Balanced Bilinear transformation in State-Space. This method computes a
continuous time equivalent of a discrete time system, such that

    G_c(z) = z2s[G_d(z)]

in a manner which accomplishes the following
  (i)   Preserves the infinity L-infinity norm over the transformation
  (ii)  Finds a system which balances B and C, in the sense that ||B||_2=||C||_2
  (iii) Satisfies G_d(z) = s2z[z2s[G_d(z)]] for some map s2z[]
"""
function bilineard2c(
    Ad::AbstractArray,
    Bd::AbstractArray,
    Cd::AbstractArray,
    Dd::AbstractArray,
    Ts::Number;
    tolerance = 1e-12,
)

    Id = Matrix{Float64}(I, size(Ad, 1), size(Ad, 2))

    Pd = Ad - Id
    Qd = Ad + Id
    ialpha = 2 / Ts #Should be this, but the nyquist frequency doesnt add up unless
    ialpha = 1 / Ts

    Ac = ialpha * (Pd / Qd)
    Bc = Qd \ Bd
    Cc = 2 * ialpha * (Cd / Qd)
    Dc = Dd - (Cd / Qd) * Bd

    # Scaling for improved numerical stability
    σB = maximum(svd(Bc).S)
    σC = maximum(svd(Cc).S)
    if σB > tolerance && σC > tolerance
        λd = sqrt(σB / σC)
    else
        λd = 1
        error(
            "Warning, the problem is poorly cnditioned. Consider an alternate discretization scheme.",
        )
    end
    Bc /= λd
    Cc *= λd

    return Ac, Bc, Cc, Dc
end

"""
    bilineard2c(sys::StateSpace)

Applies a Balanced Bilinear transformation to continuous-time statespace object
"""
function bilineard2c(sys::StateSpace)
    Ad, Bd, Cd, Dd = ssdata(sys)
    Ts = sys.Ts

    if Ts <= 0
        error("Error, the input must be a discrete time system.")
    end

    Ac, Bc, Cc, Dc = bilineard2c(Ad, Bd, Cd, Dd, Ts)
    return ss(Ac, Bc, Cc, Dc)
end

"""
    bilineard2c(sys::ExtendedStateSpace)

Applies a Balanced Bilinear transformation to continuous-time extended statespace object
"""
function bilineard2c(sys::ExtendedStateSpace)
    Ad = sys.A
    Bd = sys.B
    Cd = sys.C
    Dd = sys.D
    Ts = sys.Ts

    m1 = size(sys.B1, 2)
    m2 = size(sys.B2, 2)
    p1 = size(sys.C1, 1)
    p2 = size(sys.C2, 1)

    if Ts <= 0
        error("Error, the input must be a discrete time system.")
    end

    Ac, Bc, Cc, Dc = bilineard2c(Ad, Bd, Cd, Dd, Ts)

    A = Ac
    B1 = Bc[:, 1:m1]
    B2 = Bc[:, (m1+1):(m1+m2)]
    C1 = Cc[1:p1, :]
    C2 = Cc[(p1+1):(p1+p2), :]
    D11 = Dc[1:p1, 1:m1]
    D12 = Dc[1:p1, (m1+1):(m1+m2)]
    D21 = Dc[(p1+1):(p1+p2), 1:m1]
    D22 = Dc[(p1+1):(p1+p2), (m1+1):(m1+m2)]

    return ss(A, B1, B2, C1, C2, D11, D12, D21, D22)
end

"""
    bilinearc2d(Ac::AbstractArray, Bc::AbstractArray, Cc::AbstractArray, Dc::AbstractArray, Ts::Number; tolerance=1e-12)

Balanced Bilinear transformation in State-Space. This method computes a
discrete time equivalent of a continuous-time system, such that

    G_d(z) = s2z[G_c(s)]

in a manner which accomplishes the following
  (i)   Preserves the infinity L-infinity norm over the transformation
  (ii)  Finds a system which balances B and C, in the sense that ||B||_2=||C||_2
  (iii) Satisfies G_c(s) = z2s[s2z[G_c(s)]] for some map z2s[]
"""
function bilinearc2d(
    Ac::AbstractArray,
    Bc::AbstractArray,
    Cc::AbstractArray,
    Dc::AbstractArray,
    Ts::Number;
    tolerance = 1e-12,
)

    Id = Matrix{Float64}(I, size(Ac, 1), size(Ac, 2))
    alpha = Ts / 2 #Should be this, but the nyquist frequency doesnt add up
    alpha = Ts

    # Check that the bilinear tranformation is possible
    if minimum(svd(Id - alpha * Ac).S) < 1e-12
        error(
            "The transformation is extremely poorly conditioned, with min(svd(Id - alpha * Ac).S) < 1e-12. Consider an alternate discretization scheme.",
        )
    end

    PP = Id - alpha * Ac
    QQ = Id + alpha * Ac

    Ad = PP \ QQ
    Bd = (PP \ Bc)
    Cd = 2 * alpha * (Cc / PP)

    # Scaling for improved numerical stability
    σB = maximum(svd(Bd).S)
    σC = maximum(svd(Cd).S)
    if σB > tolerance && σC > tolerance
        λc = sqrt(σB / σC)
    else
        λc = 1
        error(
            "Warning, the problem is poorly cnditioned. Consider an alternate discretization scheme.",
        )
    end

    Bd /= λc
    Cd *= λc

    Dd = alpha * Cc / PP * Bc + Dc
    return Ad, Bd, Cd, Dd, Ts
end

"""
    bilinearc2d(sys::StateSpace, Ts::Number)

Applies a Balanced Bilinear transformation to a discrete-time statespace object
"""
function bilinearc2d(sys::StateSpace{Continuous}, Ts::Number)
    Ac, Bc, Cc, Dc = ssdata(sys)

    if Ts <= 0
        throw(ArgumentError("Error, the the discretization time Ts must be positive."))
    end

    Ad, Bd, Cd, Dd = bilinearc2d(Ac, Bc, Cc, Dc, Ts)
    return ss(Ad, Bd, Cd, Dd, Ts)
end

"""
    bilinearc2d(sys::ExtendedStateSpace, Ts::Number)

Applies a Balanced Bilinear transformation to a discrete-time extended statespace object
"""
function bilinearc2d(sys::ExtendedStateSpace{Continuous}, Ts::Number)
    Ac = sys.A
    Bc = sys.B
    Cc = sys.C
    Dc = sys.D

    m1 = size(sys.B1, 2)
    m2 = size(sys.B2, 2)
    p1 = size(sys.C1, 1)
    p2 = size(sys.C2, 1)

    if Ts <= 0
        error("Error, the the discretization time Ts must be positive.")
    end

    Ad, Bd, Cd, Dd = bilinearc2d(Ac, Bc, Cc, Dc, Ts)

    A = Ad
    B1 = Bd[:, 1:m1]
    B2 = Bd[:, (m1+1):(m1+m2)]
    C1 = Cd[1:p1, :]
    C2 = Cd[(p1+1):(p1+p2), :]
    D11 = Dd[1:p1, 1:m1]
    D12 = Dd[1:p1, (m1+1):(m1+m2)]
    D21 = Dd[(p1+1):(p1+p2), 1:m1]
    D22 = Dd[(p1+1):(p1+p2), (m1+1):(m1+m2)]

    return ss(A, B1, B2, C1, C2, D11, D12, D21, D22, Ts)
end

function fudge_inv(s::AbstractStateSpace, ε = 1e-3)
    s = deepcopy(s)
    s.D .+= ε*I(size(s.D,1))
    inv(s)
end


text:
using RatPack
using Distributions
using DataFrames
using CSV
using Test

# test case from "Whos's #1", Langville and Meyer
file = "../data/test_competitions_1.csv"
(input_competitions,  player_list) = read_results( file )
summarize(input_competitions,  player_list)
input_ratings = RatingsList( player_list )
r = update_ratings(UpdateMassey(),
                   input_ratings,
                   input_competitions)

R1 = RatingsTable(input_ratings.players) # empty
R2 = RatingsTable(input_ratings)         # 1 row but now values
R3 = RatingsTable(r)                     # 1 initialised row
R4 = RatingsTable(r.players, 4)          # 4 uninitialised rows
setindex!(R4, r, 2)
R4[3] = r
@testset "RatingsTable structure" begin
    @test append!(R1,R3) == R3
    @test R1.ratings[1,:] == R1[1]
    @test R4[2:2] == RatingsTable(r).ratings
end

massey_ratings = update_ratings(UpdateMassey(),
                                input_ratings,
                                input_competitions)

# include home ground factor in the original case from  "Whos's #1", Langville and Meyer
file2 = "../data/test_competitions_2.csv"
(input_competitions2,  player_list2) = read_results( file2 )
summarize(input_competitions2,  player_list2)
input_ratings2 = RatingsList( player_list2 )

# test case from NFL season 2009 (see file for source, but also used in "Whos's #1", Langville and Meyer)
#   this is the regular season, excluding playoffs
file = "../data/nfl_2009_regular.csv"
(nfl_competitions,  nfl_player_list) = read_results( file )
nfl_ratings = RatingsList( nfl_player_list )

#   this is the whole season, including playoffs
file = "../data/nfl_2009.csv"
(nfl_ext_competitions,  nfl_ext_player_list) = read_results( file )
nfl_ext_ratings = RatingsList( nfl_ext_player_list )

# small netflix-like example from "Whos's #1", Langville and Meyer
file = "../data/small_netflix_eg.csv"
(netflix_competitions,  netflix_player_list) = read_results( file )
netflix_ratings = RatingsList( netflix_player_list )

@testset "Utilities" begin
    @testset "Line search" begin
        r = [0.0, 2.0]
        ε = 1.0e-6
        (x_star1, y_star1, error_estimate1, n_evaluations2) =
            linesearch( "Grid", x->sin(x); range=r, ε=ε)
        (x_star2, y_star2, error_estimate2, n_evaluations2) =
            linesearch( "Golden", x->sin(x); range=r, ε=ε)
        @test abs(x_star1 - x_star2) < ε
    end
end

@testset "I/O" begin
    @test_throws ArgumentError o = read_results( "test" )
    # @test_throws ErrorException SurrealFinite("1", [x1], [x0])
    # @test_throws DomainError convert(SurrealFinite, NaN )

    @test (sort(summarize(input_competitions,  player_list), :Player) ==
           DataFrame(Player=["Duke","Miami","UNC","UVA","VT"], Record=["0/4","4/4","2/4","1/4","3/4"], ScoreDiff=[-124,91,-40,-17,90])
           )

    test_out = "Tmp/test_ratings.csv"
    write_ratings( test_out,  massey_ratings)
    df,r = read_ratings( test_out )
    @test r == massey_ratings

    # test for a wide range of input errors in competitions
    error_list = ["Input DataFrame must have a 'Player A' column'",
                  "Input DataFrame 'Player A' column must have type String",
                  "Input DataFrame must have a 'Player B' column'",
                  "Input DataFrame 'Player B' column must have type String",
                  "Input DataFrame 'Outcome' column must have type Union{Missing, Integer}",
                  "Input DataFrame 'Outcome' column should be -1,0,1 (or missing)",
                  "Input DataFrame 'Result' column must have type Union{Missing, Real}",
                  "Input DataFrame 'Result' column should be >= 0 (or missing)",
                  "Input DataFrame 'No. of matches' column must have type Union{Missing, Integer}",
                  "Input DataFrame 'No. of matches' column should be >= 1 (or missing)",
                  "Input DataFrame 'Score A' column must have type Union{Missing, Integer}",
                  "Input DataFrame 'Score B' column must have type Union{Missing, Integer}",
                  "Input DataFrame only has 'Score A' column, not 'Score B'",
                  "Input DataFrame only has 'Score B' column, not 'Score A'",
                  "Input DataFrame 'Margin' column must have type Union{Missing, Integer}",
                  "Input DataFrame 'Factor A' column must have type Union{Missing, Integer}",
                  "Input DataFrame 'Factor B' column must have type Union{Missing, Integer}",
                  "Input DataFrame only has 'Factor A' column, not 'Factor B'",
                  "Input DataFrame only has 'Factor B' column, not 'Factor A'",
                  "Input DataFrame: Margins are inconsistent with Scores",
                  "Input DataFrame: Outcomes are inconsistent with Margins",
                  "",
                  "",
                  ]

    errored_file_dir = "../data/errored_files"
    files = readdir(errored_file_dir)
    i = 1
    for f in files
        # global file, df, i
        if match(r".csv$", f) !== nothing
            file = "$errored_file_dir/$f"
            # println("testing $file ($i of $(length(files)))")
            df = CSV.read(file; comment="#")
            @test_throws ErrorException(error_list[i]) check_results( df )

            i += 1
        end
    end

end

@testset "Info" begin
    for (i,u) in enumerate(update_rule_list)
        if u != "Iterate" && u != "SampleIterate"
            nm = Meta.parse("Update$u()")
            info_update_rule = :(update_info( $nm ) )
            D = eval( info_update_rule )
            # println("$u, $(D[:name])")
            @test D[:name] == u
        end
    end
end

@testset "Massey" begin
    # see \"Whos's #1\", Langville and Meyer, p.11
    required_output = Dict{String,Float64}(
        "Duke"  => -24.8,
        "Miami" =>  18.2,
        "UNC"   =>  -8.0,
        "UVA"   =>  -3.4,
        "VT"    =>  18.0)
    for k in keys(required_output)
        @test abs(required_output[k] - massey_ratings.ratings[k]) < 1.0e-1
    end

    massey_ratings2 = update_ratings(UpdateMassey(),
                                     netflix_ratings,
                                     netflix_competitions)
    # see \"Whos's #1\", Langville and Meyer, Table 3.3, p.26
    required_output = Dict{String,Float64}(
        "M1" =>  0.65,
        "M2" =>  1.01,
        "M3" => -0.55,
        "M4" => -1.11
        )
    for k in keys(required_output)
        @test abs(required_output[k] - massey_ratings2.ratings[k]) < 1.0e-2
    end

end

@testset "Colley" begin
    colley_ratings = update_ratings(UpdateColley(),
                                    input_ratings,
                                    input_competitions)
    # see \"Whos's #1\", Langville and Meyer, Tabl 3.1, p.23
    required_output = Dict{String,Float64}(
        "Duke"  => 0.21,
        "Miami" => 0.79,
        "UNC"   => 0.50,
        "UVA"   => 0.36,
        "VT"    => 0.64)
    for k in keys(required_output)
        @test abs(required_output[k] - colley_ratings.ratings[k]) < 1.0e-2
    end

    colley_ratings2 = update_ratings(UpdateColley(), netflix_ratings,
                                     netflix_competitions)
    # see \"Whos's #1\", Langville and Meyer, Table 3.3, p.26
    required_output = Dict{String,Float64}(
        "M1" =>  0.67,
        "M2" =>  0.63,
        "M3" =>  0.34,
        "M4" =>  0.35
        )
    for k in keys(required_output)
        @test abs(required_output[k] - colley_ratings2.ratings[k]) < 1.0e-2
    end

end

@testset "MasseyColley" begin
    massey_colley_ratings = update_ratings(UpdateMasseyColley(),
                                           input_ratings,
                                           input_competitions)
    # see \"Whos's #1\", Langville and Meyer, Tabl 3.2, p.25
    required_output = Dict{String,Float64}(
        "Duke"  => -17.7,
        "Miami" =>  13.0,
        "UNC"   =>  -5.7,
        "UVA"   =>  -2.4,
        "VT"    =>  12.9)
    for k in keys(required_output)
        @test abs(required_output[k] - massey_colley_ratings.ratings[k]) < 1.0e-1
    end

end

@testset "KeenerScores" begin
    keener_score_ratings = update_ratings(UpdateKeenerScores(; skew=true, norm=false),
                                          nfl_ratings,
                                          nfl_competitions)
    # see \"Whos's #1\", Langville and Meyer, p.44
    required_output = Dict{String,Float64}(
   "New Orleans Saints" => 0.0361385,
    "Green Bay Packers" => 0.0357225,
 "New England Patriots" => 0.0350508,
   "San Diego Chargers" => 0.0350264,
   "Indianapolis Colts" => 0.0348175,
    "Minnesota Vikings" => 0.0347833,
       "Dallas Cowboys" => 0.0347101,
        "New York Jets" => 0.0346833,
  "Philadelphia Eagles" => 0.0338832,
     "Baltimore Ravens" => 0.0338213,
  "Pittsburgh Steelers" => 0.0335293,
       "Houston Texans" => 0.0334149,
      "Atlanta Falcons" => 0.0326898,
    "Arizona Cardinals" => 0.0323456,
  "San Francisco 49ers" => 0.0318757,
       "Denver Broncos" => 0.0317885,
   "Cincinnati Bengals" => 0.0314828,
    "Carolina Panthers" => 0.0307855,
     "Tennessee Titans" => 0.0305376,
      "New York Giants" => 0.0304804,
       "Miami Dolphins" => 0.0298047,
        "Chicago Bears" => 0.0294099,
  "Washington Redskins" => 0.0291071,
        "Buffalo Bills" => 0.0290657,
 "Jacksonville Jaguars" => 0.0289624,
   "Kansas City Chiefs" => 0.0280055,
     "Cleveland Browns" => 0.0279232,
     "Seattle Seahawks" => 0.0272619,
      "Oakland Raiders" => 0.0262216,
 "Tampa Bay Buccaneers" => 0.0261939,
        "Detroit Lions" => 0.0255954,
       "St. Louis Rams" => 0.0248815
    )
    for k in keys(required_output)
        @test abs(required_output[k] - keener_score_ratings.ratings[k]) < 1.0e-5
    end
end

@testset "Revert" begin
    colley_ratings = update_ratings(UpdateColley(),
                                    input_ratings,
                                    input_competitions)
    r0 = 1.0
    a = 0.25
    revert_ratings = update_ratings(UpdateRevert(; r_base=r0, α=a),
                                    colley_ratings)
    required_output = Dict{String,Float64}(
        "Duke"  => 0.21*a + r0*(1-a),
        "Miami" => 0.79*a + r0*(1-a),
        "UNC"   => 0.50*a + r0*(1-a),
        "UVA"   => 0.36*a + r0*(1-a),
        "VT"    => 0.64*a + r0*(1-a))
    for k in keys(required_output)
        @test abs(required_output[k] - revert_ratings.ratings[k]) < 1.0e-2
    end

    # check it complains when you give invalid arguments
    @test_throws ArgumentError rule = UpdateRevert(α = -0.1)
    @test_throws ArgumentError rule = UpdateRevert(α =  1.1)
end

@testset "Elo" begin
    r0 = 0.0
    rule = UpdateElo(;r0 = r0, K = 32.0, θ=1000.0/log(10.0) )
    r = nfl_ext_ratings
    for i=1:size(nfl_ext_competitions,1)
        # apply Elo 1-by-1 to each result
         r = update_ratings(rule, r, nfl_ext_competitions[ i:i, :])
    end
    S = sort(collect(r.ratings), by = tuple -> last(tuple), rev=true)
    r_mean = mean(last.(S)) # should be close to r0 = 0.0
    @test abs(r0 - r_mean) < 1.0e-6
    elo_ratings = r
    # see \"Whos's #1\", Langville and Meyer, p.58
    required_output = Dict{String,Float64}(
   "New Orleans Saints" => 173.661,
   "Indianapolis Colts" => 170.331,
   "San Diego Chargers" => 127.582,
    "Minnesota Vikings" => 103.504,
       "Dallas Cowboys" => 89.1285,
  "Philadelphia Eagles" => 69.5331,
    "Green Bay Packers" => 67.8292,
    "Arizona Cardinals" => 53.227 ,
        "New York Jets" => 50.1431,
 "New England Patriots" => 39.6328,
       "Houston Texans" => 33.9024,
   "Cincinnati Bengals" => 33.0116,
     "Baltimore Ravens" => 32.0826,
      "Atlanta Falcons" => 28.1178,
  "Pittsburgh Steelers" => 27.1246,
     "Tennessee Titans" => 13.2216,
    "Carolina Panthers" => 11.4745,
  "San Francisco 49ers" => -1.28445,
      "New York Giants" => -5.3217,
       "Denver Broncos" => -11.1262,
       "Miami Dolphins" => -26.7173,
        "Chicago Bears" => -28.1416,
 "Jacksonville Jaguars" => -36.2142,
        "Buffalo Bills" => -53.3495,
     "Cleveland Browns" => -74.6639,
      "Oakland Raiders" => -83.3188,
     "Seattle Seahawks" => -88.8452,
   "Kansas City Chiefs" => -109.281,
  "Washington Redskins" => -110.212,
 "Tampa Bay Buccaneers" => -130.102,
        "Detroit Lions" => -170.809,
       "St. Louis Rams" => -194.119
                                           )
    for k in keys(required_output)
        @test abs(required_output[k] - elo_ratings.ratings[k]) < 1.0e-2
    end

    # do a test for: K=32 weeks 1-15; K=16 weeks 16-17; K=64 for playoffs
    #  p.60-61

    # check it complains when you give invalid arguments
    @test_throws ArgumentError rule = UpdateElo(K = -1)
    @test_throws ArgumentError rule = UpdateElo(θ = -1)
end

@testset "EloF" begin
    r0 = 0.0
    rule = UpdateEloF(;r0 = r0, K = 32.0, θ=1000.0/log(10.0), factor_scale = 0.0 )
    r = nfl_ext_ratings
    for i=1:size(nfl_ext_competitions,1)
        # apply Elo 1-by-1 to each result
        r = update_ratings(rule, r, nfl_ext_competitions[ i:i, :])
    end
    S = sort(collect(r.ratings), by = tuple -> last(tuple), rev=true)
    r_mean = mean(last.(S)) # should be close to r0 = 0.0
    @test abs(r0 - r_mean) < 1.0e-6
    elo_ratings = r
    # see \"Whos's #1\", Langville and Meyer, p.58
    required_output = Dict{String,Float64}(
   "New Orleans Saints" => 173.661,
   "Indianapolis Colts" => 170.331,
   "San Diego Chargers" => 127.582,
    "Minnesota Vikings" => 103.504,
       "Dallas Cowboys" => 89.1285,
  "Philadelphia Eagles" => 69.5331,
    "Green Bay Packers" => 67.8292,
    "Arizona Cardinals" => 53.227 ,
        "New York Jets" => 50.1431,
 "New England Patriots" => 39.6328,
       "Houston Texans" => 33.9024,
   "Cincinnati Bengals" => 33.0116,
     "Baltimore Ravens" => 32.0826,
      "Atlanta Falcons" => 28.1178,
  "Pittsburgh Steelers" => 27.1246,
     "Tennessee Titans" => 13.2216,
    "Carolina Panthers" => 11.4745,
  "San Francisco 49ers" => -1.28445,
      "New York Giants" => -5.3217,
       "Denver Broncos" => -11.1262,
       "Miami Dolphins" => -26.7173,
        "Chicago Bears" => -28.1416,
 "Jacksonville Jaguars" => -36.2142,
        "Buffalo Bills" => -53.3495,
     "Cleveland Browns" => -74.6639,
      "Oakland Raiders" => -83.3188,
     "Seattle Seahawks" => -88.8452,
   "Kansas City Chiefs" => -109.281,
  "Washington Redskins" => -110.212,
 "Tampa Bay Buccaneers" => -130.102,
        "Detroit Lions" => -170.809,
       "St. Louis Rams" => -194.119
                                           )
    for k in keys(required_output)
        @test abs(required_output[k] - elo_ratings.ratings[k]) < 1.0e-2
    end

    # do a test where factors, e.g., home ground are important
    rule = UpdateEloF(;r0 = r0, K = 32.0, θ=1000.0/log(10.0), factor_scale=10.0 )
    r = input_ratings2
    for i=1:10
        r = update_ratings(rule,
                           r,
                           input_competitions2)
    end
    # check the output against something???

    # check it complains when you give invalid arguments
    @test_throws ArgumentError rule = UpdateEloF(K = -1)
    @test_throws ArgumentError rule = UpdateEloF(θ = -1)
    @test_throws ArgumentError rule = UpdateEloF(factor_scale = -1.0)
end

@testset "Iterate" begin
    r0 = 0.0
    # nice way to say, do Elo line-by-line on the inputs
    rule = UpdateIterate( UpdateElo(;r0 = r0, K = 32.0, θ=1000.0/log(10.0) ), 1)
    iterate_ratings = update_ratings(rule, nfl_ext_ratings, nfl_ext_competitions)
    # see \"Whos's #1\", Langville and Meyer, p.58
    required_output = Dict{String,Float64}(
   "New Orleans Saints" => 173.661,
   "Indianapolis Colts" => 170.331,
   "San Diego Chargers" => 127.582,
    "Minnesota Vikings" => 103.504,
       "Dallas Cowboys" => 89.1285,
  "Philadelphia Eagles" => 69.5331,
    "Green Bay Packers" => 67.8292,
    "Arizona Cardinals" => 53.227 ,
        "New York Jets" => 50.1431,
 "New England Patriots" => 39.6328,
       "Houston Texans" => 33.9024,
   "Cincinnati Bengals" => 33.0116,
     "Baltimore Ravens" => 32.0826,
      "Atlanta Falcons" => 28.1178,
  "Pittsburgh Steelers" => 27.1246,
     "Tennessee Titans" => 13.2216,
    "Carolina Panthers" => 11.4745,
  "San Francisco 49ers" => -1.28445,
      "New York Giants" => -5.3217,
       "Denver Broncos" => -11.1262,
       "Miami Dolphins" => -26.7173,
        "Chicago Bears" => -28.1416,
 "Jacksonville Jaguars" => -36.2142,
        "Buffalo Bills" => -53.3495,
     "Cleveland Browns" => -74.6639,
      "Oakland Raiders" => -83.3188,
     "Seattle Seahawks" => -88.8452,
   "Kansas City Chiefs" => -109.281,
  "Washington Redskins" => -110.212,
 "Tampa Bay Buccaneers" => -130.102,
        "Detroit Lions" => -170.809,
       "St. Louis Rams" => -194.119
                                           )
    for k in keys(required_output)
        @test abs(required_output[k] - iterate_ratings.ratings[k]) < 1.0e-2
    end

    # check recording doesn't break anything
    R5 = RatingsTable( nfl_ext_ratings.players, size(nfl_ext_competitions,1) )
    iterate_ratings2 = update_ratings(rule, nfl_ext_ratings, nfl_ext_competitions; record=R5)
    @test iterate_ratings2 == iterate_ratings
    @test R5[end:end] == RatingsTable(iterate_ratings2).ratings

    # extra I/O question
    test_out2 = "Tmp/test_ratings_table_out.csv"
    write_ratingstable( test_out2, R5 )
    R7 = read_ratingstable( test_out2 )
    @test R5 == R7

    # check it complains when you give invalid arguments
    @test_throws ArgumentError rule = UpdateIterate( UpdateColley(), 1)
    @test_throws ArgumentError rule = UpdateIterate( UpdateElo(), 0)
end

r0 = 0.0
n_samples = 1000
batch_size = 41
rule = UpdateSampleIterate( rule=UpdateElo(;r0 = r0, K = 32.0, θ=1000.0/log(10.0) ),
                            batch_size = batch_size,
                            n_samples = n_samples)
R6 = RatingsTable( nfl_ext_ratings.players, Int(ceil(n_samples/batch_size)) )
sample_ratings = update_ratings(rule, nfl_ext_ratings, nfl_ext_competitions; record=R6)
@testset "SampleIterate" begin
    @test R6[end:end] == RatingsTable(sample_ratings).ratings

    # check it complains when you give invalid arguments
    @test_throws ArgumentError rule = UpdateSampleIterate( UpdateColley(), 1, 0)
    @test_throws ArgumentError rule = UpdateSampleIterate( UpdateElo(), 0, 0)
    @test_throws ArgumentError rule = UpdateSampleIterate( UpdateElo(), 0, 1)
end

@testset "Scoring" begin
    # test normalised direction and score values
    for (i,s) in enumerate(scoring_rule_names)
        @testset "$s" begin
            nm = Meta.parse("$s(;normalise=true)")
            score_dir = :(score_direction( $nm ) )
            @test eval( score_dir ) == 1
            score_one = :(scoring_function( $nm, [0.0, 1.0], 2 ) )
            @test abs( eval( score_one  ) - 1.0) < 1.0e-6
            score_half = :(scoring_function( $nm, [0.5, 0.5], 2 ) )
            @test abs( eval( score_half ) - 0.0) < 1.0e-6
        end
    end
end

@testset "Predict" begin
    @testset "   outcomes" begin
        # cases that can predict outcome
        @test all( predict_outcome(UpdateElo(), 1.0, 1.0, missing, missing) .== (0.5, 0.5, 0.0) )
        @test all( predict_outcome(UpdateEloF(), 1.0, 1.0, missing, missing) .== (0.5, 0.5, 0.0) )
        @test all( predict_outcome(UpdateColley(), 1.0, 1.0, missing, missing) .== (0.0, 0.0, 1.0) )

        # case that use a sub-rule for prediction (default is Elo)
        @test all( predict_outcome(UpdateIterate(), 1.0, 1.0, missing, missing) .== (0.5, 0.5, 0.0) )
        @test all( predict_outcome(UpdateSampleIterate(), 1.0, 1.0, missing, missing) .== (0.5, 0.5, 0.0) )

        # cases based on score, not outcome
        @test_throws ErrorException predict_outcome(UpdateMassey(), 1.0, 1.0, missing, missing)
        @test_throws ErrorException predict_outcome(UpdateMasseyColley(), 1.0, 1.0, missing, missing)
        @test_throws ErrorException predict_outcome(UpdateKeenerScores(), 1.0, 1.0, missing, missing)

        # cases that don't have a predictive model
        @test_throws ErrorException predict_outcome(UpdateRevert(), 1.0, 1.0, missing, missing)
    end

    @testset "   margin" begin
        # not implemented yet
    end
end

ext: jl


text:
using SafeTestsets
@safetestset "aI Tests" begin include("aI_tests.jl") end


text:
using LinearAlgebra
using StaticArrays

abstract type AbstractUnitCell end
abstract type AbtractPeriodicCell <: AbstractUnitCell end
abstract type AbstractOrthorombicCell <: AbtractPeriodicCell end

struct NonPeriodic <: AbstractUnitCell end

struct VariableCell <: AbtractPeriodicCell
    cell::Vector{AbtractPeriodicCell}
end

struct CubicCell <: AbstractOrthorombicCell
    a::Float64
end

struct OrthorombicCell <: AbstractOrthorombicCell
    abc::SVector{3,Float64}
    OrthorombicCell(a::Real,b::Real,c::Real) = new([a,b,c])
    OrthorombicCell(a::Real) = new([a,a,a])
    function OrthorombicCell(a::AbstractVector)
        @assert length(a) == 3
        new(a)
    end
end

struct TriclinicCell <: AbtractPeriodicCell
    abc::SMatrix{3,3,Float64}
    function TriclinicCell(abc::AbstractMatrix)
        @assert size(abc) == (3,3)
        new(abc)
    end
    function TriclinicCell(a::AbstractVector)
        @assert length(a) == 3
        new(Diagonal(a))
    end
    function TriclinicCell(a::Real, b::Real, c::Real, α::Real, β::Real, γ::Real)
        tmp = zeros(3,3)
        tmp[1,1] = sind(γ)*a
        tmp[2,1] = cosd(γ)*a
        tmp[2,2] = b
        tmp[2,3] = cosd(α)*c
        tmp[1,3] = c*(cosd(β)-cosd(γ)*cosd(α))/sind(γ)
        tmp[3,3] = sqrt(c^2-tmp[1,3]^2-tmp[2,3]^2)
        new(tmp)
    end
end

parsecell(a::Number) = CubicCell(a)
function parsecell(a::Real, b::Real, c::Real)
    a == b == c && return CubicCell(a)
    return OrthorombicCell(a,b,c)
end

function parsecell(a::Real, b::Real, c::Real, α::Real, β::Real, γ::Real)
    if α == β == γ == 90
        a == b == c && return CubicCell(a)
        return OrthorombicCell(a,b,c)
    else
        return TriclinicCell(a,b,c,α,β,γ)
    end
end

Base.show(io::IO, c::CubicCell) = print(io, "CubicCell a=",c.a)
Base.show(io::IO, c::OrthorombicCell) = print(io, "OrthorombicCell")
Base.show(io::IO, c::TriclinicCell) = print(io, "TriclinicCell")
Base.show(io::IO, c::VariableCell) = print(io, "Variable cell of length ", length(c))

Base.length(c::VariableCell) = length(c.cell)
Base.getindex(c::VariableCell, i) = c.cell[i]
Base.lastindex(c::VariableCell) = length(c)
Base.firstindex(c::VariableCell) = 1

Base.iterate(c::VariableCell, state=1) = state > length(c) ? nothing : (c[state], state+1)

celldiag(c::AbtractPeriodicCell) = diag(cellmatrix(c))

cellmatrix(c::CubicCell) = Diagonal([c.a, c.a, c.a])
cellmatrix(c::OrthorombicCell) = Diagonal(c.abc)
cellmatrix(c::TriclinicCell) = c.abc
cellmatrix(c::VariableCell) = [ cellmatrix(x) for x in c]
cellvectorlengths(c::AbtractPeriodicCell) = [norm(x) for x in eachcol(cellmatrix(c))]
cellvectorlengths(c::AbstractOrthorombicCell) = (Diagonal ∘ cellmatrix)(c)
cellvectorlengths(c::VariableCell) = [ cellvectorlengths(x) for x in c ]

volume(c::CubicCell) = c.a^3
volume(c::OrthorombicCell) = prod(c.abc)
volume(c::AbtractPeriodicCell) = (abs ∘ det ∘ cellmatrix)(c)
volume(c::VariableCell) = [ volume(x) for x in c ]


text:
begin
local color_chooser_shader = TemplateProgram(
  joinpath(shaderdir, "colorchooser.vert"), joinpath(shaderdir, "colorchooser.frag"),
  fragdatalocation=[(0, "fragment_color"), (1, "fragment_groupid")]
)

local COLOR_QUAD = genquad(Vec3(0, 0, 0), Vec3(1, 0, 0), Vec3(0, 1, 0))

#GLPlot.toopengl{T <: AbstractRGB}(colorinput::Input{T}) = toopengl(lift(x->AlphaColorValue(x, one(T)), RGBA{T}, colorinput))

function visualize{X <: AbstractAlphaColorValue}(style::Style, color::X, data)
  screen       = data[:screen]
  camera       = screen.orthographiccam

  rdata = merge(@compat(Dict(
    :vertex => GLBuffer(COLOR_QUAD[1]),
    :uv     => GLBuffer(COLOR_QUAD[2]),
    :index  => indexbuffer(COLOR_QUAD[4]),
  )), data)

  rdata[:view]       = camera.view
  rdata[:projection] = camera.projection
  rdata[:color]      = color

  obj = RenderObject(rdata, color_chooser_shader, color_chooser_boundingbox)

  prerender!(obj, glEnable, GL_DEPTH_TEST, glDepthFunc, GL_LESS, glDisable, GL_CULL_FACE, enabletransparency)#
  postrender!(obj, render, obj.vertexarray) # Render the vertexarray

  # hover is true, if mouse
  hover = lift(SELECTION[:mouse_hover]) do selection
    selection[1][1] == obj.id
  end

  all_signals = foldl((tohsva(color), false, false, Vec2(0)), SELECTION[:mouse_hover]) do v0, selection

    hsv, hue_sat0, bright_trans0, mouse0 = v0
    mouse           = screen.inputs[:mouseposition].value
    mouse_clicked   = screen.inputs[:mousebuttonspressed].value

    hue_sat = in(0, mouse_clicked) && selection[1][1] == obj.id
    bright_trans = in(1, mouse_clicked) && selection[1][1] == obj.id

    if hue_sat && hue_sat0
      diff = mouse - mouse0
      hue = mod(hsv.c.h + diff[1], 360)
      sat = max(min(hsv.c.s + (diff[2] / 30.0), 1.0), 0.0)

      return (tohsva(hue, sat, hsv.c.v, hsv.alpha), hue_sat, bright_trans, mouse)
    elseif hue_sat && !hue_sat0
      return (hsv, hue_sat, bright_trans, mouse)
    end

    if bright_trans && bright_trans0
      diff    = mouse - mouse0
      brightness  = max(min(hsv.c.v - (diff[2]/100.0), 1.0), 0.0)
      alpha     = max(min(hsv.alpha + (diff[1]/100.0), 1.0), 0.0)

      return (tohsva(hsv.c.h, hsv.c.s, brightness, alpha), hue_sat0, bright_trans, mouse)
    elseif bright_trans && !bright_trans0
      return (hsv, hue_sat0, bright_trans, mouse)
    end

    return (hsv, hue_sat, bright_trans, mouse)
  end
  color1 = lift(x -> torgba(x[1]), all_signals)
  color1 = lift(x -> Vec4(x.c.r, x.c.g, x.c.b, x.alpha), Vec4, color1)
  hue_saturation = lift(x -> x[2], all_signals)
  brightness_transparency = lift(x -> x[3], all_signals)

  obj.uniforms[:color]                    = color1
  obj.uniforms[:hover]                    = hover
  obj.uniforms[:hue_saturation]           = hue_saturation
  obj.uniforms[:brightness_transparency]  = brightness_transparency

  return obj, color1
end

function color_chooser_boundingbox(obj)
  middle      = obj[:middle]
  swatchsize  = obj[:swatchsize]
  border_size = obj[:border_size]
  model       = obj[:model]
  verts       = COLOR_QUAD[1]
  minv = model*Vec4(minimum(verts)...,0f0)
  maxv = model*Vec4(maximum(verts)...,0f0)
  AABB(Vec3(minv[1:3]), Vec3(maxv[1:3]))
end
end # local begin color chooser


text:
module EvalMetrics

import Base: show, precision
import DocStringExtensions: SIGNATURES
import Statistics: quantile
import StatsBase: RealVector
using RecipesBase
using Reexport

include("encodings/Encodings.jl")
@reexport using .Encodings

include("utilities.jl")
include("confusion_matrix.jl")

const CMVector{T<:Real} = AbstractVector{ConfusionMatrix{T}}

include("metrics.jl")
include("thresholds.jl")
include("curves.jl")

export
    # confusion matrix
    ConfusionMatrix,

    # performance metrics
    true_positive,
    true_negative,
    false_positive,
    false_negative,
    true_positive_rate,
        sensitivity,
        recall,
        hit_rate,
    true_negative_rate,
        specificity,
        selectivity,
    false_positive_rate,
        fall_out,
        type_I_error,
    false_negative_rate,
        miss_rate,
        type_II_error,
    precision,
        positive_predictive_value,
    negative_predictive_value,
    false_discovery_rate,
    false_omission_rate,
    threat_score,
        critical_success_index,
    accuracy,
    balanced_accuracy,
    error_rate,
    balanced_error_rate,
    f1_score,
    fβ_score,
    matthews_correlation_coefficient,
        mcc,
    quant,
    topquant,
    positive_likelihood_ratio,
    negative_likelihood_ratio,
    diagnostic_odds_ratio,
    prevalence,

    # threshold functions
    thresholds,
    threshold_at_tpr,
    threshold_at_tnr,
    threshold_at_fpr,
    threshold_at_fnr,
    threshold_at_k,

    #curves
    auc_trapezoidal,
    curve, auc,
    PRCurve, prcurve, au_prcurve, prplot,
    ROCCurve, roccurve, au_roccurve, rocplot,

    # utilities
    binary_eval_report,
    mergesorted

binary_eval_report(targets::AbstractVector, scores::RealVector, fpr = 0.05) =
    binary_eval_report(current_encoding(), targets, scores, fpr)

function binary_eval_report(enc::TwoClassEncoding, targets::AbstractVector, scores::RealVector, fpr = 0.05)
    t = threshold_at_fpr(enc, targets, scores, fpr)
    c = ConfusionMatrix(enc, targets, scores, t)

    return Dict(
        "accuracy@fpr$(fpr)" => accuracy(c),
        "au_prcurve" => au_prcurve(enc, targets, scores),
        "au_roccurve" => au_roccurve(enc, targets, scores),
        "precision@fpr$(fpr)" => precision(c),
        "prevalence" => prevalence(c),
        "recall@fpr$(fpr)" => recall(c),
        "samples" => length(targets),
        "true negative rate@fpr$(fpr)" => true_negative_rate(c)
    )
end

end


text:
struct eFAST <: GSAMethod
    num_harmonics::Int
end
eFAST(; num_harmonics::Int=4) = eFAST(num_harmonics)

struct eFASTResult{T1}
    first_order::T1
    total_order::T1
end

function gsa(f,method::eFAST,p_range::AbstractVector;n::Int=1000,batch=false, distributed::Val{SHARED_ARRAY} = Val(false), kwargs...) where {SHARED_ARRAY}
    @unpack num_harmonics = method
    num_params = length(p_range)
    omega = [ (n-1) ÷ (2*num_harmonics) ]
    m = omega[1] ÷ (2*num_harmonics)

    if m >= num_params-1
        append!(omega, floor.(Int, collect(range(1,stop=m,length=num_params-1))))
    else
        append!(omega, collect(range(0,stop=num_params-2)) .% m .+1 )
    end

    omega_temp = similar(omega)
    s = (2/n) * (0:n-1)
    if SHARED_ARRAY
        ps = SharedMatrix{Float64}((num_params, n*num_params))
    else
        ps = Matrix{Float64}(undef, num_params, n*num_params)
    end
    @inbounds for i in 1:num_params
        omega_temp[i] = omega[1]
        for k ∈ 1:i-1
            omega_temp[k] = omega[k+1]
        end
        for k ∈ i+1:num_params
            omega_temp[k] = omega[k]
        end
        l = (i-1)*n+1:i*n
        phi = 2rand()
        for j in 1:num_params
            ps[j,l] .= quantile.(Uniform(p_range[j][1], p_range[j][2]), 0.5 .+ (1/pi) .* (asin.(sinpi.(omega_temp[j] .* s .+ phi))))
        end
    end

    if batch
        all_y = f(ps)
        multioutput = all_y isa AbstractMatrix
        y_size = nothing
        gsa_efast_all_y_analysis(method, all_y, num_params, y_size, n, omega, Val(multioutput))
    else
        _y = [f(ps[:,j]) for j in 1:size(ps,2)]
        multioutput = !(eltype(_y) <: Number)
        if eltype(_y) <: RecursiveArrayTools.AbstractVectorOfArray
            y_size = size(_y[1])
            __y = vec.(_y)
        else
            y_size = nothing
            __y = _y
        end
        if multioutput
            gsa_efast_all_y_analysis(method, reduce(hcat,__y), num_params, y_size, n, omega, Val(true))
        else
            gsa_efast_all_y_analysis(method, __y, num_params, y_size, n, omega, Val(false))
        end
    end
end
function gsa_efast_all_y_analysis(method, all_y, num_params, y_size, n, omega, ::Val{multioutput}) where {multioutput}
    @unpack num_harmonics = method
    if multioutput
        size_ = size(all_y)
        first_order = Vector{Vector{eltype(all_y)}}(undef, num_params)
        total_order = Vector{Vector{eltype(all_y)}}(undef, num_params)
    else
        first_order = Vector{eltype(all_y)}(undef, num_params)
        total_order = Vector{eltype(all_y)}(undef, num_params)
    end
    for i in 1:num_params
        if !multioutput
            ft = (fft(all_y[(i-1)*n+1:i*n]))[2:(n ÷ 2)]
            ys = abs2.(ft .* inv(n))
            varnce = 2*sum(ys)
            first_order[i] = 2*sum(ys[(1:num_harmonics)*Int(omega[1])])/varnce
            total_order[i] = 1 .- 2*sum(ys[1:(omega[1] ÷ 2)])/varnce
        else
            ys = Vector{Vector{eltype(all_y)}}(undef, size(all_y,1))
            varnce = Vector{eltype(all_y)}(undef, size(all_y,1))
            for j ∈ eachindex(ys)
                ff = fft(all_y[j,(i-1)*n+1:i*n])[2:(n ÷ 2)]
                ys[j] = ysⱼ = abs2.(ff .* inv(n))
                varnce[j] = 2*sum(ysⱼ)
            end
            first_order[i] = map((y,var) -> 2*sum(y[(1:num_harmonics)*(omega[1])])./var, ys, varnce)
            total_order[i] = map((y,var) -> 1 .- 2*sum(y[1:(omega[1] ÷ 2)])./var, ys, varnce)
        end
    end
    if isnothing(y_size)
        _first_order = reduce(hcat,first_order)
        _total_order = reduce(hcat,total_order)
    else
        f_shape = let y_size = y_size
            x -> [reshape(x[:,i], y_size) for i in 1:size(x,2)]
        end
        _first_order = reduce(hcat,map(f_shape,first_order))
        _total_order = reduce(hcat,map(f_shape,total_order))
    end
    return eFASTResult(_first_order, _total_order)
end


text:
# test uniqvalue, unionall, discretize, topnames
using TermWin
using DataArrays
using Base.Test

@test isna( uniqvalue( @data( [ 1,2,3 ] ) ) )
@test uniqvalue( @data( [ 2,2,2 ] ) ) == 2
@test uniqvalue( @data( [ 2,2,NA ] ) ) == 2
@test isna( uniqvalue( @data( [ ] ) ) )

@test isna( uniqvalue( @data( [ "","",NA ] ) ) )
@test uniqvalue( @data( [ "a","a", "", NA ] ) ) == "a"
@test isna( uniqvalue( @data( [ "a","b", "", NA ] ) ) )

arr = unionall( @data( Array[ [1], [2,3] ] ) )
sort!( arr )
@test arr == [1,2,3]
arr = DataArray( Array{Int,1}[ Int[1], Int[2,3] ] )
push!(arr, NA)
arr = unionall(arr)
sort!(arr)
@test arr == [1,2,3]

# cut_categories
@test TermWin.cut_categories( Float64, [ 0,1,2,3] ) ==
[
"1. <0",
"2. [0,1)",
"3. [1,2)",
"4. [2,3)",
"5. ≥3",
]
@test TermWin.cut_categories( Int, [ 0,1,2,3] ) ==
[
"1. ≤-1",
"2. 0",
"3. 1",
"4. 2",
"5. 3+"
]
@test TermWin.cut_categories( Float64, [ 0,1,2,3], boundedness = :bounded ) ==
[
"1. [0,1)",
"2. [1,2)",
"3. [2,3]"
]
@test TermWin.cut_categories( Float64, [ 0,1,2,3], boundedness = :bounded ) ==
[
"1. [0,1)",
"2. [1,2)",
"3. [2,3]"
]
@test TermWin.cut_categories( Float64, [ 0,1,2,3], boundedness = :boundedbelow ) ==
[
"1. [0,1)",
"2. [1,2)",
"3. [2,3)",
"4. ≥3",
]
@test TermWin.cut_categories( Float64, [ 0,1,2,3], boundedness = :boundedabove ) ==
[
"1. <0",
"2. [0,1)",
"3. [1,2)",
"4. [2,3]"
]
@test TermWin.cut_categories( Float64, [ 0,1,2,3], boundedness = :bounded, leftequal = false) ==
[
"1. [0,1]",
"2. (1,2]",
"3. (2,3]"
]
@test TermWin.cut_categories( Float64, [ 0.0,.25,0.5,0.75,1.0], boundedness = :bounded, scale=100.0, suffix="%" ) ==
[
"1. [0%,25%)",
"2. [25%,50%)",
"3. [50%,75%)",
"4. [75%,100%]"
]

arr = [-1.0,0.0,0.5,1.,1.5,2.,3.]
result = discretize( arr, [0,1,2] )
@test result[1] == "1. <0"
@test result[2] == "2. [0,1)"
@test result[3] == "2. [0,1)"
@test result[4] == "3. [1,2)"
@test result[5] == "3. [1,2)"
@test result[6] == "4. ≥2"
@test result[7] == "4. ≥2"

result = discretize( arr, [0,1,2], boundedness = :bounded )
@test isna( result[1] )
@test result[2] == "1. [0,1)"
@test result[3] == "1. [0,1)"
@test result[4] == "2. [1,2]"
@test result[5] == "2. [1,2]"
@test result[6] == "2. [1,2]"
@test isna( result[7] )

result = discretize( arr, [0,1,2], boundedness = :boundedbelow )
@test isna( result[1] )
@test result[2] == "1. [0,1)"
@test result[3] == "1. [0,1)"
@test result[4] == "2. [1,2)"
@test result[5] == "2. [1,2)"
@test result[6] == "3. ≥2"
@test result[7] == "3. ≥2"

result = discretize( arr, [0,1,2], boundedness = :boundedabove )
@test result[1] == "1. <0"
@test result[2] == "2. [0,1)"
@test result[3] == "2. [0,1)"
@test result[4] == "3. [1,2]"
@test result[5] == "3. [1,2]"
@test result[6] == "3. [1,2]"
@test isna( result[7] )

arr = @data( [ 1,2,3, NA] )
result = discretize( arr, [ 0,1,2,3])
@test result[1] == "3. 1"
@test result[2] == "4. 2"
@test result[3] == "5. 3+"
@test isna( result[4] )

arr = @data( [ -1.0, 1.0,2.0,3.0, NA] )
result = discretize( arr, [ 0,1,2,3])
@test result[1] == "1. <0"
@test result[2] == "3. [1,2)"
@test result[3] == "4. [2,3)"
@test result[4] == "5. ≥3"
@test isna( result[5] )

arr = @data( [ -1.0, 1.0,2.0,3.0, NA] )
result = discretize( arr, [ 0,1,2,3], label="x")
@test result[1] == "1. x < 0"
@test result[2] == "3. 1 ≤ x < 2"
@test result[3] == "4. 2 ≤ x < 3"
@test result[4] == "5. 3 ≤ x"
@test isna( result[5] )

arr = @data( [ -1.0, 1.0,2.0,3.0, NA] )
result = discretize( arr, [ 0,1,2,3], label="x", absolute=true)
@test result[1] == "3. 1 ≤ |x| < 2"
@test result[2] == "3. 1 ≤ |x| < 2"
@test result[3] == "4. 2 ≤ |x| < 3"
@test result[4] == "5. 3 ≤ |x|"
@test isna( result[5] )

name=@data( [ "Alice", "Bob", "Jane", "Joe" ] )
score=@data( [ 7, 8, 9, 5 ] )
result = topnames( name, score, 2 )
@test result[ 1 ] == "Others"
@test result[ 2 ] == "2. Bob"
@test result[ 3 ] == "1. Jane"
@test result[ 4 ] == "Others"

name=@data( [ "Alice", "Bob", "Jane", "Joe" ] )
score=@data( [ 7, 9, 9, 5 ] )
result = topnames( name, score, 2 )
@test result[ 1 ] == "2. Alice"
@test result[ 2 ] == "1. Bob"
@test result[ 3 ] == "1. Jane"
@test result[ 4 ] == "Others"

# Note that in actual usage "name" is always sorted
# see CalcPivot constructor's handling of topnames' "by" argument
name=@data( [ "Alice", "Bob", "Jane", "Joe" ] )
score=@data( [ 7, 9, 9, 5 ] )
result = topnames( name, score, 2, dense=false )
@test result[ 1 ] == "Others" # Alice would have been 3rd
@test result[ 2 ] == "1. Bob"
@test result[ 3 ] == "1. Jane"
@test result[ 4 ] == "Others"

name=@data( [ "Alice", "Bob", "Jane", "Joe" ] )
score=@data( [ 7, 8, -9, 5 ] )
result = topnames( name, score, 2, absolute=true )
@test result[ 1 ] == "Others"
@test result[ 2 ] == "1. Bob" # because it's positive
@test result[ 3 ] == "2. Jane" # because it's by magnitude the largest
@test result[ 4 ] == "Others"

name=@data( [ "Alice", "Bob", "Jane", "Joe" ] )
score=@data( [ 7, 8, -9, 5 ] )
result = topnames( name, score, 2, absolute=true, parens=true )
@test result[ 1 ] == "Others"
@test result[ 2 ] == "1. Bob" # because it's positive
@test result[ 3 ] == "2. (Jane)" # because it's by magnitude the largest
@test result[ 4 ] == "Others"


text:
"""
测试Gamma4相关的功能
"""

using Test
using MARY_fRG.Fermi
using MARY_fRG.FlowEquation
using MARY_fRG.Drawers
using MARY_fRG.Fermi.Surface

@testset "普通的正方格子" begin
    qs = common_square_lattice(0.)
    Γ4 = ECGamma4(
        qs, 8.0, 16, 50
    )
    plt = draw_lines(Γ4.model.brillouin.edges)
    sf, sfpidx = const_energy_line_in_patches(Γ4.ltris, Γ4.ladjs,
    Γ4.lpats, 0., Γ4.model.dispersion[1])
    draw_lines!(plt, sf; lcidx=sfpidx)
    patcs = Γ4.patches[1]
    draw_points!(plt, patcs; text=1:1:length(patcs))
    draw_polygon!(plt, Γ4.ltris; pcidx=Γ4.lpats)
    @savepng plt "Gamma4"
    @test Γ4.k4tab[1, 1, 1, 1, 1, 12, 12] == 1
end

@testset "普通的三角格子" begin
    ts = common_triangle_lattice(0.)
    Γ4 = ECGamma4(
        ts, 6.0, 18, 50
    )
    plt = draw_lines(Γ4.model.brillouin.edges)
    sf, sfpidx = const_energy_line_in_patches(Γ4.ltris, Γ4.ladjs,
    Γ4.lpats, 0., Γ4.model.dispersion[1])
    draw_lines!(plt, sf; lcidx=sfpidx)
    patcs = Γ4.patches[1]
    draw_points!(plt, patcs; text=1:1:length(patcs))
    draw_polygon!(plt, Γ4.ltris; pcidx=Γ4.lpats)
    @savepng plt "Gamma42"
    @test Γ4.k4tab[1, 1, 1, 1, 1, 12, 12] == 1
end

@testset "两个能带的正方格子" begin
    disp1(x, y) = -2*(cos(x) + cos(y))
    disp2(x, y) = -2*(cos(x) + cos(y)) + 1.0
    qs = QuadrateSystem{:DOUBLE}([disp1, disp2])
    Γ4 = ECGamma4(
        qs, 8.0, 16, 50
    )
    plt = draw_lines(Γ4.model.brillouin.edges)
    #
    sf, sfpidx = const_energy_line_in_patches(Γ4.ltris, Γ4.ladjs,
    Γ4.lpats, 0., Γ4.model.dispersion[1])
    draw_lines!(plt, sf; lcidx=sfpidx)
    #
    sf2, sf2pidx = const_energy_line_in_patches(Γ4.ltris, Γ4.ladjs,
    Γ4.lpats, 0., Γ4.model.dispersion[2])
    draw_lines!(plt, sf2; lcidx=sf2pidx)
    #
    patcs = Γ4.patches[1]
    draw_points!(plt, patcs; text=1:1:length(patcs))
    patcs2 = Γ4.patches[2]
    draw_points!(plt, patcs2; text=1:1:length(patcs))
    draw_polygon!(plt, Γ4.ltris; pcidx=Γ4.lpats)
    @savepng plt "Gamma43"
    @test Γ4.k4tab[1, 1, 2, 2, 1, 9, 1] == 9
end


text:
# Text provided under a Creative Commons Attribution license, CC-BY.  All code is made available under the FSF-approved BSD-3 license.  (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784.
# [@LorenaABarba](https://twitter.com/LorenaABarba)
# Translation to Julia by [@miguelraz_](https://twitter.com/miguelraz_)

# 12 steps to Navier–Stokes
# =====
# ***

# You should have completed your own code for [Step 5](./07_Step_5.ipynb) before continuing to this lesson. As with Steps 1 to 4, we will build incrementally, so it's important to complete the previous step!

# We continue ...

# Step 6: 2-D Convection
# ----
# ***

# Now we solve 2D Convection, represented by the pair of coupled partial differential equations below:

# $$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = 0$$

# $$\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = 0$$

# Discretizing these equations using the methods we've applied previously yields:

# $$\frac{u_{i,j}^{n+1}-u_{i,j}^n}{\Delta t} + u_{i,j}^n \frac{u_{i,j}^n-u_{i-1,j}^n}{\Delta x} + v_{i,j}^n \frac{u_{i,j}^n-u_{i,j-1}^n}{\Delta y} = 0$$

# $$\frac{v_{i,j}^{n+1}-v_{i,j}^n}{\Delta t} + u_{i,j}^n \frac{v_{i,j}^n-v_{i-1,j}^n}{\Delta x} + v_{i,j}^n \frac{v_{i,j}^n-v_{i,j-1}^n}{\Delta y} = 0$$

# Rearranging both equations, we solve for $u_{i,j}^{n+1}$ and $v_{i,j}^{n+1}$, respectively.  Note that these equations are also coupled.

# $$u_{i,j}^{n+1} = u_{i,j}^n - u_{i,j} \frac{\Delta t}{\Delta x} (u_{i,j}^n-u_{i-1,j}^n) - v_{i,j}^n \frac{\Delta t}{\Delta y} (u_{i,j}^n-u_{i,j-1}^n)$$

# $$v_{i,j}^{n+1} = v_{i,j}^n - u_{i,j} \frac{\Delta t}{\Delta x} (v_{i,j}^n-v_{i-1,j}^n) - v_{i,j}^n \frac{\Delta t}{\Delta y} (v_{i,j}^n-v_{i,j-1}^n)$$

# ### Initial Conditions

# The initial conditions are the same that we used for 1D convection, applied in both the x and y directions.

# $$u,\ v\ = \begin{cases}\begin{matrix}
# 2 & \text{for } x,y \in (0.5, 1)\times(0.5,1) \cr
# 1 & \text{everywhere else}
# \end{matrix}\end{cases}$$

# ### Boundary Conditions

# The boundary conditions hold u and v equal to 1 along the boundaries of the grid
# .

# $$u = 1,\ v = 1 \text{ for } \begin{cases} \begin{matrix}x=0,2\cr y=0,2 \end{matrix}\end{cases}$$

using Plots

### Variable declarations

nₓ = 101
ny = 101
nₜ = 80
c = 1
Δx = 2 / (nₓ - 1)
Δy = 2 / (ny - 1)
σ = 0.2
Δt = σ * Δx

rangex = range(start = 0, stop = 2, length = nₓ)
rangey = range(start = 0, stop = 2, length = ny)
indices = CartesianIndices((rangex, rangey))

uₙ = ones(ny, nₓ) ##create a 1xn vector of 1's
vₙ = ones(ny, nₓ)
u  = ones(ny, nₓ)
v  = ones(ny, nₓ)

###Assign initial conditions
##set hat function I.C. : u(.5<=x<=1 && .5<=y<=1 ) is 2
u[indices] .= 2.0
##set hat function I.C. : v(.5<=x<=1 && .5<=y<=1 ) is 2
v[indices] .= 2.0

# TODO PLOT
# Why is [meshgrid](https://groups.google.com/g/julia-users/c/83Pfg9HGhGQ/m/9G_0wi-GBQAJ?pli=1) not an efficient idea?
# Summary: Just use a 2D array comprehension
fig = pyplot.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')
X, Y = numpy.meshgrid(x, y)

ax.plot_surface(X, Y, u, cmap=cm.viridis, rstride=2, cstride=2)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$');

# TODO BIKESHED Syntax
# Possible OffsetArrays.jl example?
for n in  1:nₜ
    un = copy(u)
    vn = copy(v)
    @views u[begin+1:end, begin+1:end] = (un[begin+1:end, begin+1:end] -
                 (un[begin+1:end, begin+1:end] * c * Δt / Δx * (un[begin+1:end, begin+1:end] - un[begin+1:end, begin:end-1])) -
                  vn[begin+1:end, begin+1:end] * c * Δt / Δy * (un[begin+1:end, begin+1:end] - un[:-1, begin+1:end]))
    @views v[begin+1:end, begin+1:end] = (vn[begin+1):end, begin+1:end] -
                 (un[begin+1:end, begin+1:end] * c * Δt / Δx * (vn[begin+1:end, begin+1:end] - vn[begin+1:end, begin:end-1])) -
                  vn[begin+1:end, begin+1:end] * c * Δt / Δy * (vn[begin+1:end, begin+1:end] - vn[begin:end-1, begin+1:end]))

    u[begin, :] = 1
    u[end, :] = 1
    u[:, begin] = 1
    u[:, end] = 1

    v[begin, :] = 1
    v[end, :] = 1
    v[:, begin] = 1
    v[:, end] = 1
end

# TODO PLOT
fig = pyplot.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')
X, Y = numpy.meshgrid(x, y)

ax.plot_surface(X, Y, u, cmap=cm.viridis, rstride=2, cstride=2)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$');

# TODO PLOT
fig = pyplot.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')
X, Y = numpy.meshgrid(x, y)
ax.plot_surface(X, Y, v, cmap=cm.viridis, rstride=2, cstride=2)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$');

## Learn More

# The video lesson that walks you through the details for Steps 5 to 8 is **Video Lesson 6** on You Tube:

    # TODO YOUTUBE
    # src="https://www.youtube.com/embed/tUg_dE3NXoY"

max_stars_repo_name: mkitti/CFDJulia

max_issues_repo_name: mkitti/CFDJulia

max_forks_repo_name: mkitti/CFDJulia


text:
function eracompile(
    init::Dict, eroot::Dict;
    modID::AbstractString, parID::AbstractString,
    regID::AbstractString="GLB", timeID::Union{Integer,Vector}=0,
    gres::Real=0, plvls::Union{AbstractString,Integer,Vector{<:Real}}
)

    emod,epar,ereg,etime = erainitialize(
        init;
        modID=modID,parID=parID,regID=regID,timeID=timeID,
        gres=gres
    );

    if typeof(plvls) <: Array
          for p in plvls; epar["level"] = p; eracompile(emod,epar,ereg,etime,eroot); end
    else; epar["level"] = plvls; eracompile(emod,epar,ereg,etime,eroot);
    end

end

function eracompile(
    emod::AbstractDict, epar::AbstractDict, ereg::AbstractDict, etime::AbstractDict,
    eroot::AbstractDict
)

    @info "$(Dates.now()) - Preallocating arrays ..."
    nlon,nlat = ereg["size"]; nt = etime["End"] + 1 - etime["Begin"]
    eavg = zeros(nlon,nlat);
    emax = zeros(nlon,nlat); emin = zeros(nlon,nlat); erng = zeros(nlon,nlat);
    edhr = zeros(nlon,nlat); eitr = zeros(nlon,nlat); esea = zeros(nlon,nlat);

    @info "$(Dates.now()) - Extracting preliminarily-analyzed reanalysis data ..."
    for yr = etime["Begin"] : etime["End"]

        eds,evar = eraanaread(
            "domain_yearly_mean_climatology",
            emod,epar,ereg,eroot,Date(yr)
        )
        eavg += evar[:]
        if yr == etime["Begin"]
            emax += evar[:]; emin += evar[:]
        else
            emax .= max.(evar[:],emax)
            emin .= min.(evar[:],emin)
        end
        close(eds)

        ds1,rmax = eraanaread(
            "domain_yearly_maximum_climatology",
            emod,epar,ereg,eroot,Date(yr)
        )
        ds2,rmin = eraanaread(
            "domain_yearly_minimum_climatology",
            emod,epar,ereg,eroot,Date(yr)
        )
        erng += rmax[:] .- rmin[:]
        close(ds1); close(ds2)

        eds,evar = eraanaread(
            "domain_monthly_mean_climatology",
            emod,epar,ereg,eroot,Date(yr)
        )
        esea += dropdims(maximum(evar[:],dims=3) .- minimum(evar[:],dims=3),dims=3)
        close(eds)

        ds1,rmax = eraanaread(
            "domain_monthly_maximum_climatology",
            emod,epar,ereg,eroot,Date(yr)
        )
        ds2,rmin = eraanaread(
            "domain_monthly_minimum_climatology",
            emod,epar,ereg,eroot,Date(yr)
        )
        eitr += dropdims(mean(rmax[:] .- rmin[:],dims=3),dims=3)
        close(ds1); close(ds2)

        eds,evar = eraanaread(
            "domain_yearly_mean_hourly",
            emod,epar,ereg,eroot,Date(yr)
        )
        edhr += dropdims(maximum(evar[:],dims=3) .- minimum(evar[:],dims=3),dims=3)
        close(eds)

    end

    @info "$(Dates.now()) - Calculating yearly mean, and diurnal, seasonal and interannual variability ..."
    eavg = eavg/nt; eian = emax .- emin
    esea = esea/nt; eitr = eitr/nt; edhr = edhr/nt
    erng = erng/nt .- (esea .+ eitr)
    eracmpsave(eavg,edhr,eitr,esea,eian,emod,epar,ereg,eroot)

end

function eracmpsave(
    eavg::Array{<:Real,2}, edhr::Array{<:Real,2},
    eitr::Array{<:Real,2}, esea::Array{<:Real,2}, eian::Array{<:Real,2},
    emod::Dict, epar::Dict, ereg::Dict, eroot::Dict
)

    @info "$(Dates.now()) - Saving compiled $(uppercase(emod["dataset"])) $(epar["name"]) data in $(gregionfullname(ereg["region"])) (Horizontal Resolution: $(ereg["step"])) ..."

    cfol = eraanafolder(epar,ereg,eroot); fcmp = eracmpname(emod,epar,ereg);
    cfnc = joinpath(cfol,fcmp);

    if isfile(cfnc)
        @info "$(Dates.now()) - Stale NetCDF file $(cfnc) detected.  Overwriting ..."
        rm(cfnc);
    end

    @debug "$(Dates.now()) - Creating NetCDF file $(afnc) for compiled $(emod["dataset"]) $(epar["name"]) data ..."

    ds = Dataset(cfnc,"c");
    ds.dim["longitude"] = ereg["size"][1]
    ds.dim["latitude"]  = ereg["size"][2]

    nclon = defVar(ds,"longitude",Float64,("longitude",),attrib = Dict(
        "units"     => "degrees_east",
        "long_name" => "longitude",
    ))

    nclat = defVar(ds,"latitude",Float64,("latitude",),attrib = Dict(
        "units"     => "degrees_north",
        "long_name" => "latitude",
    ))

    nclon[:] = ereg["lon"];
    nclat[:] = ereg["lat"];

    ncavg = defVar(ds,"average",Float32,("longitude","latitude"),
        attrib = Dict(
            "long_name" => epar["era5"],
            "full_name" => epar["name"],
            "units"     => epar["unit"],
            "level"     => epar["level"],
    ))

    ncian = defVar(ds,"variability_interannual",Float32,("longitude","latitude"),
        attrib = Dict(
            "long_name" => epar["era5"],
            "full_name" => epar["name"],
            "units"     => epar["unit"],
            "level"     => epar["level"],
    ))

    ncsea = defVar(ds,"variability_seasonal",Float32,("longitude","latitude"),
        attrib = Dict(
            "long_name" => epar["era5"],
            "full_name" => epar["name"],
            "units"     => epar["unit"],
            "level"     => epar["level"],
    ))

    ncitr = defVar(ds,"variability_intraseasonal",Float32,("longitude","latitude"),
        attrib = Dict(
            "long_name" => epar["era5"],
            "full_name" => epar["name"],
            "units"     => epar["unit"],
            "level"     => epar["level"],
    ))

    ncdhr = defVar(ds,"variability_diurnal",Float32,("longitude","latitude"),
        attrib = Dict(
            "long_name" => epar["era5"],
            "full_name" => epar["name"],
            "units"     => epar["unit"],
            "level"     => epar["level"],
    ))

    ncavg[:] = eavg; ncian[:] = eian; ncsea[:] = esea; ncitr[:] = eitr; ncdhr[:] = edhr

    @info "$(Dates.now()) - Compiled $(uppercase(emod["dataset"])) $(epar["name"]) in $(gregionfullname(ereg["region"])) (Horizontal Resolution: $(ereg["step"])) has been saved into file $(cfnc) and moved to the data directory $(cfol)."

end


text:
export Evaluation,ivp,bvp,Dirichlet,Neumann

## Evaluation constructors

abstract type Evaluation{T}<:Operator{T} end

@functional Evaluation

# M = leftendpoint/rightendpoint if endpoint
struct ConcreteEvaluation{S,M,OT,T} <: Evaluation{T}
    space::S
    x::M
    order::OT
end

ConcreteEvaluation(sp::Space,x,o::Number) =
    ConcreteEvaluation{typeof(sp),typeof(x),typeof(o),rangetype(sp)}(sp,x,o)

Evaluation(::Type{T},sp::Space,x,order) where {T} =
    ConcreteEvaluation{typeof(sp),typeof(x),typeof(order),T}(sp,x,order)
# TODO: This seems like a bad idea: if you are specifying x, just go with the generic version
function Evaluation(::Type{T},sp::UnivariateSpace,x::Number,order) where {T}
    d=domain(sp)
    if isa(d,IntervalOrSegment) && isapprox(leftendpoint(d),x)
        Evaluation(T,sp,leftendpoint,order)
    elseif isa(d,IntervalOrSegment) && isapprox(rightendpoint(d),x)
        Evaluation(T,sp,rightendpoint,order)
    else
        ConcreteEvaluation{typeof(sp),typeof(x),typeof(order),T}(sp,x,order)
    end
end

Evaluation(sp::Space,x,order) = Evaluation(rangetype(sp),sp,x,order)

Evaluation(d::Space,x::Union{Number,typeof(leftendpoint),typeof(rightendpoint)}) = Evaluation(d,x,0)
Evaluation(::Type{T},d::Space,n...) where {T} = error("Override Evaluation for $(typeof(d))")
Evaluation(::Type{T},d,n...) where {T} = Evaluation(T,Space(d),n...)
Evaluation(S::Space,n...) = error("Override Evaluation for $(typeof(S))")
Evaluation(d,n...) = Evaluation(Space(d),n...)
Evaluation(x::Union{Number,typeof(leftendpoint),typeof(rightendpoint)}) = Evaluation(UnsetSpace(),x,0)
Evaluation(x::Union{Number,typeof(leftendpoint),typeof(rightendpoint)},k::Integer) =
    Evaluation(UnsetSpace(),x,k)

rangespace(E::ConcreteEvaluation{<:AmbiguousSpace}) = ConstantSpace()
rangespace(E::ConcreteEvaluation) = ConstantSpace(Point(E.x))

function convert(::Type{Operator{T}},E::ConcreteEvaluation) where T
    if T == eltype(E)
        E
    else
        ConcreteEvaluation{typeof(E.space),typeof(E.x),typeof(E.order),T}(E.space,E.x,E.order)
    end
end

## default getindex
getindex(D::ConcreteEvaluation,k::Integer) =
    eltype(D)(differentiate(Fun(D.space,[zeros(eltype(D),k-1);one(eltype(D))]),D.order)(D.x))

#special leftendpoint/rightendpoint overrides
for (dop, fop) in ((:leftendpoint,:first), (:rightendpoint,:last))
    @eval begin
        rangespace(E::ConcreteEvaluation{<:AmbiguousSpace,typeof($dop)}) = UnsetSpace()
        function rangespace(E::ConcreteEvaluation{S,typeof($dop)}) where {S}
            d = domain(domainspace(E))
            isambiguous(d) && return ConstantSpace()
            return ConstantSpace(Point($dop(d)))
        end
        function getindex(D::ConcreteEvaluation{S,typeof($dop)},k::Integer) where {S}
            T=eltype(D)
            T($fop(differentiate(Fun(D.space,[zeros(T,k-1);one(T)]),D.order)))
        end
    end
end

## EvaluationWrapper

struct EvaluationWrapper{S<:Space,M,FS<:Operator,OT,T<:Number} <: Evaluation{T}
    space::S
    x::M
    order::OT
    op::FS
end

@wrapper EvaluationWrapper
EvaluationWrapper(sp::Space,x,order,func::Operator) =
    EvaluationWrapper{typeof(sp),typeof(x),typeof(func),typeof(order),eltype(func)}(sp,x,order,func)

domainspace(E::Evaluation) = E.space
domain(E::Evaluation) = domain(E.space)
promotedomainspace(E::Evaluation,sp::Space) = Evaluation(sp,E.x,E.order)

function convert(::Type{Operator{T}},E::EvaluationWrapper) where T
    if T == eltype(E)
        E
    else
        EvaluationWrapper(E.space,E.x,E.order,convert(Operator{T},E.op))::Operator{T}
    end
end

## Convenience routines

evaluate(d::Domain,x) = Evaluation(d,x)
ldiffbc(d,k) = Evaluation(d,leftendpoint,k)
rdiffbc(d,k) = Evaluation(d,rightendpoint,k)

ldirichlet(d) = ldiffbc(d,0)
rdirichlet(d) = rdiffbc(d,0)
lneumann(d) = ldiffbc(d,1)
rneumann(d) = rdiffbc(d,1)

ivp(d,k) = Operator{prectype(d)}[ldiffbc(d,i) for i=0:k-1]
bvp(d,k) = vcat(Operator{prectype(d)}[ldiffbc(d,i) for i=0:div(k,2)-1],
                Operator{prectype(d)}[rdiffbc(d,i) for i=0:div(k,2)-1])

periodic(d,k) = Operator{prectype(d)}[Evaluation(d,leftendpoint,i)-Evaluation(d,rightendpoint,i) for i=0:k]

# shorthand for second order
ivp(d) = ivp(d,2)
bvp(d) = bvp(d,2)

for op in (:rdirichlet,:ldirichlet,:lneumann,:rneumann,:ivp,:bvp)
    @eval begin
        $op() = $op(UnsetSpace())
    end
end

for op in (:ldiffbc,:rdiffbc,:ivp,:bvp,:periodic)
    @eval $op(k::Integer) = $op(UnsetSpace(),k)
end

abstract type Dirichlet{S,T} <: Operator{T} end

struct ConcreteDirichlet{S,V,T} <: Dirichlet{S,T}
    domainspace::S
    rangespace::V
    order::Int
end

ConcreteDirichlet(sp::Space,rs::Space,order) =
    ConcreteDirichlet{typeof(sp),typeof(rs),rangetype(sp)}(sp,rs,order)
ConcreteDirichlet(sp::Space,order) = ConcreteDirichlet(sp,Space(∂(domain(sp))),order)
ConcreteDirichlet(sp::Space) = ConcreteDirichlet(sp,0)

convert(::Type{Operator{T}},B::ConcreteDirichlet{S,V}) where {S,V,T} =
    ConcreteDirichlet{S,V,T}(B.domainspace,B.rangespace,B.order)

struct DirichletWrapper{S,T} <: Dirichlet{S,T}
    op::S
    order::Int
end

@wrapper DirichletWrapper

DirichletWrapper(B::Operator,λ=0) = DirichletWrapper{typeof(B),eltype(B)}(B,λ)

convert(::Type{Operator{T}},B::DirichletWrapper) where {T} =
    DirichletWrapper(Operator{T}(B.op),B.order)::Operator{T}

# Default is to use diffbca
default_Dirichlet(sp::Space,λ) = DirichletWrapper([ldiffbc(sp,λ);rdiffbc(sp,λ)],λ)
Dirichlet(sp::Space,λ) = default_Dirichlet(sp,λ)
Dirichlet(sp::Space) = Dirichlet(sp,0)
Dirichlet() = Dirichlet(UnsetSpace())

Dirichlet(d::Domain,λ...) = Dirichlet(Space(d),λ...)
Neumann(sp::Space) = Dirichlet(sp,1)
Neumann(d::Domain) = Neumann(Space(d))
Neumann() = Dirichlet(UnsetSpace(),1)

domainspace(B::ConcreteDirichlet) = B.domainspace
rangespace(B::ConcreteDirichlet) = B.rangespace

promotedomainspace(E::Dirichlet,sp::Space) = Dirichlet(sp,E.order)

"""
`Evaluation(sp,x,k)` is the functional associated with evaluating the
`k`-th derivative at a point `x` for the space `sp`.
"""
Evaluation(sp::Space,x,k)

"""
`Evaluation(sp,x)` is the functional associated with evaluating
at a point `x` for the space `sp`.
"""
Evaluation(sp::Space,x)

"""
`Evaluation(x)` is the functional associated with evaluating
at a point `x`.
"""
Evaluation(x)

"""
`Dirichlet(sp,k)` is the operator associated with restricting the
`k`-th derivative on the boundary for the space `sp`.
"""
Dirichlet(sp::Space,k)

"""
`Dirichlet(sp)` is the operator associated with restricting the
 the boundary for the space `sp`.
"""
Dirichlet(sp::Space)

"""
`Dirichlet()` is the operator associated with restricting on the
 the boundary.
"""
Dirichlet()

"""
`Neumann(sp)` is the operator associated with restricting the
normal derivative on the boundary for the space `sp`.
At the moment it is implemented as `Dirichlet(sp,1)`.
"""
Neumann(sp::Space)

"""
`Neumann( is the operator associated with restricting the
normal derivative on the boundary.
"""
Neumann()


text:
####################  Wrapper Mask  ####################

struct FlipMask{M} <: AbstractWrapperMask
    mask::M
end

Base.:!(m::AbstractAttenMask) = FlipMask(m)
Base.:!(m::FlipMask) = m.mask

@init @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
    Adapt.adapt(to::CUDA.Adaptor, m::FlipMask) = Indexer{typeof(m)}((mask = adapt(to, m.mask),))
end

adapt_structure(to, x::FlipMask) = FlipMask(adapt(to, x.mask))
GetIndexer(m::FlipMask) = Indexer{typeof(m)}((mask = GetIndexer(m.mask),))

Base.@propagate_inbounds Base.getindex(m::Indexer{<:FlipMask}, I::Integer...) = !m.mask[I...]

check_constraint(m::FlipMask, x) = check_constraint(m.mask, x)

AxesConstraint(m::FlipMask) = AxesConstraint(m.mask)

Base.show(io::IO, m::FlipMask) = (print(io, '!'); show(io, m.mask); io)

Broadcast.broadcasted(::MaskStyle, !, m::AbstractAttenMask) = FlipMask(m)

struct CombinedMask{C, Ts<:Tuple} <: AbstractWrapperMask
    f::C
    masks::Ts
end

Base.:|(m1::AbstractAttenMask, m2::AbstractAttenMask) = CombinedMask(|, (m1, m2))
Base.:|(m1::CombinedMask{typeof(|)}, m2::AbstractAttenMask) = CombinedMask(|, (m1.masks..., m2))
Base.:|(m1::AbstractAttenMask, m2::CombinedMask{typeof(|)}) = CombinedMask(|, (m1, m2.masks...))
Base.:|(m1::CombinedMask{typeof(|)}, m2::CombinedMask{typeof(|)}) = CombinedMask(|, (m1.masks..., m2.masks...))
Base.:&(m1::AbstractAttenMask, m2::AbstractAttenMask) = CombinedMask(&, (m1, m2))
Base.:&(m1::CombinedMask{typeof(&)}, m2::AbstractAttenMask) = CombinedMask(&, (m1.masks..., m2))
Base.:&(m1::AbstractAttenMask, m2::CombinedMask{typeof(&)}) = CombinedMask(&, (m1, m2.masks...))
Base.:&(m1::CombinedMask{typeof(&)}, m2::CombinedMask{typeof(&)}) = CombinedMask(&, (m1.masks..., m2.masks...))

@init @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
    Adapt.adapt(to::CUDA.Adaptor, m::CombinedMask) = Indexer{typeof(m)}((f = adapt(to, m.f),
                                                                         masks = map(Base.Fix1(adapt, to), m.masks)))
end

adapt_structure(to, x::CombinedMask) = CombinedMask(x.f, adapt(to, x.masks))
GetIndexer(m::CombinedMask) = Indexer{typeof(m)}((m.f, masks = map(GetIndexer, m.masks)))

@inline function _combine_getmask(f, masks, I)
    if length(masks) == 2
        m1 = masks[1][I]
        m2 = masks[2][I]
        return f(m1, m2)
    else
        m1 = masks[1][I]
        m2 = _combine_getmask(f, Base.tail(masks), I)
        return f(m1, m2)
    end
end

Base.@propagate_inbounds Base.getindex(m::Indexer{M}, I::Integer...) where M <: CombinedMask = m[I]
Base.@propagate_inbounds function Base.getindex(m::Indexer{M}, I::Tuple) where M <: CombinedMask
    return _combine_getmask(m.f, m.masks, I)
end

check_constraint(m::CombinedMask, x) = check_constraint(m.masks, x)

function AxesConstraint(m::CombinedMask)
    merge_constraint(map(AxesConstraint, m.masks)...)
end

Broadcast.broadcasted(::MaskStyle, f, m1::AbstractAttenMask, m2::AbstractAttenMask) = CombinedMask(f, (m1, m2))

function Base.show(io::IO, m::CombinedMask)
    print(io, '(')
    show(io, first(m.masks))
    for mask in Base.tail(m.masks)
        print(io, ' ', m.f, ' ')
        show(io, mask)
    end
    print(io, ')')
    io
end

struct BatchedMask{M, S<:StaticInt} <: AbstractWrapperMask
    mask::M
    batch_dim::S
end

compute_batch_dim(::Tuple{}) = 0
compute_batch_dim(cs::Tuple{NDimConstraint}) = 0
compute_batch_dim(cs::Tuple{NDimConstraint, All1Constraint}) = 0
compute_batch_dim(cs::Tuple{NDimConstraint, Vararg{DimConstraint}}) = length(cs) - 1

function BatchedMask(mask)
    batch_dim = compute_batch_dim(AxesConstraint(mask))
    return BatchedMask(mask, static(batch_dim))
end

@init @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
    Adapt.adapt(to::CUDA.Adaptor, m::BatchedMask) = Indexer{typeof(m)}((mask = adapt(to, m.mask), batch_dim = m.batch_dim))
end

adapt_structure(to, x::BatchedMask) = BatchedMask(adapt(to, x.mask), x.batch_dim)

GetIndexer(m::BatchedMask) = Indexer{typeof(m)}((mask = GetIndexer(m.mask), batch_dim = m.batch_dim))

@inline function _tailtuples(I, dim)
    offset = length(I) - dim
    return ntuple(i->I[i+offset], dim)
end

Base.@propagate_inbounds function Base.getindex(m::Indexer{M}, I::Integer...) where M <: BatchedMask
    i = I[1]
    j = I[2]
    J = _tailtuples(I, m.batch_dim)
    return m.mask[(i, j, J...)]
end

batch_constraint(::Tuple{}) = ()
batch_constraint(cs::Tuple{NDimConstraint}) = (NDimConstraint(cs[1].n, true),)
batch_constraint(cs::Tuple{NDimConstraint, All1Constraint}) = (NDimConstraint(cs[1].n, true),)
function batch_constraint(cs::Tuple{NDimConstraint, Vararg{DimConstraint}})
    dcs = Base.tail(cs)
    n = length(dcs)
    return (NDimConstraint(cs[1].n, true), ntuple(i->DimConstraint(i-n-1, dcs[i].val), n)...)
end

AxesConstraint(m::BatchedMask) = batch_constraint(AxesConstraint(m.mask))

struct RepeatMask{M} <: AbstractWrapperMask
    mask::M
    num::Int
end

@init @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
    Adapt.adapt(to::CUDA.Adaptor, m::RepeatMask) = Indexer{typeof(m)}((mask = adapt(to, m.mask), num = m.num))
end

adapt_structure(to, x::RepeatMask) = RepeatMask(adapt(to, x.mask), x.num)

GetIndexer(m::RepeatMask) = Indexer{typeof(m)}((mask = GetIndexer(m.mask), num = m.num))

Base.@propagate_inbounds function Base.getindex(m::Indexer{M}, I::Integer...) where M <: RepeatMask
    dim = length(I)
    b = I[dim]
    J = Base.setindex(I,  fld1(b, m.num), dim)
    return m.mask[J]
end

@inline head(t::Tuple) = Base.reverse(Base.tail(Base.reverse(t)))

AxesConstraint(m::RepeatMask) = multiply_constraint(AxesConstraint(m.mask), m.num)

@inline multiply_constraint(::Tuple{}, _) = ()
@inline multiply_constraint(cs::Tuple{NDimConstraint}, n) = cs
@inline function multiply_constraint(cs::Tuple{NDimConstraint, All1Constraint}, n)
    c = cs[2]
    return (NDimConstraint(c.n), ntuple(c.n - c.from + 1) do i
        dim = i + c.from - 1
        DimConstraint(dim, dim != c.n ? 1 : n)
    end...)
end

@inline function multiply_constraint(cs::Tuple{NDimConstraint, Vararg{DimConstraint}}, n)
    h = head(cs)
    c = cs[end]
    return (h..., DimConstraint(c.dim, c.val * n))
end


text:
using StaticArrays

const Point2d=MVector{2,Float64}
const Point3d=MVector{3,Float64}
const Point4d=MVector{4,Float64}

Point2d() = Point2d(0.0,0.0)
Point3d() = Point3d(0.0,0.0,0.0)
Point4d() = Point4d(0.0,0.0,0.0,0.0)

"""
	norm(p::Point2d)

Return the norm of a `Point2d` object
# Example
```
julia> GL.norm(GL.Point2d(1,1))
1.4142135623730951
```
"""
function norm(p::Point2d)
	return sqrt(p[1]*p[1] + p[2]*p[2])
end

"""
	norm(p::Point3d)

Return the norm of a `Point3d` object
# Example
```
julia> GL.norm(GL.Point3d(1,1,1))
1.7320508075688772
```
"""
function norm(p::Point3d)
	return sqrt(p[1]*p[1] + p[2]*p[2] + p[3]*p[3])
end

"""
norm(p::Point4d)

Return the norm of a `Point4d` object
# Example
```
julia> GL.norm(GL.Point4d(1,1,1,1))
2.0
```
"""
function norm(p::Point4d)
	return sqrt(p[1]*p[1] + p[2]*p[2] + p[3]*p[3] + p[4]*p[4])
end

"""
	normalized(p::Point2d)

Return the unit vector parallel to `p::Point2d`.
# Example
```
julia> GL.normalized(GL.Point2d(1,1))
2-element MArray{Tuple{2},Float64,1,2}:
 0.7071067811865475
 0.7071067811865475
```
"""
function normalized(p::Point2d)
	len=norm(p)
	return Point2d(p[1] / len, p[2] / len)
end

"""
normalized(p::Point3d)

Return the unit vector parallel to `p::Point3d`.
# Example
```
julia> GL.normalized(GL.Point3d(1,1,1))
3-element MArray{Tuple{3},Float64,1,3}:
 0.5773502691896258
 0.5773502691896258
 0.5773502691896258

julia> GL.norm(GL.normalized(GL.Point3d(1,1,1)))
1.0
```
"""
function normalized(p::Point3d)
	len=norm(p)
	return Point3d(p[1] / len, p[2] / len, p[3] / len)
end

"""
normalized(p::Point4d)

Return the unit vector parallel to `p::Point4d`.
# Example
```
julia> GL.normalized(GL.Point4d(1,1,1,1))
4-element MArray{Tuple{4},Float64,1,4}:
 0.5
 0.5
 0.5
 0.5

julia> GL.norm(GL.normalized(GL.Point4d(1,1,1)))
1.0
```
"""
function normalized(p::Point4d)
	len=norm(p)
	return Point4d(p[1] / len, p[2] / len, p[3] / len, p[4] / len)
end

"""
	dropW(p::Point4d)

Divide `p` for the fourth coordinate, and return the first three ones.
# Example
```
julia> GL.dropW(GL.Point4d(2,2,2,2))
3-element MArray{Tuple{3},Float64,1,3}:
 1.0
 1.0
 1.0
"""
function dropW(p::Point4d)
	return Point3d(p[1] / p[4], p[2] / p[4], p[3] / p[4])
end