text:
Vcheckbounds(A::Operator,kr::Colon) = nothing

checkbounds(A::Operator,kr) =
    (maximum(kr) > length(A) || minimum(kr) < 1) && throw(BoundsError(A,kr))

checkbounds(A::Operator,kr::Union{Colon,InfRanges},jr::Union{Colon,InfRanges}) = nothing

checkbounds(A::Operator,kr::Union{Colon,InfRanges},jr) =
    (maximum(jr) > size(A,2) || minimum(jr) < 1) && throw(BoundsError(A,(kr,jr)))

checkbounds(A::Operator,kr,jr::Union{Colon,InfRanges}) =
    (maximum(kr) > size(A,1)  || minimum(kr) < 1 ) && throw(BoundsError(A,(kr,jr)))

checkbounds(A::Operator,kr,jr) =
    (!isempty(kr) && (maximum(kr) > size(A,1) || minimum(kr) < 1)) ||
    (!isempty(jr) && (maximum(jr) > size(A,2) || minimum(jr) < 1)) &&
    throw(BoundsError(A,(kr,jr)))

checkbounds(A::Operator,K::Block,J::Block) =
     1 ≤ first(K.n[1]) ≤ length(blocklengths(rangespace(A))) &&
     1 ≤ first(J.n[1]) ≤ length(blocklengths(domainspace(A)))

checkbounds(A::Operator,K::BlockRange{1},J::BlockRange{1}) =
    isempty(K) || isempty(J) ||
        checkbounds(A, Block(maximum(K.indices[1])), Block(maximum(J.indices[1])))

## SubOperator

struct SubOperator{T,B,I,DI,BI} <: Operator{T}
    parent::B
    indexes::I
    dims::DI
    bandwidths::BI
end

function SubOperator(A,inds,dims,lu)
    checkbounds(A,inds...)
    SubOperator{eltype(A),typeof(A),typeof(inds),
                typeof(dims),typeof(lu)}(A,inds,dims,lu)
end

# work around strange bug with bool size
SubOperator(A,inds,dims::Tuple{Bool,Bool},lu) = SubOperator(A,inds,Int.(dims),lu)

function SubOperator(A,inds::Tuple{Block,Block},lu)
    checkbounds(A,inds...)
    SubOperator(A,inds,(blocklengths(rangespace(A))[inds[1].n[1]],blocklengths(domainspace(A))[inds[2].n[1]]),lu)
end

SubOperator(A, inds::Tuple{Block,Block}) = SubOperator(A,inds,subblockbandwidths(A))
function SubOperator(A, inds::Tuple{BlockRange{1,R},BlockRange{1,R}}) where R
    checkbounds(A,inds...)
    dims = (sum(blocklengths(rangespace(A))[inds[1].indices[1]]),
            sum(blocklengths(domainspace(A))[inds[2].indices[1]]))
    SubOperator(A,inds,dims,(dims[1]-1,dims[2]-1))
end

# cannot infer ranges
SubOperator(A,inds,dims) = SubOperator(A,inds,dims,(dims[1]-1,dims[2]-1))
SubOperator(A,inds) = SubOperator(A,inds,map(length,inds))

convert(::Type{Operator{T}},SO::SubOperator) where {T} =
    SubOperator(Operator{T}(SO.parent),SO.indexes,SO.dims,SO.bandwidths)::Operator{T}

function view(A::Operator,kr::InfRanges,jr::InfRanges)
    @assert isinf(size(A,1)) && isinf(size(A,2))
    st=step(kr)
    if isbanded(A) && st==step(jr)  # Otherwise, its not a banded operator
        kr1=first(kr)
        jr1=first(jr)
        l,u=(bandwidth(A,1)+jr1-kr1)÷st,(bandwidth(A,2)+kr1-jr1)÷st
    else
        l,u=∞,∞
    end
    SubOperator(A,(kr,jr),size(A),(l,u))
end

view(V::SubOperator, kr::AbstractRange, jr::AbstractRange) =
    view(V.parent,reindex(V,parentindices(V),(kr,jr))...)

function view(A::Operator, kr::AbstractRange, jr::AbstractRange)
    st=step(kr)
    if isbanded(A) && st == step(jr)
        kr1=first(kr)
        jr1=first(jr)
        l,u=(bandwidth(A,1)+jr1-kr1)÷st,(bandwidth(A,2)+kr1-jr1)÷st
        SubOperator(A,(kr,jr),(length(kr),length(jr)),(l,u))
    else
        SubOperator(A,(kr,jr))
    end
end

function view(A::Operator,kr::UnitRange,jr::UnitRange)
    if isbanded(A)
        shft=first(kr)-first(jr)
        l,u=bandwidth(A,1)-shft,bandwidth(A,2)+shft
        SubOperator(A,(kr,jr),(length(kr),length(jr)),(l,u))
    else
        SubOperator(A,(kr,jr))
    end
end

view(A::Operator,::Colon,::Colon) = view(A,1:size(A,1),1:size(A,2))
view(A::Operator,::Colon,jr) = view(A,1:size(A,1),jr)
view(A::Operator,kr,::Colon) = view(A,kr,1:size(A,2))

view(A::Operator,K::Block,J::Block) = SubOperator(A,(K,J))
view(A::Operator,K::Block,j::Colon) = view(A,blockrows(A,K),j)
view(A::Operator,k::Colon,J::Block) = view(A,k,blockcols(A,J))
view(A::Operator, K::Block, j) = view(A,blockrows(A,Int(K)),j)
view(A::Operator, k, J::Block) = view(A,k,blockcols(A,Int(J))) #TODO: fix view
view(A::Operator,KR::BlockRange,JR::BlockRange) = SubOperator(A,(KR,JR))

view(A::Operator,k,j) = SubOperator(A,(k,j))

defaultgetindex(B::Operator,k::InfRanges, j::InfRanges) = view(B, k, j)
defaultgetindex(B::Operator,k::AbstractRange, j::InfRanges) = view(B, k, j)
defaultgetindex(B::Operator,k::InfRanges, j::AbstractRange) = view(B, k, j)

reindex(A::Operator, B::Tuple{Block,Any}, kj::Tuple{Any,Any}) =
    (reindex(rangespace(A),(B[1],), (kj[1],))[1], reindex(domainspace(A),tail(B), tail(kj))[1])
# always reindex left-to-right, so if we have only a single tuple, then
# we must be the domainspace
reindex(A::Operator, B::Tuple{Block{1}}, kj::Tuple{Any}) = reindex(domainspace(A),B,kj)

reindex(A::Operator, B::Tuple{BlockRange1,Any}, kj::Tuple{Any,Any}) =
    (reindex(rangespace(A),(B[1],), (kj[1],))[1], reindex(domainspace(A),tail(B), tail(kj))[1])
# always reindex left-to-right, so if we have only a single tuple, then
# we must be the domainspace
reindex(A::Operator, B::Tuple{BlockRange1}, kj::Tuple{Any}) =
    reindex(domainspace(A),B,kj)
# Blocks are preserved under ranges
for TYP in (:Block,:BlockRange1,:(AbstractVector{Block{1}}))
    @eval begin
        reindex(A::Operator, B::Tuple{AbstractVector{Int},Any}, kj::Tuple{$TYP,Any}) =
            (reindex(rangespace(A), (B[1],), (kj[1],))[1], reindex(domainspace(A),tail(B), tail(kj))[1])
        reindex(A::Operator, B::Tuple{AbstractVector{Int}}, kj::Tuple{$TYP}) =
            reindex(domainspace(A),B,kj)
    end
end

view(V::SubOperator,kr::UnitRange,jr::UnitRange) = view(V.parent,reindex(V,parentindices(V),(kr,jr))...)
view(V::SubOperator,K::Block,J::Block) = view(V.parent,reindex(V,parentindices(V),(K,J))...)
view(V::SubOperator,KR::BlockRange,JR::BlockRange) = view(V.parent, reindex(V,parentindices(V),(KR,JR))...)
function view(V::SubOperator,::Type{FiniteRange},jr::AbstractVector{Int})
    cs = (isbanded(V) || isblockbandedbelow(V)) ? colstop(V,maximum(jr)) : mapreduce(j->colstop(V,j),max,jr)
    view(V,1:cs,jr)
end

view(V::SubOperator, kr, jr) = view(V.parent,reindex(V,parentindices(V),(kr,jr))...)
view(V::SubOperator,kr::InfRanges,jr::InfRanges) = view(V.parent,reindex(V,parentindices(V),(kr,jr))...)

bandwidths(S::SubOperator) = S.bandwidths
function colstop(S::SubOperator{T,OP,Tuple{UnitRange{Int},UnitRange{Int}}},j::Integer) where {T,OP}
    cs = colstop(parent(S),parentindices(S)[2][j])
    kr = parentindices(S)[1]
    n = size(S,1)
    if cs < first(kr)
        0
    elseif cs ≥ last(kr)
        n
    else
        min(n,findfirst(isequal(cs),kr))
    end
end
colstart(S::SubOperator{T,OP,Tuple{UnitRange{Int},UnitRange{Int}}},j::Integer) where {T,OP} =
    max(findfirst(parentindices(S)[1],colstart(parent(S),parentindices(S)[2][j])),1)
rowstart(S::SubOperator{T,OP,Tuple{UnitRange{Int},UnitRange{Int}}},j::Integer) where {T,OP} =
    max(1,findfirst(parentindices(S)[2],rowstart(parent(S),parentindices(S)[1][j])))
rowstop(S::SubOperator{T,OP,Tuple{UnitRange{Int},UnitRange{Int}}},j::Integer) where {T,OP} =
        findfirst(parentindices(S)[2],rowstop(parent(S),parentindices(S)[1][j]))

# blocks don't change
blockcolstop(S::SubOperator{T,OP,Tuple{II,JJ}},J::Integer) where {T,OP,II<:AbstractRange{Int},JJ<:AbstractRange{Int}} =
    blockcolstop(parent(S),J)

israggedbelow(S::SubOperator) = israggedbelow(parent(S))

# since blocks don't change with indexex, neither do blockbandwidths
blockbandwidths(S::SubOperator{T,OP,Tuple{II,JJ}}) where {T,OP,II<:AbstractRange{Int},JJ<:AbstractRange{Int}} =
    blockbandwidths(parent(S))
function blockbandwidths(S::SubOperator{T,B,Tuple{BlockRange1,BlockRange1}}) where {T,B}
    KR,JR = parentindices(S)
    l,u = blockbandwidths(parent(S))
    sh = first(KR).n[1]-first(JR).n[1]
    l-sh,u+sh
end

isblockbanded(S::SubOperator{T,B,Tuple{Block,Block}}) where {T,B} = false
isbanded(S::SubOperator{T,B,Tuple{Block,Block}}) where {T,B} = isbandedblockbanded(parent(S))
bandwidths(S::SubOperator{T,B,Tuple{Block,Block}}) where {T,B} = subblockbandwidths(parent(S))
blockbandwidths(S::SubOperator{T,B,Tuple{Block,Block}}) where {T,B} = 0,0

function BandedBlockBandedMatrix(::Type{Zeros}, S::SubOperator)
    kr,jr=parentindices(S)
    KO=parent(S)
    l,u=blockbandwidths(KO)
    λ,μ=subblockbandwidths(KO)

    rt=rangespace(KO)
    dt=domainspace(KO)
    k1,j1=isempty(kr) || isempty(jr) ? (first(kr),first(jr)) :
                                        reindex(S,parentindices(S),(1,1))

    # each row/column that we differ from the the block start shifts
    # the sub block inds
    J = block(dt,j1)
    K = block(rt,k1)
    jsh=j1-blockstart(dt,J)
    ksh=k1-blockstart(rt,K)

    rows,cols = blocklengths(rangespace(S)), blocklengths(domainspace(S))

    BandedBlockBandedMatrix(Zeros{eltype(KO)}(sum(rows),sum(cols)),
                                (rows,cols), (l,u), (λ-jsh,μ+ksh))
end

function BandedBlockBandedMatrix(::Type{Zeros}, S::SubOperator{T,B,Tuple{BlockRange1,BlockRange1}}) where {T,B}
    KR,JR = parentindices(S)
    KO = parent(S)
    l,u = blockbandwidths(KO)::Tuple{Int,Int}
    λ,μ = subblockbandwidths(KO)::Tuple{Int,Int}

    rt = rangespace(KO)
    dt = domainspace(KO)
    J = first(JR)
    K = first(KR)
    bl_sh = Int(J) - Int(K)

    KBR = blocklengthrange(rt,KR)
    KJR = blocklengthrange(dt,JR)

    BandedBlockBandedMatrix(Zeros{eltype(KO)}(sum(KBR),sum(KJR)),
                                (AbstractVector{Int}(KBR),AbstractVector{Int}(KJR)), (l+bl_sh,u-bl_sh), (λ,μ))
end

function domainspace(S::SubOperator)
    P =parent(S)
    sp=domainspace(P)
    kr=parentindices(S)[2]

    SubSpace{typeof(sp),typeof(kr),domaintype(sp),rangetype(sp)}(sp,kr)
end
function rangespace(S::SubOperator)
    P =parent(S)
    sp=rangespace(P)
    kr=parentindices(S)[1]

    SubSpace{typeof(sp),typeof(kr),domaintype(sp),rangetype(sp)}(sp,kr)
end

size(V::SubOperator) = V.dims
size(V::SubOperator,k::Int) = V.dims[k]

unsafe_getindex(V::SubOperator,k::Integer,j::Integer) = V.parent[reindex(V,parentindices(V),(k,j))...]
getindex(V::SubOperator,k::Integer,j::Integer) = V.parent[reindex(V,parentindices(V),(k,j))...]
getindex(V::SubOperator,k::Integer,j::AbstractRange) = V.parent[reindex(V,parentindices(V),(k,j))...]
getindex(V::SubOperator,k::AbstractRange,j::Integer) = V.parent[reindex(V,parentindices(V),(k,j))...]
getindex(V::SubOperator,k::AbstractRange,j::AbstractRange) = V.parent[reindex(V,parentindices(V),(k,j))...]
Base.parent(S::SubOperator) = S.parent
Base.parentindices(S::SubOperator) = S.indexes

for OP in (:isblockbanded,:isblockbandedabove,:isblockbandedbelow,
                :isbandedblockbanded,:isbandedblockbandedabove,
                :isbandedblockbandedbelow)
    @eval $OP(S::SubOperator) = $OP(parent(S))
end

# TODO: These should be removed as the general purpose case will work,
# once the notion of bandedness of finite dimensional operators is made sense of

_colstops(V) = Int[max(0,colstop(V,j)) for j=1:size(V,2)]

for TYP in (:RaggedMatrix, :Matrix)
    def_TYP = Meta.parse("default_" * string(TYP))
    @eval begin
        function $TYP(V::SubOperator)
            if isinf(size(V,1)) || isinf(size(V,2))
                error("Cannot convert $V to a $TYP")
            end
            A = parent(V)
            if isbanded(A)
                $TYP(BandedMatrix(V))
            else
                $def_TYP(V)
            end
        end

        function $TYP(V::SubOperator{T,BB,NTuple{2,UnitRange{Int}}}) where {T,BB}
            if isinf(size(V,1)) || isinf(size(V,2))
                error("Cannot convert $V to a $TYP")
            end
            A = parent(V)
            if isbanded(A)
                $TYP(BandedMatrix(V))
            elseif isbandedblockbanded(A)
                N = block(rangespace(A), last(parentindices(V)[1]))
                M = block(domainspace(A), last(parentindices(V)[2]))
                B = A[Block(1):N, Block(1):M]
                RaggedMatrix(view(B, parentindices(V)...), _colstops(V))
            else
                $def_TYP(V)
            end
        end
    end
end

# fast converts to banded matrices would be based on indices, not blocks
function BandedMatrix(S::SubOperator{T,B,Tuple{BlockRange1,BlockRange1}}) where {T,B}
    A = parent(S)
    ds = domainspace(A)
    rs = rangespace(A)
    KR,JR = parentindices(S)
    BandedMatrix(view(A,
                      blockstart(rs,first(KR)):blockstop(rs,last(KR)),
                      blockstart(ds,first(JR)):blockstop(ds,last(JR))))
end

function mul_coefficients(A::SubOperator{T,B,Tuple{UnitRange{Int},UnitRange{Int}}},b) where {T,B}
    if size(A,2) == length(b)
        AbstractMatrix(A)*b
    else
        view(A,:,1:length(b))*b
    end
end


text:
"
    colwise_dot!(r::AbstractArray, a::AbstractMatrix, b::AbstractMatrix)

In-place dot product between columns of matrices `a` and `b`
"
function colwise_dot!(r::AbstractArray, a::AbstractMatrix, b::AbstractMatrix)
    n = length(r)
    @assert n == size(a, 2) == size(b, 2) && size(a, 1) == size(b, 1) "Inconsistent argument dimensions."
    for j = 1:n
        v = zero(promote_type(eltype(a), eltype(b)))
        @simd for i = 1:size(a, 1)
            @inbounds v += a[i, j]*b[i, j]
        end
        r[j] = v
    end
    return r
end

function pdf!(r, X, μ, Σ, Z)
    broadcast!(-, Z, X, μ)
    c = 1/(2π)^(length(μ)/2)*sqrt(det(Σ))
    colwise_dot!(r, Z, inv(Σ) * Z)
    broadcast!(e->c*exp(-e/2), r, r)
end

"Calculate responsibilities (or posterior probability)"
function logpdf!(r, X::AbstractMatrix, μ::AbstractVector, L::LowerTriangular,
                 Z::AbstractMatrix=similar(X))
    broadcast!(-, Z, X, μ)
    lmul!(L,Z)
    sum!(r', Z.*Z)
    r .+= size(X,1)*log2π - 2*cliplogdet(L)
    r ./= -2
end

function logpdf!(r, X::AbstractMatrix{T}, μ::AbstractVector{T}, D::Diagonal{T},
                 Z::AbstractMatrix{T}=zeros(T,0,0)) where {T<:AbstractFloat}
    m, n = size(X)
    @inbounds @simd for j in 1:n
        s = zero(T)
        for i in 1:m
            s += abs2(D[i,i] * (X[i,j] - μ[i]))
        end
        r[j] = s
    end
    r .+= m*log2π - 2*cliplogdet(D)
    r ./= -2
end

function logpdf!(r, X::AbstractMatrix{T}, μ::AbstractVector{T}, D::UniformScaling{T},
                 Z::AbstractMatrix{T}=zeros(T,0,0)) where {T<:AbstractFloat}
    m, n = size(X)
    J = D[1,1]
    logdetD = -2m*log(J)
    @inbounds @simd for j in 1:n
        s = zero(T)
        for i in 1:m
            s += abs2(J*(X[i,j] - μ[i]))
        end
        r[j] = s
    end
    r .+= m*log2π + logdetD
    r ./= -2
end

function logpdf!(r, X::AbstractMatrix, μ::AbstractVector, Σ::Hermitian, Z=similar(X))
    broadcast!(-, Z, X, μ)
    colwise_dot!(r, Z, Σ \ Z)
    r ./= -2
    r .+= -(log(2π)*length(μ) + logdet(Σ))/2
end

function logpdf!(r, X::AbstractMatrix, μ::AbstractVector, Ch::Factorization, Z=similar(X))
    broadcast!(-, Z, X, μ)
    # define aux vars
    tmp = similar(Z)
    # r .= Z.*(inv(Ch)*Z)
    ldiv!(tmp, Ch, Z)
    sum!(r', broadcast!(*, tmp, Z, tmp))
    r ./= -2
    r .+= -(log(2π)*length(μ) + logdet(Ch))/2
end

"""
    logsumexp!(logΣexp, w)

In-place row-wise calculation of `log( ∑ exp(rᵢ) )` of row `rᵢ` in the matrix `w`.
Modifies the original weight vector to `w = exp(w-offset)`.
Uses a numerically stable algorithm with offset to control for overflow and `log1p` to control for underflow.

References:
- https://arxiv.org/pdf/1412.8695.pdf eq 3.8 for p(y)
- https://discourse.julialang.org/t/fast-logsumexp/22827/7?u=baggepinnen for stable logsumexp
"""
function logsumexp!(Σ::AbstractVector{T}, w::AbstractMatrix{T};
                    tmpind = Array{CartesianIndex{2}}(undef, length(Σ)),
                    offset = similar(Σ)) where {T <: AbstractFloat}
    findmax!(offset, tmpind, w)
    w .= exp.(w .- offset)
    sumexept!(Σ, w, tmpind) # Σ = ∑wₑ-1
    Σ .= log1p.(Σ) .+ offset
end

logsumexp(w::AbstractMatrix{T})  where {T <: AbstractFloat} =
    logsumexp!(zeros(T, size(w,1)), copy(w))

function sumexept!(s, w, idx)
    n, d = size(w)
    @inbounds for i in 1:n
        l, m = idx[i][1], idx[i][2]
        s[l] = 0
        w[l, m] -= 1
        @simd for j in 1:d
            s[l] += w[l, j]
        end
        w[l, m] += 1
    end
    s
end

function logsumexp!(Σ::AbstractVector{T}, w::AbstractMatrix{T}, offset::AbstractVector{T}) where {T <: AbstractFloat}
    n, d = size(w)
    maximum!(offset, w)
    @inbounds @simd for i in 1:n
        o = offset[i]
        wr = view(w, i, :)
        Σ[i] = -1
        for e in wr
            Σ[i] += exp(e - o)
        end
        Σ[i] = log1p(Σ[i]) + o
        wr .= exp.(wr .- Σ[i])
    end
end

# generate clusters
function clusters(t, n, d=2)
    ct = cos(t)
    st = sin(t)

    x0 = randn(n) .+ 5.0
    y0 = randn(n) .* 0.3

    x = x0 .* ct .- y0 .* st
    y = x0 .* st .+ y0 .* ct
    if d <= 2
        [x y]'
    else
        vcat([x y]', rand(d-2, n))
    end
end

function initialize_kmeans(X::AbstractMatrix{T}, k::Int) where {T <: AbstractFloat}
    d, n = size(X)
    C = kmeans(X, k)
    A = assignments(C)
    [A[i] == j ? one(T) : zero(T) for i in 1:n, j in 1:k]
end

function initialize_random(X::AbstractMatrix{T}, k::Int) where {T <: AbstractFloat}
    d, n = size(X)
    Rₙₖ = rand(T, n, k)
    Rₙₖ ./ sum(Rₙₖ, dims=2)
end

function stats(X::AbstractMatrix{T}, Rₙₖ::AbstractMatrix{T}) where {T<:AbstractFloat}
    nₖ = sum(Rₙₖ, dims=1) .+ 10eps(T)
    nₖ, X*Rₙₖ./nₖ
end

function stats!(nₖ::AbstractVector{T}, μₖ::AbstractMatrix{T},
                X::AbstractMatrix{T},  Rₙₖ::AbstractMatrix{T}) where {T<:AbstractFloat}
    sum!(nₖ', Rₙₖ)
    nₖ .+= 10eps(T)
    mul!(μₖ,X,Rₙₖ)
    μₖ ./= nₖ'
end

function covariances!(Σₖ::AbstractArray{T,3}, nₖ::AbstractVector{T}, μₖ::AbstractMatrix{T},
                      X::AbstractMatrix{T}, Rₙₖ::AbstractMatrix{T}, ::Type{FullNormal};
                      covreg::Real=1e-6, aux::AbstractMatrix{T}=similar(X)) where {T<:AbstractFloat}
    d, k = size(μₖ)
    for j in 1:k
        broadcast!(-, aux, X, @view μₖ[:,j])
        Σ = view(Σₖ, :, :, j)
        Σ .= (Rₙₖ[:,j]'.*aux)*aux'/nₖ[j]
        Σ[1:(d+1):end] .+= covreg
    end
    Σₖ
end

function covariances!(Σₖ::AbstractArray{T,3}, nₖ::AbstractVector{T}, μₖ::AbstractMatrix{T},
                      X::AbstractMatrix{T}, Rₙₖ::AbstractMatrix{T}, ::Type{DiagNormal};
                      covreg::Real=1e-6, aux::AbstractMatrix{T}=similar(μₖ)) where {T<:AbstractFloat}
    d, k = size(μₖ)
    S = reshape(Σₖ, d, k)
    mul!(S, (X.*X), Rₙₖ)
    S ./= nₖ'
    broadcast!(*, aux, μₖ, (X*Rₙₖ))
    aux ./= nₖ'
    S .-= 2*aux - μₖ.^2 .- covreg
    Σₖ
end

function covariances!(Σₖ::AbstractArray{T,3}, nₖ::AbstractVector{T}, μₖ::AbstractMatrix{T},
                      X::AbstractMatrix{T}, Rₙₖ::AbstractMatrix{T}, ::Type{IsoNormal};
                      covreg::Real=1e-6, aux::AbstractMatrix{T}=similar(μₖ)) where {T<:AbstractFloat}
    d, k = size(μₖ)
    Σs = reshape(aux, d, 1, k)
    covariances!(Σs, nₖ, μₖ, X, Rₙₖ, DiagNormal, covreg=covreg)
    Σₖ .= mean(Σs, dims=1)
end

new_params(::Type{FullNormal}, ::Type{T}, d, k) where {T<:AbstractFloat} =
    (zeros(T, k), zeros(T, d, k), zeros(T, d, d, k))

new_params(::Type{DiagNormal}, ::Type{T}, d, k) where {T<:AbstractFloat} =
    (zeros(T, k), zeros(T, d, k), zeros(T, d, 1, k))

new_params(::Type{IsoNormal}, ::Type{T}, d, k) where {T<:AbstractFloat} =
    (zeros(T, k), zeros(T, d, k), zeros(T, 1, 1, k))

auxiliary(::Type{FullNormal}, X::AbstractMatrix{T}, k::Int) where {T<:AbstractFloat} =
    (zeros(T, size(X,2)), similar(X))

auxiliary(::Type{DiagNormal}, X::AbstractMatrix{T}, k::Int) where {T<:AbstractFloat} =
    (zeros(T, size(X,2)), zeros(T, size(X,1), k))

auxiliary(::Type{IsoNormal}, X::AbstractMatrix{T}, k::Int) where {T<:AbstractFloat} =
    auxiliary(DiagNormal, X, k)

"""
    initialize(X::AbstractMatrix{T}, k::Int; init=:kmenas)

Initialization of the Gaussian mixture parameters.
"""
function initialize(::Type{MV}, X::AbstractMatrix{T}, k::Int;
                    init=:kmeans) where {T<:AbstractFloat, MV<:MultivariateNormal}
    Rₙₖ = if init == :kmeans
        initialize_kmeans(X, k)
    elseif init == :random
        initialize_random(X, k)
    else
        error("Unimplemented initialization algorithm: $init")
    end

    d, n = size(X)
    nₖ, μₖ, Σₖ = new_params(MV, T, d, k)
    stats!(nₖ, μₖ, X, Rₙₖ)
    covariances!(Σₖ, nₖ, μₖ, X, Rₙₖ, MV)
    nₖ/n, μₖ, Σₖ, Rₙₖ
end

function precision(::Type{FullNormal}, Σ::AbstractMatrix)
    Ch = cholesky!(Hermitian(Σ, :L), check=false)
    !issuccess(Ch) && @debug "Cholesky factorization failed" err=Ch.info Σ
    if Ch.info<0
        error("Cholesky factorization failed ($(Ch.info)).\\
              Try to decrease the number of components or increase `covreg`: $covreg.")
    end
    LinearAlgebra.inv!(Ch.L)
end

function precision(::Type{DiagNormal}, Σ::AbstractMatrix)
    any(iszero, Σ) && error("Factorization failed. Try to decrease the number of components")
    Diagonal(1 ./ sqrt.(vec(Σ)))
end

function precision(::Type{IsoNormal}, Σ::AbstractMatrix)
    any(iszero, Σ) && error("Factorization failed. Try to decrease the number of components")
    (1/sqrt(Σ[]))I
end

function repaircov!(Σ)
    F = eigen(Σ)
    V = F.values
    # replace negative eigenvalues by zero
    V .= max.(V, 0)
    # reconstruct covariance matrix
    Σ .= F.vectors * Diagonal(V) * F.vectors'
end

cliplog(e::T) where {T<:AbstractFloat} = ifelse(abs(e)>eps(T), e, eps(T)) |> log
cliplogdet(X::AbstractMatrix) = sum(cliplog, diag(X))

ext: jl


text:
import .AbstractOps: AbstractOp, weight

abstract type AbstractTerm{W, CoeffT} end

Base.copy(pt::AbstractTerm) = typeof(pt)(copy(op_string(pt)), copy(pt.coeff))

function _show_abstract_term(io::IO, term::AbstractTerm)
    print(io, op_string(term))
    print(io, " * ")
    if term.coeff isa Real  # could use CoeffT here.
        print(io, term.coeff)
    else
        print(io, "(", term.coeff, ")")
    end
end

function Base.show(io::IO, term::AbstractTerm)
    m = length(term)
    print(io, m, "-factor", " ", typeof(term), ":\n")
    _show_abstract_term(io, term)
end

function Base.show(io::IO, ps::AbstractTerm{T,Z4Group}) where {T}
    print(io, ps.coeff, " ")
    print(io, op_string(ps))
end

# This is not used to show output at the REPL. Don't know why.
Base.show(io::IO, ::MIME{Symbol("text/plain")}, v::AbstractVector{T}) where {T <: AbstractTerm} = show(io, v)
function Base.show(io::IO, v::AbstractVector{T}) where {T <: AbstractTerm}
    print(io, length(v), "-element ", typeof(v), ":\n")
    for term in v
        _show_abstract_term(io, term)
        println(io)
    end
    return nothing
end

Base.show(m::MIME{Symbol("text/input")}, term::AbstractTerm) = show(stdout, m, term)
function Base.show(io::IO, mime::MIME{Symbol("text/input")}, term::AbstractTerm)
    print(io, strip_typeof(term), "(")
    print(io, "\"", op_string(term), "\"")
    print(io, ", ")
    show(io, mime, term.coeff)
    print(io, ")")
end

Base.one(term::AbstractTerm{W}) where {W} = typeof(term)(fill(one(W), length(term)), one(term.coeff))

function Base.isone(term::AbstractTerm)
    return all(isone, op_string(term)) && isone(term.coeff)
end

Base.:(==)(op1::AbstractTerm, op2::AbstractTerm) = op1.coeff == op2.coeff && op_string(op1) == op_string(op2)

_isless(x, y) = isless(x, y)
_isless(x::Complex, y::Complex) = isless(abs2(x), abs2(y))

function Base.isless(op1::AbstractTerm, op2::AbstractTerm)
    if op_string(op1) == op_string(op2)
        return _isless(op1.coeff, op2.coeff)
    end
    return isless(op_string(op1), op_string(op2))
end

####
#### Container interface
####

# :popat!
for func in (:length, :size, :eltype, :eachindex, :axes, :splice!, :getindex,
             :setindex!, :iterate, :pop!, :popfirst!, :lastindex, :firstindex, :first, :last)
    @eval begin
        Base.$func(ps::AbstractTerm, args...) = $func(op_string(ps), args...)
    end
end

# Remove ambiguity found by Aqua
Base.last(ps::AbstractTerm, n::Integer) = last(op_string(ps), n)
Base.first(ps::AbstractTerm, n::Integer) = first(op_string(ps), n)

## TODO: Is there a way to combine this with methods for Base functions above ?
for func in (:weight, )
    @eval begin
        $func(ps::AbstractTerm, args...) = $func(op_string(ps), args...)
    end
end

for func in (:push!, :pushfirst!, :insert!)
    @eval begin
        Base.$func(ps::AbstractTerm, args...) = ($func(op_string(ps), args...); ps)
    end
end

Base.iszero(term::AbstractTerm) = iszero(term.coeff) || any(iszero, op_string(term))

"""
    reverse!(pt::AbstractTerm)

Reverse the order of the factors in `pt` in place.
See `reverse`.
"""
Base.reverse!(pt::AbstractTerm) = (reverse!(op_string(pt)); pt)

"""
    reverse(pt::AbstractTerm)

Reverse the order of the factors in `pt`.
See `reverse!`.
"""
Base.reverse(pt::AbstractTerm) = reverse!(copy(pt))

####
#### Algebra
####

Base.:*(z::Number, term::AbstractTerm) = strip_typeof(term)(op_string(term), term.coeff * z)

Base.:*(ps::AbstractTerm, z::Number) = z * ps

Base.:*(z::Number, p::AbstractOp) = term_type(typeof(p))([p], z)
Base.:*(p::AbstractOp, z::Number) = z * p

Base.:-(p::AbstractTerm) = typeof(p)(op_string(p), -p.coeff)

Base.:/(p::AbstractOp, z::Number) = term_type(typeof(p))([p], inv(z))

Base.:/(ps::AbstractTerm, z::Number) = typeof(ps)(op_string(ps), ps.coeff / z)


text:
# Copyright © 2021 Alexander L. Hayes
# Apache 2.0 License

module TestConvert

using Test
using RelationalDatasets

@testset "Test converting propositional datasets" begin

    @test RelationalDatasets.from_vector(
        [0 1 1; 1 0 2; 2 2 0],
        [0, 0, 1],
    ) == (RelationalDatasets.RelationalDataset((
            pos=["v4(id3)."],
            neg=["v4(id1).", "v4(id2)."],
            facts=[
                "v1(id1,0).",
                "v1(id2,1).",
                "v1(id3,2).",
                "v2(id1,1).",
                "v2(id2,0).",
                "v2(id3,2).",
                "v3(id1,1).",
                "v3(id2,2).",
                "v3(id3,0).",
            ])),
        [
            "v1(+id,#varv1).",
            "v2(+id,#varv2).",
            "v3(+id,#varv3).",
            "v4(+id).",
        ]
    )

    @test RelationalDatasets.from_vector(
        [0 1 1; 1 0 2; 2 2 0],
        [0.1, 0.2, 0.3],
    ) == (RelationalDatasets.RelationalDataset((
        pos=[
            "regressionExample(v4(id1),0.1).",
            "regressionExample(v4(id2),0.2).",
            "regressionExample(v4(id3),0.3).",
        ],
        neg=[],
        facts=[
            "v1(id1,0).",
            "v1(id2,1).",
            "v1(id3,2).",
            "v2(id1,1).",
            "v2(id2,0).",
            "v2(id3,2).",
            "v3(id1,1).",
            "v3(id2,2).",
            "v3(id3,0).",
        ])),
        [
            "v1(+id,#varv1).",
            "v2(+id,#varv2).",
            "v3(+id,#varv3).",
            "v4(+id).",
        ]
    )

    @test RelationalDatasets.from_vector(
        [1 1; 0 1],
        [0, 1],
        ["a", "b", "c"],
    ) == (RelationalDatasets.RelationalDataset((
        pos=["c(id2)."],
        neg=["c(id1)."],
        facts=[
            "a(id1,1).",
            "a(id2,0).",
            "b(id1,1).",
            "b(id2,1).",
        ])),
        [
            "a(+id,#vara).",
            "b(+id,#varb).",
            "c(+id).",
        ]
    )

    @test RelationalDatasets.from_vector(
        [1 1; 0 1],
        [0.5, 1.0],
        ["a", "b", "c"],
    ) == (RelationalDatasets.RelationalDataset((
        pos=["regressionExample(c(id1),0.5).", "regressionExample(c(id2),1.0)."],
        neg=[],
        facts=[
            "a(id1,1).",
            "a(id2,0).",
            "b(id1,1).",
            "b(id2,1).",
        ])),
        [
            "a(+id,#vara).",
            "b(+id,#varb).",
            "c(+id).",
        ]
    )

    @test RelationalDatasets.from_vector(
        [0 1 1; 1 0 2; 2 2 0],
        [0, 1, 2],
    ) == (RelationalDatasets.RelationalDataset((
            pos=["v4(id1,0).", "v4(id2,1).", "v4(id3,2)."],
            neg=[],
            facts=[
                "v1(id1,0).",
                "v1(id2,1).",
                "v1(id3,2).",
                "v2(id1,1).",
                "v2(id2,0).",
                "v2(id3,2).",
                "v3(id1,1).",
                "v3(id2,2).",
                "v3(id3,0).",
            ])),
        [
            "v1(+id,#varv1).",
            "v2(+id,#varv2).",
            "v3(+id,#varv3).",
            "v4(+id,#classlabel).",
        ]
    )

end

@testset "Misaligned Classification and Regression" begin
    @test_throws AssertionError RelationalDatasets.from_vector([1 2], [1])
    @test_throws AssertionError RelationalDatasets.from_vector([1 2], [0.5])
end

@testset "Misaligned Custom Variable Names" begin
    @test_throws AssertionError RelationalDatasets.from_vector(
        [1 2],
        [1],
        ["a", "b", "c", "d"],
    )

    @test_throws AssertionError RelationalDatasets.from_vector(
        [1 2],
        [0.5],
        ["a", "b", "c", "d"]
    )
end

end # module


text:
#####################################################################
##                      Data Type Declarations                     ##
#####################################################################

struct ExtendedStateSpace{TE,T} <: AbstractStateSpace{TE}
    A::Matrix{T}
    B1::Matrix{T}
    B2::Matrix{T}
    C1::Matrix{T}
    C2::Matrix{T}
    D11::Matrix{T}
    D12::Matrix{T}
    D21::Matrix{T}
    D22::Matrix{T}
    timeevol::TE
    function ExtendedStateSpace{TE,T}(
        A,
        B1,
        B2,
        C1,
        C2,
        D11,
        D12,
        D21,
        D22,
        timeevol::TE,
    ) where {TE,T}
        nx = size(A, 1)
        nw = size(B1, 2)
        nu = size(B2, 2)
        nz = size(C1, 1)
        ny = size(C2, 1)

        if size(A, 2) != nx && nx != 0
            error("A must be square")
        elseif size(B1, 1) != nx
            error("B1 must have the same row size as A")
        elseif size(B2, 1) != nx
            error("B2 must have the same row size as A")
        elseif size(C1, 2) != nx
            error("C1 must have the same column size as A")
        elseif size(C2, 2) != nx
            error("C2 must have the same column size as A")
        elseif nw != size(D11, 2)
            error("D11 must have the same column size as B1")
        elseif nw != size(D21, 2)
            error("D21 must have the same column size as B1")
        elseif nu != size(D12, 2)
            error("D12 must have the same column size as B2")
        elseif nu != size(D22, 2)
            error("D22 must have the same column size as B2")
        elseif nz != size(D11, 1)
            error("D11 must have the same row size as C1")
        elseif nz != size(D12, 1)
            error("D12 must have the same row size as C1")
        elseif ny != size(D21, 1)
            error("D21 must have the same row size as C2")
        elseif ny != size(D22, 1)
            error("D22 must have the same row size as C2")
        end

        new{TE,T}(A, B1, B2, C1, C2, D11, D12, D21, D22, timeevol)
    end
end

function ExtendedStateSpace(
    A::AbstractArray{T},
    B1::AbstractArray,
    B2::AbstractArray,
    C1::AbstractArray,
    C2::AbstractArray,
    D11::AbstractArray = zeros(T, size(C1, 1), size(B1, 2)),
    D12::AbstractArray = zeros(T, size(C1, 1), size(B2, 2)),
    D21::AbstractArray = zeros(T, size(C2, 1), size(B1, 2)),
    D22::AbstractArray = zeros(T, size(C2, 1), size(B2, 2)),
    Ts = nothing,
) where T
        # Validate sampling time
    if (Ts isa Real) && Ts <= 0
        error("Ts must be either a positive number or nothing
                (continuous system)")
    end
    if Ts isa Real
        Ts = Discrete(Ts)
    elseif Ts === nothing
        Ts = Continuous()
    end
    return ExtendedStateSpace{typeof(Ts), T}(
        to_matrix(T, A),
        to_matrix(T, B1),
        to_matrix(T, B2),
        to_matrix(T, C1),
        to_matrix(T, C2),
        to_matrix(T, D11),
        to_matrix(T, D12),
        to_matrix(T, D21),
        to_matrix(T, D22),
        Ts,
    )
end

"""
    ss(A, B1, B2, C1, C2, D11, D12, D21, D22 [, Ts])
    ss(A, B1, B2, C1, C2; D11, D12, D21, D22 [, Ts])
Create [`ExtendedStateSpace`](@ref)
"""
function ss(
    A::AbstractArray,
    B1::AbstractArray,
    B2::AbstractArray,
    C1::AbstractArray,
    C2::AbstractArray,
    D11::AbstractArray,
    D12::AbstractArray,
    D21::AbstractArray,
    D22::AbstractArray,
    Ts = nothing,
)
    return ExtendedStateSpace(A, B1, B2, C1, C2, D11, D12, D21, D22, Ts)
end

function ss(
    A::AbstractArray{T},
    B1::AbstractArray,
    B2::AbstractArray,
    C1::AbstractArray,
    C2::AbstractArray;
    D11::AbstractArray = zeros(T, size(C1, 1), size(B1, 2)),
    D12::AbstractArray = zeros(T, size(C1, 1), size(B2, 2)),
    D21::AbstractArray = zeros(T, size(C2, 1), size(B1, 2)),
    D22::AbstractArray = zeros(T, size(C2, 1), size(B2, 2)),
    Ts = nothing,
) where T
    return ExtendedStateSpace(A, B1, B2, C1, C2, D11, D12, D21, D22, Ts)
end

function Base.promote_rule(::Type{StateSpace{TE, F1}}, ::Type{ExtendedStateSpace{TE, F2}}) where {TE, F1, F2}
    ExtendedStateSpace{TE, promote_type(F1, F2)}
end

function Base.convert(::Type{ExtendedStateSpace{TE, F2}}, s::StateSpace{TE, F1})where {TE, F1, F2}
    partition(s, 0, 0)
end

function Base.getproperty(sys::ExtendedStateSpace, s::Symbol)
    if s === :Ts
        # if !isdiscrete(sys) # NOTE this line seems to be breaking inference of isdiscrete (is there a test for this?)
        if isdiscrete(sys)
            return timeevol(sys).Ts
        else
            @warn "Getting time 0.0 for non-discrete systems is deprecated. Check `isdiscrete` before trying to access time."
            return 0.0
        end
    elseif s === :nx
        return nstates(sys)
    elseif s === :nu
        return size(sys.B2, 2)
    elseif s === :ny # TODO: now size(sys.C, 1) is not always the same as sys.ny
        return size(sys.C2, 1)
    elseif s === :nw
        return size(sys.B1, 2)
    elseif s === :nz
        return size(sys.C1, 1)
    elseif s === :B
        [sys.B1 sys.B2]
    elseif s === :C
        [sys.C1; sys.C2]
    elseif s === :D
        [sys.D11 sys.D12; sys.D21 sys.D22]
    elseif s === :zinds
        return 1:size(sys.C1, 1)
    elseif s === :yinds
        return size(sys.C1, 1) .+ (1:size(sys.C2, 1))
    elseif s === :winds
        return 1:size(sys.B1, 2)
    elseif s === :uinds
        return size(sys.B1, 2) .+ (1:size(sys.B2, 2))
    else
        return getfield(sys, s)
    end
end

ControlSystems.StateSpace(s::ExtendedStateSpace) = ss(ssdata(s)..., s.timeevol)

ssdata_e(sys::ExtendedStateSpace) = sys.A,
sys.B1,
sys.B2,
sys.C1,
sys.C2,
sys.D11,
sys.D12,
sys.D21,
sys.D22

# Funtions for number of intputs, outputs and states
ninputs(sys::ExtendedStateSpace) = size(sys.D, 2)
noutputs(sys::ExtendedStateSpace) = size(sys.D, 1)
nstates(sys::ExtendedStateSpace) = size(sys.A, 1)

#####################################################################
##                         Math Operators                          ##
#####################################################################

## EQUALITY ##
function Base.:(==)(sys1::ExtendedStateSpace, sys2::ExtendedStateSpace)
    return all(
        getfield(sys1, f) == getfield(sys2, f) for f in fieldnames(ExtendedStateSpace)
    )
end

## Approximate ##
function Base.isapprox(sys1::ExtendedStateSpace, sys2::ExtendedStateSpace)
    return all(
        getfield(sys1, f) ≈ getfield(sys2, f) for f in fieldnames(ExtendedStateSpace)
    )
end

## ADDITION ##
function Base.:+(s1::ExtendedStateSpace, s2::ExtendedStateSpace)
    timeevol = common_timeevol(s1,s2)

    A = blockdiag(s1.A, s2.A)

    B = [[s1.B1; s2.B1] blockdiag(s1.B2, s2.B2)]

    C = [[s1.C1 s2.C1];
    blockdiag(s1.C2, s2.C2)]

    D = [(s1.D11 + s2.D11) s1.D12 s2.D12;
    [s1.D21; s2.D21] blockdiag(s1.D22, s2.D22)]

    P = StateSpace(A, B, C, D, timeevol) # How to handle discrete?
    partition(P, s1.nw, s1.nz)
end

## SUBTRACTION ##
Base.:-(s1::ExtendedStateSpace, s2::ExtendedStateSpace) = +(s1, -s2)

## MULTIPLICATION ##
# NOTE: this series connects the upper parts, i.e., the performance connecmapping
# If you want to series connect the system mapping, use invert_mappings(invert_mappings(sys1)*invert_mappings(sys2))
function Base.:*(s1::ExtendedStateSpace, s2::ExtendedStateSpace)
    timeevol = common_timeevol(s1,s2)

    A = [s1.A                           s1.B1*s2.C1;
    zeros(size(s2.A,1),size(s1.A,2))      s2.A]

    B = [s1.B1*s2.D11                         s1.B2           s1.B1*s2.D12;
    s2.B1              zeros(size(s2.B2,1),size(s1.B2,2))          s2.B2]

    C = [s1.C1                       s1.D11*s2.C1;
    s1.C2                        s1.D21*s2.C1;
    zeros(size(s2.C2,1),size(s1.C2,2))         s2.C2]

    D = [s1.D11*s2.D11           s1.D12        s1.D11*s2.D12;
    s1.D21*s2.D11           s1.D22        s1.D21*s2.D12;
    s2.D21          zeros(size(s2.D22,1),size(s1.D22,2))          s2.D22        ]

    P = StateSpace(A, B, C, D, timeevol)
    partition(P, s2.nw, s1.nz)
end

function Base.:*(s1::ExtendedStateSpace, s2::Number)
    A, B1, B2, C1, C2, D11, D12, D21, D22 = ssdata_e(s1)
    ss(A, B1, B2, s2*C1, C2, s2*D11, s2*D12, D21, D22, s1.timeevol)
end

function invert_mappings(s::ExtendedStateSpace)
    A, B1, B2, C1, C2, D11, D12, D21, D22 = ssdata_e(s)
    ss(A, B2, B1, C2, C1, D22, D21, D12, D11, s.timeevol)
end

## DIVISION ##
# function Base.:/(n::Number, sys::ExtendedStateSpace)
#     # Ensure s.D is invertible
#     A, B, C, D = ssdata(sys)
#     Dinv = try
#         inv(D)
#     catch
#         error("D isn't invertible")
#     end
#     partition(ss(A - B*Dinv*C, B*Dinv, -n*Dinv*C, n*Dinv, sys.timeevol), sys.nz, sys.nw) # nz, nw reversed (inv)
# end

## NEGATION ##
function Base.:-(sys::ST) where ST <: ExtendedStateSpace
    A, B1, B2, C1, C2, D11, D12, D21, D22 = ssdata_e(sys)
    ST(A, B1, B2, -C1, -C2, -D11, -D12, -D21, -D22, sys.timeevol)
end

#####################################################################
##                       Indexing Functions                        ##
#####################################################################
Base.ndims(::ExtendedStateSpace) = 2 # NOTE: Also for SISO systems?
Base.size(sys::ExtendedStateSpace) = (noutputs(sys), ninputs(sys)) # NOTE: or just size(sys.D)
Base.size(sys::ExtendedStateSpace, d::Integer) = d <= 2 ? size(sys)[d] : 1
Base.eltype(::Type{S}) where {S<:ExtendedStateSpace} = S
ControlSystems.numeric_type(sys::ExtendedStateSpace) = eltype(sys.A)

function Base.getindex(sys::ExtendedStateSpace, inds...)
    if size(inds, 1) != 2
        error("Must specify 2 indices to index statespace model")
    end
    rows, cols = ControlSystems.index2range(inds...) # FIXME: ControlSystems.index2range(inds...)
    return ss(
        copy(sys.A),
        sys.B[:, cols],
        sys.C[rows, :],
        sys.D[rows, cols],
        sys.timeevol,
    )
end

#####################################################################
##                        Display Functions                        ##
#####################################################################

Base.print(io::IO, sys::ExtendedStateSpace) = show(io, sys)

function Base.show(io::IO, sys::ExtendedStateSpace)
    # Compose the name vectors
    println(io, typeof(sys))
    iszero(sys.A) || println(io, "A = \n", _string_mat_with_headers(sys.A))
    iszero(sys.B1) || println(io, "B1 = \n", _string_mat_with_headers(sys.B1))
    iszero(sys.B2) || println(io, "B2 = \n", _string_mat_with_headers(sys.B2))
    iszero(sys.C1) || println(io, "C1 = \n", _string_mat_with_headers(sys.C1))
    iszero(sys.C2) || println(io, "C2 = \n", _string_mat_with_headers(sys.C2))
    iszero(sys.D11) || println(io, "D11 = \n", _string_mat_with_headers(sys.D11))
    iszero(sys.D12) || println(io, "D12 = \n", _string_mat_with_headers(sys.D12))
    iszero(sys.D21) || println(io, "D21 = \n", _string_mat_with_headers(sys.D21))
    iszero(sys.D22) || println(io, "D22 = \n", _string_mat_with_headers(sys.D22))

    # Print sample time
    if isdiscrete(sys) > 0
        println(io, "Sample Time: ", sys.timeevol.Ts, " (seconds)")
    end

    # Print model type
    if nstates(sys) == 0
        print(io, "Static gain")
    elseif iscontinuous(sys)
        print(io, "Continuous-time extended state-space model")
    else
        print(io, "Discrete-time extended state-space model")
    end
end

function ControlSystems.lft(G::ExtendedStateSpace, Δ, type=:l)

    # if !(G.nu > Δ.ny && G.ny > Δ.nu)
    #     error("Must have G.nu > Δ.ny and G.ny > Δ.nu for lower/upper lft")
    # end

    if type === :l
        lft(ss(G), Δ, :l)
    else
        error("Invalid type of lft ($type), specify type=:l")
    end
end

macro sizecompat(a,b)
    n1 = string(a)
    n2 = string(b)
    quote
        s1 = size($(esc(a)), 2)
        s2 = size($(esc(b)), 1)
        if !(s1 == s2)
            throw(ArgumentError(string("size(", $n1, ", 2) and size(", $n2, ", 1) must be compatible, got ($(s1), $(s2))")))
        end
    end
end

function ControlSystems.feedback(s1::ExtendedStateSpace, s2::ExtendedStateSpace;
    Wperm=:, Zperm=:)

    timeevol = common_timeevol(s1,s2)
    s1_B1 = s1.B1
    s1_B2 = s1.B2
    s1_C1 = s1.C1
    s1_C2 = s1.C2
    s1_D11 = s1.D11
    s1_D12 = s1.D12
    s1_D21 = s1.D21
    s1_D22 = s1.D22

    # s2_B1 = s2.B2 # These are reversed
    # s2_B2 = s2.B1
    # s2_C1 = s2.C2
    # s2_C2 = s2.C1
    # s2_D11 = s2.D22
    # s2_D12 = s2.D21
    # s2_D21 = s2.D12
    # s2_D22 = s2.D11

    s2_B1 = s2.B1
    s2_B2 = s2.B2
    s2_C1 = s2.C1
    s2_C2 = s2.C2
    s2_D11 = s2.D11
    s2_D12 = s2.D12
    s2_D21 = s2.D21
    s2_D22 = s2.D22

    @sizecompat s1_B2 s2_D22
    @sizecompat s1_B2 s2_C2
    @sizecompat s2_D22 s1_D21
    @sizecompat s2_D22 s1_C2
    @sizecompat s1_D12 s2_C2
    @sizecompat s2_D12 s1_D22
    @sizecompat s1_D22 s2_C2

    if iszero(s1_D22) || iszero(s2_D22)
        A = [s1.A + s1_B2*s2_D22*s1_C2        s1_B2*s2_C2;
                s2_B2*s1_C2            s2.A + s2_B2*s1_D22*s2_C2]

        B1 = [
            s1_B1 + s1_B2*s2_D22*s1_D21
                    s2_B2*s1_D21
        ]
        B2  = [
            s1_B2*s2_D21
            s2_B1 + s2_B2*s1_D22*s2_D21
        ]
        C1 = [s1_C1+s1_D12*s2_D22*s1_C2        s1_D12*s2_C2]
        C2 = [s2_D12*s1_C2           s2_C1+s2_D12*s1_D22*s2_C2]
        D11 = s1_D11 + s1_D12*s2_D22*s1_D21
        D12 = s1_D12*s2_D21
        D21 = s2_D12*s1_D21
        D22 = s2_D11 + s2_D12*s1_D22*s2_D21
    else
        R1 = try
            inv(I - s2_D22*s1_D22)
        catch
            error("Ill-posed feedback interconnection,  I - s2_D22*s1_D22 or I - s2_D22*s1_D22 not invertible")
        end

        R2 = try
            inv(I - s1_D22*s2_D22)
        catch
            error("Ill-posed feedback interconnection,  I - s2_D22*s1_D22 or I - s2_D22*s1_D22 not invertible")
        end

        A = [s1.A + s1_B2*R1*s2_D22*s1_C2        s1_B2*R1*s2_C2;
                s2_B2*R2*s1_C2            s2.A + s2_B2*R2*s1_D22*s2_C2]

        B1 = [s1_B1+s1_B2*R1*s2_D22*s1_D21;
                    s2_B2*R2*s1_D21]
        B2 = [s1_B2*R1*s2_D21; s2_B1 + s2_B2*R2*s1_D22*s2_D21]
        C1 = [s1_C1 + s1_D12*R1*s2_D22*s1_C2        s1_D12*R1*s2_C2]
        C2 = [s2_D12*R2*s1_C2           s2_C1+s2_D12*R2*s1_D22*s2_C2]
        D11 = [s1_D11 + s1_D12*R1*s2_D22*s1_D21]
        D12 = s1_D12*R1*s2_D21
        D21 = s2_D12*R2*s1_D21
        D22 = s2_D11 + s2_D12*R2*s1_D22*s2_D21
    end

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

ExtendedStateSpace(s::AbstractStateSpace) = ss(s.A, I(s.nx), s.B, I(s.nx), s.C; D22=s.D, Ts = s.timeevol)

"""
    partition(P::AbstractStateSpace; u, y, w=!u, z=!y)

Partition `P` into an [`ExtendedStateSpace`](@ref).
"""
function partition(P::AbstractStateSpace; u=nothing, y=nothing,
    w = nothing,
    z = nothing
)
    if w === nothing
        w = setdiff(1:P.nu, u)
    end
    if z === nothing
        z = setdiff(1:P.ny, y)
    end
    if u === nothing
        u = setdiff(1:P.nu, w)
    end
    if y === nothing
        y = setdiff(1:P.ny, z)
    end
    u = vcat(u)
    y = vcat(y)
    z = vcat(z)
    w = vcat(w)
    ss(P.A, P.B[:, w], P.B[:, u], P.C[z, :], P.C[y, :],
    P.D[z, w], P.D[z, u], P.D[y, w], P.D[y, u], P.timeevol)
end

"""
    partition(P::AbstractStateSpace, nw::Int, nz::Int)
"""
partition(P::AbstractStateSpace, nu1::Int, ny1::Int) = partition(P, w=1:nu1, z=1:ny1)

"""
    system_mapping(P::ExtendedStateSpace)

Return the system from u -> y
"""
function system_mapping(P::ExtendedStateSpace)
    ss(P.A, P.B2, P.C2, P.D22, P.timeevol)
end

"""
    performance_mapping(P::ExtendedStateSpace)

Return the system from w -> z
"""
function performance_mapping(P::ExtendedStateSpace)
    ss(P.A, P.B1, P.C1, P.D11, P.timeevol)
end

"""
    noise_mapping(P::ExtendedStateSpace)

Return the system from w -> y
"""
function noise_mapping(P::ExtendedStateSpace)
    ss(P.A, P.B1, P.C2, P.D12, P.timeevol)
end


text:
@discord_object struct EmbedFooter
    text::String
    icon_url::String
    proxy_icon_url::String
end


text:
# In this tutorial, we will learn
#  - How to solve nonlinear multi-field PDEs in Gridap
#  - How to build FE spaces whose functions have zero mean value
#
# ## Problem statement
#
# The goal of this last tutorial is to solve a nonlinear multi-field PDE. As a model problem, we consider a well known benchmark in computational fluid dynamics, the lid-driven cavity for the incompressible Navier-Stokes equations. Formally, the PDE we want to solve is: find the velocity vector $u$ and the pressure $p$ such that
#
# ```math
# \left\lbrace
# \begin{aligned}
# -\Delta u + \mathit{Re}\ (u\cdot \nabla)\ u + \nabla p = 0 &\text{ in }\Omega,\\
# \nabla\cdot u = 0 &\text{ in } \Omega,\\
# u = g &\text{ on } \partial\Omega,
# \end{aligned}
# \right.
# ```
#
# where the computational domain is the unit square $\Omega \doteq (0,1)^d$, $d=2$, $\mathit{Re}$ is the Reynolds number (here, we take $\mathit{Re}=10$), and $(w \cdot \nabla)\ u = (\nabla u)^t w$  is the well known convection operator. In this example, the driving force is the Dirichlet boundary velocity $g$, which is a non-zero horizontal velocity with a value of $g = (1,0)^t$ on the top side of the cavity, namely the boundary $(0,1)\times\{1\}$, and $g=0$ elsewhere on $\partial\Omega$. Since we impose Dirichlet boundary conditions on the entire boundary $\partial\Omega$, the mean value of the pressure is constrained to zero in order have a well posed problem,
#
# ```math
# \int_\Omega q \ {\rm d}\Omega = 0.
# ```
#
# ## Numerical Scheme
#
# In order to approximate this problem we chose a formulation based on inf-sub stable $Q_k/P_{k-1}$ elements with continuous velocities and discontinuous pressures (see, e.g., [1] for specific details). The interpolation spaces are defined as follows.  The velocity interpolation space is
#
# ```math
# V \doteq \{ v \in [C^0(\Omega)]^d:\ v|_T\in [Q_k(T)]^d \text{ for all } T\in\mathcal{T} \},
# ```
# where $T$ denotes an arbitrary cell of the FE mesh $\mathcal{T}$, and $Q_k(T)$ is the local polynomial space in cell $T$ defined as the multi-variate polynomials in $T$ of order less or equal to $k$ in each spatial coordinate. Note that, this is the usual continuous vector-valued Lagrangian FE space of order $k$ defined on a mesh of quadrilaterals or hexahedra.  On the other hand, the space for the pressure is
#
# ```math
# \begin{aligned}
# Q_0 &\doteq \{ q \in Q: \  \int_\Omega q \ {\rm d}\Omega = 0\}, \text{ with}\\
# Q &\doteq \{ q \in L^2(\Omega):\ q|_T\in P_{k-1}(T) \text{ for all } T\in\mathcal{T}\},
# \end{aligned}
# ```
# where $P_{k-1}(T)$ is the polynomial space of multi-variate polynomials in $T$ of degree less or equal to $k-1$. Note that functions in $Q_0$ are strongly constrained to have zero mean value. This is achieved in the code by removing one degree of freedom from the (unconstrained) interpolation space $Q$ and  adding a constant to the computed pressure so that the resulting function has zero mean value.
#
# The weak form associated to these interpolation spaces reads: find $(u,p)\in U_g \times Q_0$ such that $[r(u,p)](v,q)=0$ for all $(v,q)\in V_0 \times Q_0$
# where $U_g$ and $V_0$ are the set of functions in $V$ fulfilling the Dirichlet boundary condition $g$ and $0$  on $\partial\Omega$ respectively. The weak residual $r$ evaluated at a given pair $(u,p)$ is the linear form defined as
#
# ```math
# [r(u,p)](v,q) \doteq a((u,p),(v,q))+ [c(u)](v),
# ```
# with
# ```math
# \begin{aligned}
# a((u,p),(v,q)) &\doteq \int_{\Omega} \nabla v \cdot \nabla u \ {\rm d}\Omega - \int_{\Omega} (\nabla\cdot v) \ p \ {\rm d}\Omega + \int_{\Omega} q \ (\nabla \cdot u) \ {\rm d}\Omega,\\
# [c(u)](v) &\doteq \int_{\Omega} v 	\cdot \left( (u\cdot\nabla)\ u \right)\ {\rm d}\Omega.\\
# \end{aligned}
# ```
# Note that the bilinear form $a$ is associated with the linear part of the PDE, whereas $c$ is the contribution to the residual resulting from the convective term.
#
# In order to solve this nonlinear weak equation with a Newton-Raphson method, one needs to compute the Jacobian associated with the residual $r$. In this case, the Jacobian $j$ evaluated at a pair $(u,p)$ is the bilinear form defined as
#
# ```math
# [j(u,p)]((\delta u, \delta p),(v,q)) \doteq a((\delta u,\delta p),(v,q))  + [{\rm d}c(u)](\delta u,v),
# ```
# where ${\rm d}c$ results from the linearization of the convective term, namely
# ```math
# [{\rm d}c(u)](\delta u,v) \doteq \int_{\Omega} v \cdot \left( (u\cdot\nabla)\ \delta u \right) \ {\rm d}\Omega + \int_{\Omega} v \cdot \left( (\delta u\cdot\nabla)\ u \right)  \ {\rm d}\Omega.
# ```
# The implementation of this numerical scheme is done in Gridap by combining the concepts previously seen for single-field nonlinear PDEs  and linear multi-field problems.
#
# ## Discrete model
#
# We start with the discretization of the computational domain. We consider a $100\times100$ Cartesian mesh of the unit square.

using Gridap
n = 100
domain = (0,1,0,1)
partition = (n,n)
model = CartesianDiscreteModel(domain,partition)

# For convenience, we create two new boundary tags,  namely `"diri1"` and `"diri0"`, one for the top side of the square (where the velocity is non-zero), and another for the rest of the boundary (where the velocity is zero).

labels = get_face_labeling(model)
add_tag_from_tags!(labels,"diri1",[6,])
add_tag_from_tags!(labels,"diri0",[1,2,3,4,5,7,8])

# ## FE spaces
#
# For the velocities, we need to create a conventional vector-valued continuous Lagrangian FE space. In this example, we select a second order interpolation.

D = 2
order = 2
reffeᵤ = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
V = TestFESpace(model,reffeᵤ,conformity=:H1,labels=labels,dirichlet_tags=["diri0","diri1"])

# The interpolation space for the pressure is built as follows

reffeₚ = ReferenceFE(lagrangian,Float64,order-1;space=:P)
Q = TestFESpace(model,reffeₚ,conformity=:L2,constraint=:zeromean)

# With the options `:Lagrangian`, `space=:P`, `valuetype=Float64`, and `order=order-1`, we select the local polynomial space $P_{k-1}(T)$ on the cells $T\in\mathcal{T}$. With the symbol `space=:P` we specifically chose a local Lagrangian interpolation of type "P". Without using `space=:P`, would lead to a local Lagrangian of type "Q" since this is the default for quadrilateral or hexahedral elements. On the other hand, `constraint=:zeromean` leads to a FE space, whose functions are constrained to have mean value equal to zero, which is just what we need for the pressure space. With these objects, we build the trial multi-field FE spaces

uD0 = VectorValue(0,0)
uD1 = VectorValue(1,0)
U = TrialFESpace(V,[uD0,uD1])
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V, Q])
X = MultiFieldFESpace([U, P])

# ## Triangulation and integration quadrature
#
# From the discrete model we can define the triangulation and integration measure

degree = order
Ωₕ = Triangulation(model)
dΩ = Measure(Ωₕ,degree)

# ## Nonlinear weak form
#
# The different terms of the nonlinear weak form for this example are defined following an approach similar to the one discussed for the $p$-Laplacian equation, but this time using the notation for multi-field problems.

const Re = 10.0
conv(u,∇u) = Re*(∇u')⋅u
dconv(du,∇du,u,∇u) = conv(u,∇du)+conv(du,∇u)

# The bilinear form reads
a((u,p),(v,q)) = ∫( ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u) )dΩ

# The nonlinear term and its Jacobian are given by
c(u,v) = ∫( v⊙(conv∘(u,∇(u))) )dΩ
dc(u,du,v) = ∫( v⊙(dconv∘(du,∇(du),u,∇(u))) )dΩ

# Finally, the Navier-Stokes weak form residual and Jacobian can be defined as
res((u,p),(v,q)) = a((u,p),(v,q)) + c(u,v)
jac((u,p),(du,dp),(v,q)) = a((du,dp),(v,q)) + dc(u,du,v)

# With the functions `res`, and `jac` representing the weak residual and the Jacobian, we build the nonlinear FE problem:
op = FEOperator(res,jac,X,Y)

# ## Nonlinear solver phase
#
# To finally solve the problem, we consider the same nonlinear solver as previously considered for the  $p$-Laplacian equation.

using LineSearches: BackTracking
nls = NLSolver(
  show_trace=true, method=:newton, linesearch=BackTracking())
solver = FESolver(nls)

# In this example, we solve the problem without providing an initial guess (a default one equal to zero will be generated internally)

uh, ph = solve(solver,op)

# Finally, we write the results for visualization (see next figure).

writevtk(Ωₕ,"ins-results",cellfields=["uh"=>uh,"ph"=>ph])

# ![](../assets/inc_navier_stokes/ins_solution.png)
#
#  ## References
#
#  [1] H. C. Elman, D. J. Silvester, and A. J. Wathen. *Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics*. Oxford University Press, 2005.


text:
"""
    function Triangulate.triangulate(triin::TriangulateIO, maxarea, minangle=20)

Convenience routine to avoid writing out `@sprintf` each time. Returns a `TriangulateIO` object.
"""
function Triangulate.triangulate(triin::TriangulateIO, maxarea, minangle=20)
    angle = @sprintf("%.15f", minangle)
    area  = @sprintf("%.15f", maxarea)
    triout, _ = triangulate("pa$(area)q$(angle)Q", triin)
    return triout
end

"""
    function triangulateIO_to_VXYEToV(triout::TriangulateIO)

Computes `VX`,`VY`,`EToV` from a `TriangulateIO` object.
"""
function triangulateIO_to_VXYEToV(triout::TriangulateIO)
    VX,VY = (triout.pointlist[i,:] for i = 1:size(triout.pointlist,1))
    EToV = permutedims(triout.trianglelist)
    Base.swapcols!(EToV,2,3) # to match MeshData ordering
    return (VX, VY), Matrix{Int}(EToV)
end

"""
    function get_boundary_face_labels(triout::TriangulateIO, md::MeshData{2})

Find Triangle segment labels of boundary faces. Returns two arguments:
- `boundary_face_tags`: tags of faces on the boundary
- `boundary_faces`: list of faces on the boundary of the domain
"""
function get_boundary_face_labels(triout::TriangulateIO, rd::RefElemData{2, Tri}, md::MeshData{2})
    segmentlist = sort(triout.segmentlist,dims=1)
    boundary_faces = findall(vec(md.FToF) .== 1:length(md.FToF))
    boundary_face_tags = zeros(Int,length(boundary_faces))
    for (f,boundary_face) in enumerate(boundary_faces)
        element = (boundary_face - 1) ÷ rd.Nfaces + 1
        face    = (boundary_face - 1) % rd.Nfaces + 1
        vertices_on_face = sort(md.EToV[element,rd.fv[face]])
        tag_id = findfirst(c->view(segmentlist,:,c)==vertices_on_face,axes(segmentlist,2))
        boundary_face_tags[f] = triout.segmentmarkerlist[tag_id]
    end
    return boundary_face_tags, boundary_faces
end

"""
    function get_node_boundary_tags(triout::TriangulateIO,md::MeshData{2},rd::RefElemData{2,Tri})

Computes `node_tags` = `Nfp` x `Nfaces * num_elements` array where each entry is a Triangulate.jl tag number.
"""
function get_node_boundary_tags(triout::TriangulateIO, rd::RefElemData{2, Tri}, md::MeshData{2})
    boundary_face_tags, boundary_faces = get_boundary_face_labels(triout, rd, md)
    node_tags = zeros(Int, size(md.xf, 1) ÷ rd.Nfaces, md.K * rd.Nfaces) # make Nfp x Nfaces*num_elements
    for (i, boundary_face) in enumerate(boundary_faces)
        node_tags[:, boundary_face] .= boundary_face_tags[i]
    end
    node_tags = reshape(node_tags, size(md.xf)...)
end

function tag_boundary_faces(triout::TriangulateIO,
                            rd::RefElemData{2, Tri}, md::MeshData{2},
                            boundary_list::Dict{Symbol, Int})

    boundary_face_ids_list = _tag_boundary_faces(triout,rd,md,boundary_list)
    return Dict(Pair.(keys(boundary_list), boundary_face_ids_list))
end

"""
    function tag_boundary_faces(triout::TriangulateIO,
                                rd::RefElemData{2,Tri}, md::MeshData{2},
                                boundary_list::Union{NamedTuple,Dict{Symbol,Int}})

Here, `boundary_list` is a `Dict` (or `NamedTuple`) whose values are the boundary tags for a
`TriangulateIO` mesh format. The output is a `Dict` or `NamedTuple` with keys given by
`boundary_list` and `values` equal to vectors of faces on that given boundary.

Example usage:
```julia
julia> using Triangulate, StartUpDG
julia> triout = scramjet()
julia> rd = RefElemData(Tri(),N=1)
julia> md = MeshData(triangulateIO_to_VXYEToV(triout)...,rd)
julia> tag_boundary_faces(triout,rd,md, Dict(:wall=>1, :inflow=>2, :outflow=>3))
```
"""
function tag_boundary_faces(triout::TriangulateIO,
                            rd::RefElemData{2, Tri}, md::MeshData{2},
                            boundary_list::NamedTuple)

    boundary_face_ids_list = _tag_boundary_faces(triout,rd,md,boundary_list)
    return NamedTuple(Pair.(keys(boundary_list), boundary_face_ids_list))
end

# this version works for both boundary_list::Union{NamedTuple, Dict{Symbol,Int}}
function _tag_boundary_faces(triout::TriangulateIO,
                            rd::RefElemData{2, Tri}, md::MeshData{2},
                            boundary_list)
    boundary_face_tags, boundary_faces = get_boundary_face_labels(triout, rd, md)
    boundary_face_ids_list = Vector{Int}[]
    for boundary_face_flag in values(boundary_list)
        push!(boundary_face_ids_list, boundary_faces[findall(@. boundary_face_tags == boundary_face_flag)])
    end
    return boundary_face_ids_list
end

"""
    function refine(triout, h, href = h/2)

Refinement of a previous mesh given the current mesh size `h`. Preserves boundary/volume tags.
"""
function refine(triout, h, href = h/2)
    angle = @sprintf("%.15f",20)
    area  = @sprintf("%.15f",href^2)
    triout2,_ = triangulate("rpa$(area)q$(angle)Q", triout)
    return triout2
end

VertexMeshPlotter(triout::TriangulateIO) = VertexMeshPlotter(triangulateIO_to_VXYEToV(triout)..., face_vertices(Tri()))

"""
    BoundaryTagPlotter(triout::TriangulateIO)

Plot recipe to visualize boundary tags by color. Usage: `plot(BoundaryTagPlotter(triout))`
"""
struct BoundaryTagPlotter
    triout::TriangulateIO
end

RecipesBase.@recipe function f(m::BoundaryTagPlotter)
    triout = m.triout
    tags = unique(triout.segmentmarkerlist)
    num_colors = length(tags)
    colors = distinguishable_colors(num_colors)
    xseg = zeros(2, size(triout.segmentlist,2))
    yseg = zeros(2, size(triout.segmentlist,2))
    segcolor = HSV{Float32}[]
    for (col,segment) in enumerate(eachcol(triout.segmentlist))
        xseg[:,col] .= triout.pointlist[1,segment]
        yseg[:,col] .= triout.pointlist[2,segment]
        push!(segcolor, colors[triout.segmentmarkerlist[col]])
    end
    for i = 1:num_colors
        color_ids = findall(triout.segmentmarkerlist .== tags[i])

        # NaN separators for distinct lines
        x_i = vec([xseg[:, color_ids]; fill(NaN, length(color_ids))'])
        y_i = vec([yseg[:, color_ids]; fill(NaN, length(color_ids))'])

        @series begin
            marker --> :circle
            seriescolor --> permutedims(segcolor[color_ids]),
            ratio --> 1
            label --> string(tags[i])
            x_i,y_i
        end
    end
    # we already added all shapes in @series so we don't want to return a series
    # here (we are returning an empty series which is not added to the legend.)
    primary := false
    ()
end


text:
function calc_total_risk(data)
    risk = 0.0
    risk += calc_bus_risk(data)
    risk += calc_gen_risk(data)
    risk += calc_branch_risk(data)
    risk += calc_load_risk(data)

    return risk
end

function calc_bus_risk(data)
    calc_risk(data,"bus")
end

function calc_gen_risk(data)
    calc_risk(data,"gen")
end

function calc_branch_risk(data)
    calc_risk(data,"branch")
end

function calc_load_risk(data)
    calc_risk(data,"load")
end

function calc_risk(data, comp_type)
    risk = 0.0
    comp_status = _PM.pm_component_status[comp_type]
    for (comp_id, comp) in get(data,comp_type, Dict())
        if comp_type=="bus"
            risk += get(comp, "power_risk",0.0)* (comp[comp_status]==4 ? 0 : 1)
        else
            risk += get(comp, "power_risk",0.0)*comp[comp_status]
        end
        risk += get(comp, "base_risk", 0.0)
    end

    return risk
end

function calc_load(data)
    return sum(load["pd"]*load["status"] for (load_id, load) in data["load"])
end

"disable a bus that is disconnected from generators and branches"
function disable_isolated_buses!(network)
    connected_bus = Dict(id => false for (id,bus) in network["bus"])
    for (comp_id, comp) in network["gen"]
        comp_status = _PM.pm_component_status["gen"]
        if comp[comp_status]!= _PM.pm_component_status_inactive["gen"]
            connected_bus["$(comp["gen_bus"])"] = true
        end
    end
    for (comp_id, comp) in network["branch"]
        comp_status = _PM.pm_component_status["branch"]
        if comp[comp_status]!= _PM.pm_component_status_inactive["branch"]
            connected_bus["$(comp["f_bus"])"] = true
            connected_bus["$(comp["t_bus"])"] = true
        end
    end

    for (id, connected) in connected_bus
        if !connected
            network["bus"][id]["bus_type"]=4
        end
    end
    return network
end


text:
tomlpath = joinpath(@__DIR__, "hbv_config.toml")
config = Wflow.Config(tomlpath)

model = Wflow.initialize_hbv_model(config)

model = Wflow.update(model)

# test if the first timestep was written to the CSV file
flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
@testset "CSV output" begin
    csv = CSV.File(model.writer.csv_path)
    row = csv[1]

    @test row.time == DateTime("2000-01-01T00:00:00")
    @test row.Q ≈ 517.9761666607884
    @test row.temp_bycoord ≈ 2.965437173843384
    @test row.temp_byindex ≈ 1.1716821193695068
    @test row.Q_1 ≈ 501.6252231440609
    @test row.perc_33 ≈ 2.308000087738037
    @test row.perc_34 ≈ 1.8980000019073486
    @test row.perc_35 ≈ 2.7100000381469727
    @test row.perc_36 ≈ 3.818000078201294
    @test row.perc_37 ≈ 2.1440000534057617
end

@testset "first timestep" begin
    hbv = model.vertical
    @test hbv.tt[1] ≈ 0.0
    @test model.clock.iteration == 2
    @test hbv.altitude[1] == 484.2149963378906
    @test hbv.soilmoisture[1] == 134.35299682617188
    @test hbv.runoff[1] == 7.406898120121746
    @test hbv.soilevap[1] == 0.0
    @test hbv.snow[1] == 0.0
end

# run the second timestep
model = Wflow.update(model)

@testset "second timestep" begin
    hbv = model.vertical
    @test hbv.altitude[1] == 484.2149963378906
    @test hbv.soilmoisture[1] == 134.35299682617188
    @test hbv.runoff[1] == 4.353345933656929
    @test hbv.soilevap[1] == 0.0
    @test hbv.snow[1] == 0.0
end

@testset "overland flow" begin
    q = model.lateral.land.q_av
    @test sum(q) ≈ 3278.9423039643552
    @test q[500] ≈ 0.2348646912841589
    @test q[1566] ≈ 0.030463805623037462
    @test q[network.land.order[end]] ≈ 0.08683128190194403
end

@testset "river flow" begin
    q = model.lateral.river.q_av
    @test sum(q) ≈ 51941.417623244706
    @test q[500] ≈ 5.6934140250393055
    @test q[88] ≈ 8.839325386868788
    @test q[network.river.order[end]] ≈ 2.913948136664789
end

benchmark = @benchmark Wflow.update(model)
trialmin = BenchmarkTools.minimum(benchmark)

println("HBV Model update")
print_benchmark(trialmin)

Wflow.close_files(model, delete_output = false)


text:
#Type to hold Twitter API credentials

struct OAuth2 <: AbstractOAuth
    bearer_token::String
end

oauth_version_string(::OAuth2) = "2.0"

function oauth_header(oauth::OAuth2, httpmethod, baseurl, options; kwargs...)
    "Bearer $(oauth.bearer_token)"
end


text:
function prepare_path( s )
    if isempty(s)
        return "."
    else
        if isabspath( s )
            return s
        else
            if startswith( s , "@")
                return Base.load_path_expand(s)
            else
                return abspath( joinpath( @__DIR__, "..", s ) )
            end
        end
    end
end

max_stars_repo_name: manuelbb-upb/MorbitWrapper

max_issues_repo_name: manuelbb-upb/MorbitWrapper

max_forks_repo_name: manuelbb-upb/MorbitWrapper


text:
using CollectiveSpins

a = 0.5
b = 0.7
N = 500

@time CollectiveSpins.effective_interaction.cubiclattice_orthogonal(a, 4)
(om, ga) = @time CollectiveSpins.effective_interaction.cubiclattice_orthogonal(a, N)

@time CollectiveSpins.effective_interaction.hexagonallattice3d_orthogonal(a, b, 4)
(om, ga) = @time CollectiveSpins.effective_interaction.hexagonallattice3d_orthogonal(a, b, N)


text:
using CMPUtils
using Test
recodeaudio2(".")
@testset "CMPUtils.jl" begin
    #recodeimage(".",20)
    #recodeaudio("tmp",10)
    #recodetext(".",20)
    #unit_standardize(rand(10) .*10)
end


text:
using Test
using FlightMechanicsSimulator

@testset "ConstantInput" begin
    # Constant input
    ci_one = ConstantInput(1.0)
    ci_minus_one = ConstantInput(-1.0)

    @test isapprox(get_value(ci_one, 10.), 1.0)
    @test isapprox(get_value(ci_minus_one, 10.), -1.0)
    @test isapprox(get_value(vcat(ci_one, ci_minus_one), 10.), 0.0, atol=1e-15)
end

@testset "StepInput" begin
    # StepInput
    si1 = StepInput(1.0, 5.0, 10.0)
    si2 = StepInput(-1.0, 10.0, 15.0)

    @test isapprox(get_value(si1, 0.0), 0.0)
    @test isapprox(get_value(si2, 0.0), 0.0)
    @test isapprox(get_value(si1, 5.0), 1.0)
    @test isapprox(get_value(si1, 10.0), 0.0)
    @test isapprox(get_value(si2, 10.0), -1.0)
    @test isapprox(get_value(si2, 15.0), 0.0)

    @test isapprox(get_value(vcat(si1, si2), 8), 1.0)
    @test isapprox(get_value(vcat(si1, si2), 10), -1.0)
    @test isapprox(get_value(vcat(si1, si2), 12), -1.0)
end

@testset "DoubletInput" begin
    # Doublet input
    d = DoubletInput(2, 0, 10)
    id = DoubletInput(-2, 0, 10)

    @test isapprox(get_value(d, 2.5), 1)
    @test isapprox(get_value(d, 5), -1)
    @test isapprox(get_value(d, 7.5), -1)
    @test isapprox(get_value(d, 10.), 0)

    @test isapprox(get_value(id, 2.5), -1)
    @test isapprox(get_value(id, 5), 1)
    @test isapprox(get_value(id, 7.5), 1)
    @test isapprox(get_value(id, 10.), 0)

    @test isapprox(get_value(vcat(d, id), 2.5), 0)
    @test isapprox(get_value(vcat(d, id), 5), 0)
    @test isapprox(get_value(vcat(d, id), 7.5), 0)
    @test isapprox(get_value(vcat(d, id), 10.), 0)
end

@testset "RampInput" begin
    # Ramp input
    r1 = RampInput(10, 0, 10)
    r2 = RampInput(-10, 10, 15)

    @test isapprox(get_value(r1, 0.0), 0.0)
    @test isapprox(get_value(r1, 10.0), 0.0)
    @test isapprox(get_value(r1, 5.0), 5.0)
    @test isapprox(get_value(r2, 10.0), 0.0)
    @test isapprox(get_value(r2, 15.0), 0.0)
    @test isapprox(get_value(r2, 12.5), -5.0)

    @test isapprox(get_value(vcat(r1, r2, StepInput(10, 10, 15)), 10), 10.0)
    @test isapprox(get_value(vcat(r1, r2, StepInput(10, 10, 15)), 5), 5.0)
    @test isapprox(get_value(vcat(r1, r2, StepInput(10, 10, 15)), 12.5), 5.0)
end

@testset "Sinusoidal" begin
    # Sinusoidal input
    s = SinusoidalInput(2, 1, 0, 0, Inf)

    @test isapprox(get_value(s, 0), 0)
    @test isapprox(get_value(s, 1), 0, atol=1e-15)
    @test isapprox(get_value(s, 0.5), 0, atol=1e-15)
    @test isapprox(get_value(s, 0.25), 1)
    @test isapprox(get_value(s, 0.75), -1)

    sA = SinusoidalInput(20, 1, 0, 0, Inf)

    @test isapprox(get_value(sA, 0.25), 10)
    @test isapprox(get_value(sA, 0.75), -10)

    sf = SinusoidalInput(2, 2, 0, 0, Inf)

    @test isapprox(get_value(sf, 0.5), 0, atol=1e-15)
    @test isapprox(get_value(sf, 0.25), 0, atol=1e-15)

    sϕ = SinusoidalInput(2, 1, π/2, 0, Inf)

    @test isapprox(get_value(sϕ, 0), 1.0)
    @test isapprox(get_value(sϕ, 1), 1.0, atol=1e-15)
end


text:
using QuantumLab
using Base.Test

info("$(now())  SETTING UP VARIABLES FOR FURTHER TESTING")
indent="=> "
include("setup.jl")

# test DocumentationModule
info("$(now())  TESTING:   DocumentationModule")
docModTester() = return 1
@doc """
some markdown documentation
""" docModTester
@test typeof(@doc(docModTester)) == Base.Markdown.MD
@add_doc GenericCitation("generic citation") docModTester
@test typeof(@doc(docModTester)) == Vector{DocumentationModule.Documentation}
@add_doc BookCitation(["C. Darwin"], "On the Origin of Species", "978-0451529060") docModTester
@test typeof(@doc(docModTester)) == Vector{DocumentationModule.Documentation}
@add_doc JournalCitation(["M. Mustermann"],"J. Stup. Mistakes",1,12,2016) docModTester
@test typeof(@doc(docModTester)) == Vector{DocumentationModule.Documentation}
@add_doc Citation(JournalCitation(["M. Mustermann"],"J. Stup. Mistakes",1,12,2016)) docModTester
@test typeof(@doc(docModTester)) == Vector{DocumentationModule.Documentation}
@test isa([Base.Markdown.parse("Hellau!"),BookCitation(["C. Darwin"], "On the Origin of Species", "978-0451529060")],Vector{DocumentationModule.Documentation})
@test isa([Base.Markdown.parse("Hellau!"),JournalCitation(["M. Mustermann"],"J. Stup. Mistakes",1,12,2016)],Vector{DocumentationModule.Documentation})
@test isa([Base.Markdown.parse("Hellau!"),DocumentationModule.Documentation(BookCitation(["C. Darwin"], "On the Origin of Species", "978-0451529060"))],Vector{DocumentationModule.Documentation})
@test isa([GenericCitation("Hellau!"),DocumentationModule.Documentation(BookCitation(["C. Darwin"], "On the Origin of Species", "978-0451529060"))],Vector{DocumentationModule.Documentation})

info("$(now())  TESTING:   Basic Functionality")
# test AtomModule
@test Element("C") == Element("C")
el1 = Element("C")
el2 = Element("C")
@test isequal(el1,el2)
@test !isequal(el1,Element("Na"))
# test BaseModule
mqn = MQuantumNumber(0,0,0)
for (mqn in MQuantumNumbers(LQuantumNumber("D"))) end
@test mqn == MQuantumNumber(0,0,2)
@test distance(Position(0,0,1),Position(1,0,0)) == sqrt(2)
@test length(MQuantumNumbers(LQuantumNumber("P"))) == 3
@test 2*Position(1,2,3) == Position(1,2,3)*2
@test isless(LQuantumNumber("S"),LQuantumNumber("P"))
@test_approx_eq 0. BaseModule.evaluateFunction(Position(1.,2.,3.), x->x.x+x.y-x.z)
@test h2o.atoms[2].element.symbol == "O"
@test_approx_eq_eps -1.4191843 readGeometryXYZ("h2o.xyz").atoms[3].position.x 1e-7
# test readBasisSetTX93
@test sto3g.definitions[Element("C")][1].lQuantumNumber.symbol == "S"
@test sto3g.definitions[Element("C")][1].primitives[1].exponent == 71.616837
# test BasisModule
@test_approx_eq_eps -0.7332137 bas.contractedBFs[3].primitiveBFs[1].center.y 1e-7
@test_approx_eq_eps 0.35175381 BasisFunctionsModule.evaluateFunction(origin, bas.contractedBFs[3]) 1e-8

# test IntegralsModule
info("$(now())  TESTING:   IntegralsModule")
@test_approx_eq IntegralsModule.GaussianIntegral1D_Valeev(1,.4) 0.
@test_approx_eq 0.0 matrixOverlap[6,1]
@test_approx_eq_eps 0.3129324434238492 matrixOverlap[4,1] 1e-8
@test_approx_eq_eps 29.00 matrixKinetic[2,2] 1e-2
@test_approx_eq_eps 0.16175843 IntegralsModule.FIntegral(2,0.3) 1e-8
@test_approx_eq -π IntegralsModule.computeIntegralNuclearAttraction(PrimitiveGaussianBasisFunction(Position(0,0,0),1.,MQuantumNumber(1,0,0)),PrimitiveGaussianBasisFunction(Position(0,0,0),1.,MQuantumNumber(1,0,0)),Atom(Element("C"),Position(0,0,0)))
@test_approx_eq IntegralsModule.GaussianIntegral1D_Mathematica(6,1.2)  IntegralsModule.GaussianIntegral1D_Valeev(6,1.2)
@test_approx_eq IntegralsModule.OverlapFundamental(bas.contractedBFs[1].primitiveBFs[1],bas.contractedBFs[1].primitiveBFs[2]) computeIntegralOverlap(bas.contractedBFs[1].primitiveBFs[1],bas.contractedBFs[1].primitiveBFs[2])

# test InitialGuessModule
info("$(now())  TESTING:   InitialGuessModule")
@test_approx_eq_eps -0.264689 mean(matrixSADguess)[3,2] 1e-6
@test_throws ErrorException computeDensityGuessSAD("NotImplemented","STO-3G",h2o)

# test HartreeFockModule
info("$(now())  TESTING:   HartreeFock (bas, h2o, density)")
@test_approx_eq_eps -74.96178985 (evaluateSCF(bas,h2o,density,5,info=false,detailedinfo=false)[2] + computeEnergyInteratomicRepulsion(h2o)) 1e-7 # checked against FermiONs++
#@test_approx_eq -4.473355520007 mean(HartreeFockModule.evaluateSCFStep(bas,h2o,mean(matrixSADguess),matrixOverlap,5)[1])
#@test_approx_eq_eps -74.96178985 (computeEnergyHartreeFock(bas,h2o,density) + computeEnergyInteratomicRepulsion(h2o)) 1e-7 # checked against FermiONs++

# test Shells: ShellModule and LibInt2Module
info("$(now())  TESTING:   Shells/LibInt2Shells/Basis equivalence")
shell_native  = Shell(Position(0.,0.,0.),LQuantumNumber("S"),[1.,2.,3.],[.1,.2,.3],renorm=false)
shell_libint2 = LibInt2Shell([0.,0.,0.],0,3,[1.,2.,3.],[.1,.2,.3])
shell_nativefromlibint2 = Shell(shell_libint2)
@test_approx_eq_eps shell_nativefromlibint2.coefficients[2] 0.41030724 1e-8
@test_approx_eq shell_native.coefficients[2] .2
@test_approx_eq computeMatrixBlockOverlap(shell_nativefromlibint2,shell_nativefromlibint2) computeMatrixBlockOverlap(shell_libint2,shell_libint2)
@test_approx_eq computeTensorBlockElectronRepulsionIntegrals(shell_nativefromlibint2,shell_nativefromlibint2,shell_nativefromlibint2,shell_nativefromlibint2) computeTensorBlockElectronRepulsionIntegrals(shell_libint2,shell_libint2,shell_libint2,shell_libint2)
@test_approx_eq LibInt2Module.nbf(shell_libint2) 1
import QuantumLab.libint2_available
if (libint2_available)
  @test_throws ErrorException LibInt2Shell([0.,0.,0.],100,1,[.1],[.5])
end

# test LibInt2Module
info("$(now())  TESTING:   LibInt2Module")
libInt2Finalize()
libInt2Initialize()
tmp_lib = LibInt2Shell([0.,1.,2.],1,2,[1.,2.],[0.5,0.5];renorm=false)
tmp = Shell(tmp_lib).coefficients
@test_approx_eq tmp[1] tmp[2]
destroy!(tmp_lib)
S = computeMatrixOverlap(shells)
@test_approx_eq S computeMatrixOverlap(shells_native)
@test_approx_eq S matrixOverlap
T = computeMatrixKinetic(shells)
@test_approx_eq T computeMatrixKinetic(shells_native)
@test_approx_eq T matrixKinetic
J = computeMatrixCoulomb(shells,density)
@test_approx_eq_eps J computeMatrixCoulomb(bas,density) 1e-8
@test_approx_eq computeMatrixBlockOverlap(shells[1],shells[1]) S[1,1]
@test_approx_eq computeMatrixBlockKinetic(shells[1],shells[1]) T[1,1]
@test_approx_eq computeTensorBlockElectronRepulsionIntegrals(shells[1],shells[5],shells[4],shells[4])[1,1,2,2] computeElectronRepulsionIntegral(bas.contractedBFs[1],bas.contractedBFs[7],bas.contractedBFs[5],bas.contractedBFs[5])

# test LaplaceModule
info("$(now())  TESTING:   LaplaceModule")
if(!isdir("hackbusch_pretables"))
  downloadLaplacePointsHackbusch("hackbusch_pretables")
end
R = transformRangeToIdealLaplace(0.5,3.)[2]
lp = transformLaplacePointFromIdealLaplace( findLaplacePointsHackbuschPretableLarger(15,R,"hackbusch_pretables")[1], 0.5)
@test_throws ErrorException findLaplacePointsHackbuschPretableSmaller(15,R,"hackbusch_pretables")
@test_approx_eq_eps LaplaceModule.computeInverseByLaplaceApproximation(2.3,lp) 1./2.3 1e-7
@test_approx_eq_eps LaplaceModule.computeRPADenominatorByDoubleLaplace(1.2,2.3,lp) 1./(1.2^2 + 2.3^2) 1e-3
println("Expecting warning below:")
@test findLaplacePointsHackbuschPretableLarger(2,100.,"hackbusch_pretables") == findLaplacePointsHackbuschPretableLarger(2,50.,"hackbusch_pretables")
@test transformLaplacePointFromIdealLaplace( findLaplacePointsHackbuschPretableSmaller(15,transformRangeToIdealLaplace(0.5,6.)[2],"hackbusch_pretables")[1], 0.5) == lp

# test RI-Module
info("$(now())  TESTING:   RIModule")
@test_approx_eq 0.3950752513027109 mean(computeMatrixExchangeRIK(bas,bas,matrixSADguess[1]))
tensorRICoulomb = computeTensorElectronRepulsionIntegralsRICoulomb(bas,bas)
@test_approx_eq 0.0429373705056905 mean(tensorRICoulomb)
@test_approx_eq 0.0012436634924295 tensorRICoulomb[1,3,4,5]
tensorRIOverlap = computeTensorElectronRepulsionIntegralsRIOverlap(bas,bas)
@test_approx_eq 0.0144597945326691 tensorRIOverlap[1,3,4,5]
@test_approx_eq 0.0537036598506078 mean(tensorRIOverlap)

# test SpecialMatricesModule
info("$(now())  SpecialMatricesModule")
@test_approx_eq_eps -0.785008186026 mean(computeMatrixFock(bas,h2o,matrixSADguess[1])) 1e-10                  #### !!!! ####

### display functions ###
# These are simply here, so that the test coverage isn't limited by the display functions.
# They do not however test for anything really, so do not trust them. If you do add a *real*
# test for a display function, please add it to the corresponding module block above instead
# of here.
info("$(now())  DISPLAY FUNCTIONS")
display(bas)
display(bas.contractedBFs[3])
display(JournalCitation(["M. Mustermann"],"J. Stup. Mistakes",1,12,2016))
display(GenericCitation("M. Mustermann personal note"))
display(BookCitation(["C. Darwin"], "On the Origin of Species", "978-0451529060"))
display([BookCitation(["C. Darwin"], "On the Origin of Species", "978-0451529060"), JournalCitation(["M. Mustermann"],"J. Stup. Mistakes",1,12,2016)])
display([Base.Markdown.parse("**Why, god, why??**"),BookCitation(["C. Darwin"], "On the Origin of Species", "978-0451529060"), JournalCitation(["M. Mustermann"],"J. Stup. Mistakes",1,12,2016)])
display([JournalCitation(["M. Mustermann"],"J. Stup. Mistakes",1,12,2016),JournalCitation(["M. Mustermann"],"J. Stup. Mistakes",1,12,2016)])
display(@doc(IntegralsModule.GaussianIntegral1D_Valeev))
display(shells[1])
display(shell_native)
display(lp)
summarize(BasisModule.Basis)
summarize(GaussianBasis)

info("\n\n     CONGRATULATIONS!! ALL TESTS COMPLETED!")


text:
macro cuda(lib,fun,x...)
    if Libdl.find_library(["lib$lib"], []) == ""
        msg = "Cannot find lib$lib, please install it and rerun Pkg.build(\"Knet\")."
        :(error($msg))
    else
        f2 = ("$fun","lib$lib")
        fx = Expr(:ccall, f2, :UInt32, x...)
        err = "$lib.$fun error "
        quote
            local _r = $fx
            if _r != 0
                warn($err, _r)
                Base.show_backtrace(STDOUT, backtrace())
            end
        end
    end
end

let GPU=-1, handles=Dict()
    global gpu, cublashandle, cudnnhandle

    gpu()=GPU  # (d=Cint[-1];@cuda(cudart,cudaGetDevice,(Ptr{Cint},),d);d[1])

    function gpu(i::Int)
        if 0 <= i < gpucount()
            @cuda(cudart,cudaSetDevice, (Cint,), i)
            cublashandle = get!(cublasCreate, handles, (:cublas,i))
            cudnnhandle  = get!(cudnnCreate, handles, (:cudnn,i))
        else
            i = -1
            cublashandle = cudnnhandle = nothing
        end
        return (GPU = i)
    end

    function gpu(b::Bool)
        if b
            pick = mem = -1
            for i=0:gpucount()-1
                @cuda(cudart,cudaSetDevice, (Cint,), i)
                imem = gpufree()
                if imem > mem
                    pick = i
                    mem = imem
                end
            end
            gpu(pick)
        else
            gpu(-1)
        end
    end
end

"""
gpu() returns the id of the active GPU device or -1 if none are active.

gpu(d::Int) activates the GPU device d if 0 <= d < gpucount().

gpu(true) activates the GPU device with the most available memory.

gpu(false) deactivates GPU devices.
""" gpu

gpucount()=(try; p=Cint[0]; @cuda(cudart,cudaGetDeviceCount,(Ptr{Cint},),p); p[1]; catch; 0; end)
gpumem()=(f=Csize_t[0];m=Csize_t[0]; @cuda(cudart,cudaMemGetInfo,(Ptr{Csize_t},Ptr{Csize_t}),f,m); (Int(f[1]),Int(m[1])))
gpufree()=gpumem()[1]
gpuinfo(msg="")=(print("$msg "); println((gpumem()...,meminfo()...)))
gpusync()=@cuda(cudart,cudaDeviceSynchronize,())

typealias Cptr Ptr{Void}

function cublasCreate()
    handleP = Cptr[0]
    @cuda(cublas,cublasCreate_v2, (Ptr{Cptr},), handleP)
    handle = handleP[1]
    atexit(()->@cuda(cublas,cublasDestroy_v2, (Cptr,), handle))
    return handle
end

function cudnnCreate()
    handleP = Cptr[0]
    @cuda(cudnn,cudnnCreate,(Ptr{Cptr},), handleP)
    handle = handleP[1]
    atexit(()->@cuda(cudnn,cudnnDestroy,(Cptr,), handle))
    return handle
end


text:
"""
    struct GlobalMatrices

Stores global matrices and some other things.
Everything that is very expensive to compute which stays
constant throughout the simulation should be stored here.

    GlobalMatrices(basis::Basis, filter::Filter, dimensions)

Initialize GlobalMatrices for `basis`, `filter` and `dimensions`.
"""
 struct GlobalMatrices
     normalsigns::Dict{Face, Int64}
     normalidxs::Dict{Face, Int64}
     oppositefaces::Dict{Face, Face}
     project_dofs_to_face::Dict{Face,Array{Float64,2}}
     project_flux_to_face::Dict{Face,Array{Float64,2}}
     project_dofs_from_face::Dict{Face,Array{Float64,2}}
     basissize_1d::Int64
     quadweights_nd::Array{Float64, 1}
     reference_massmatrix_face::Array{Float64,2}
     reference_massmatrix_cell::Array{Float64,2}
     reference_derivative_matrix::Array{Float64,2}

     filter_matrix::Array{Float64,2}

     function GlobalMatrices(basis::Basis, filter::Filter, dimensions)
        faces = instances(TerraDG.Face)

        normalsigns = Dict(left => -1, right => 1, top => 1, bottom => -1)
        normalidxs = Dict(left => 1, right => 1, top => 2, bottom => 2)
        oppositefaces = Dict(left => right, right => left, top => bottom, bottom => top)

        project_dofs_to_face = Dict(f => face_projection_matrix(basis, f) for f in faces)
        project_flux_to_face = Dict(f => face_projection_matrix(basis, f) for f in faces)
        project_dofs_from_face = Dict(f => transpose(face_projection_matrix(basis, f)) for f in faces)

        basissize_1d = size(basis,1)
        basissize_nd = length(basis)
        quadweights_nd = basis.quadweights;

        reference_massmatrix_cell = massmatrix(basis, dimensions)
        reference_massmatrix_face = massmatrix(basis, dimensions - 1)
        reference_derivative_matrix = derivativematrix(basis)

        filter_matrix = evaluate_filter_matrix(filter, basis)

         new(normalsigns, normalidxs, oppositefaces, project_dofs_to_face,
         project_flux_to_face, project_dofs_from_face, basissize_1d, quadweights_nd,
         reference_massmatrix_face,
         reference_massmatrix_cell, reference_derivative_matrix,
         filter_matrix)
     end
 end


text:
import Base: _throw_dmrs
import DiskArrays: splittuple
#Reshaping is really not trivial, because the access pattern would completely change for reshaped arrays,
#rectangles would not remain rectangles in the parent array. However, we can support the case where only
#singleton dimensions are added, later we could allow more special cases like joining two dimensions to one
struct ReshapedDiskArray{T,N,P<:AbstractArray{T},M} <: AbstractDiskArray{T,N}
    parent::P
    keepdim::NTuple{M,Int}
    newsize::NTuple{N,Int}
end
Base.size(r::ReshapedDiskArray) = r.newsize
haschunks(a::ReshapedDiskArray) = haschunks(a.parent)
function eachchunk(a::ReshapedDiskArray{<:Any,N}) where N
  pchunks = eachchunk(a.parent)
  inow::Int=0
  outchunks = ntuple(N) do idim
    if in(idim,a.keepdim)
      inow+=1
      pchunks.chunksize[inow],pchunks.offset[inow]
    else
      (1,0)
    end
  end
  GridChunks(a,first.(outchunks), offset=last.(outchunks))
end

function DiskArrays.readblock!(a::ReshapedDiskArray,aout,i::OrdinalRange...)
  inew = tuple_tuple_getindex(i,a.keepdim)
  DiskArrays.readblock!(a.parent,reshape(aout,map(length,inew)),inew...)
  nothing
end
tuple_tuple_getindex(t,i) = _ttgi((),t,i...)
_ttgi(o,t,i1,irest...) = _ttgi((o...,t[i1]),t,irest...)
_ttgi(o,t,i1) = (o...,t[i1])
function DiskArrays.writeblock!(a::ReshapedDiskArray,v,i::OrdinalRange...)
  inew = tuple_tuple_getindex(i,a.keepdim)
  DiskArrays.writeblock!(a.parent,reshape(v,map(length,inew)),inew...)
  nothing
end
function reshape_disk(parent, dims)
    n = length(parent)
    ndims(parent) > length(dims) && error("For DiskArrays, reshape is restricted to adding singleton dimensions")
    prod(dims) == n || _throw_dmrs(n, "size", dims)
    ipassed::Int=0
    s = size(parent)
    keepdim = map(s) do snow
        while true
            ipassed += 1
            d = dims[ipassed]
            if d>1
                d != snow && error("For DiskArrays, reshape is restricted to adding singleton dimensions")
                return ipassed
            end
        end
    end
    ReshapedDiskArray{eltype(parent),length(dims),typeof(parent),ndims(parent)}(parent,keepdim,dims)
end

import Base: _throw_dmrs
import Base.PermutedDimsArrays: genperm
struct PermutedDiskArray{T,N,P<:PermutedDimsArray{T,N}} <: AbstractDiskArray{T,N}
    a::P
end
function permutedims_disk(a,perm)
  pd = PermutedDimsArray(a,perm)
  PermutedDiskArray{eltype(a),ndims(a),typeof(pd)}(pd)
end
Base.size(r::PermutedDiskArray) = size(r.a)
haschunks(a::PermutedDiskArray) = haschunks(a.a.parent)
_getperm(a::PermutedDiskArray) = _getperm(a.a)
_getperm(::PermutedDimsArray{<:Any,<:Any,perm}) where perm = perm
_getiperm(a::PermutedDiskArray) = _getiperm(a.a)
_getiperm(::PermutedDimsArray{<:Any,<:Any,<:Any,iperm}) where iperm = iperm
function eachchunk(a::PermutedDiskArray)
  cc = eachchunk(a.a.parent)
  perm = _getperm(a.a)
  GridChunks(a,genperm(cc.chunksize,perm),offset=genperm(cc.offset,perm))
end
function DiskArrays.readblock!(a::PermutedDiskArray,aout,i::OrdinalRange...)
  iperm = _getiperm(a)
  inew = genperm(i, iperm)
  DiskArrays.readblock!(a.a.parent,PermutedDimsArray(aout,iperm),inew...)
  nothing
end
function DiskArrays.writeblock!(a::PermutedDiskArray,v,i::OrdinalRange...)
  iperm = _getiperm(a)
  inew = genperm(i, iperm)
  DiskArrays.writeblock!(a.a.parent,PermutedDimsArray(v,iperm),inew...)
  nothing
end

macro implement_reshape(t)
  quote
    Base._reshape(parent::$t, dims::NTuple{N,Int}) where N =
      reshape_disk(parent, dims)
  end
end

macro implement_permutedims(t)
  quote
    Base.permutedims(parent::$t, perm) =
      permutedims_disk(parent, perm)
  end
end


text:
@testset "Standalone Executable Integration" begin

    testpath = pwd()
    scratch = tempdir()
    cd(scratch)

    # --- Times table, file IO, mallocarray

    # Compile...
    # We have to start a new Julia process to get around the fact that Pkg.test
    # disables `@inbounds`, but ironically we can use `--compile=min` to make that
    # faster.
    status = run(`julia --compile=min $testpath/scripts/times_table.jl`)
    @test isa(status, Base.Process)
    @test status.exitcode == 0

    # Attempt to run
    println("5x5 times table:")
    status = run(`./times_table 5 5`)
    @test isa(status, Base.Process)
    @test status.exitcode == 0
    @test parsedlm(Int64, c"table.tsv", '\t') == (1:5)*(1:5)'

    # --- Random number generation

    # Compile...
    status = run(`julia --compile=min $testpath/scripts/rand_matrix.jl`)
    @test isa(status, Base.Process)
    @test status.exitcode == 0

    # Run...
    println("5x5 random matrix:")
    status = run(`./rand_matrix 5 5`)
    @test isa(status, Base.Process)
    @test status.exitcode == 0

    # --- Test LoopVectorization integration

@static if LoopVectorization.VectorizationBase.has_feature(Val{:x86_64_avx2})
    # Compile...
    status = run(`julia --compile=min $testpath/scripts/loopvec_product.jl`)
    @test isa(status, Base.Process)
    @test status.exitcode == 0

    # Run...
    println("10x10 table sum:")
    status = run(`./loopvec_product 10 10`)
    @test isa(status, Base.Process)
    @test status.exitcode == 0
    @test parsedlm(c"product.tsv",'\t')[] == 3025
end

    # Compile...
    status = run(`julia --compile=min $testpath/scripts/loopvec_matrix.jl`)
    @test isa(status, Base.Process)
    @test status.exitcode == 0

    # Run...
    println("10x3 matrix product:")
    status = run(`./loopvec_matrix 10 3`)
    @test isa(status, Base.Process)
    @test status.exitcode == 0
    A = (1:10) * (1:3)'
    @test parsedlm(c"table.tsv",'\t') == A' * A

    # --- Test string handling

    # Compile...
    status = run(`julia --compile=min $testpath/scripts/print_args.jl`)
    @test isa(status, Base.Process)
    @test status.exitcode == 0

    # Run...
    println("String indexing and handling:")
    status = run(`./print_args foo bar`)
    @test isa(status, Base.Process)
    @test status.exitcode == 0

    # --- Clean up
    cd(testpath)
end


text:
##### Moisture tendencies

#####
##### First order fluxes
#####

function flux(::Advect{Tracers{N}}, atmos, args) where {N}
    @unpack state = args
    u = state.ρu / state.ρ
    return (state.tracers.ρχ .* u')'
end

#####
##### Second order fluxes
#####

function flux(::Diffusion{Tracers{N}}, atmos, args) where {N}
    @unpack state, aux, t, diffusive = args
    ν, D_t, τ = turbulence_tensors(atmos, state, diffusive, aux, t)
    d_χ = (-D_t) * aux.tracers.δ_χ' .* diffusive.tracers.∇χ
    return d_χ * state.ρ
end

max_stars_repo_name: jm-c/CLIMA

max_issues_repo_name: jm-c/CLIMA

max_forks_repo_name: jm-c/CLIMA