text:
# ============================================================
# Structs of Control Volume
# Solver stores data in arrays of struct (AoS)
# ============================================================

# ------------------------------------------------------------
# 1D
# ------------------------------------------------------------

"""
    mutable struct ControlVolume1D{A,B} <: AbstractControlVolume1D
        w::A
        prim::B
        sw::A
    end

1D control volume with no distribution function

"""
mutable struct ControlVolume1D{A,B} <: AbstractControlVolume1D
    w::A
    prim::B
    sw::A
end

function ControlVolume1D(W, PRIM)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = zero(W)

    return ControlVolume1D{typeof(w),typeof(prim)}(w, prim, sw)
end

function Base.show(io::IO, ctr::ControlVolume1D{A,B}) where {A,B}
    print(
        io,
        "ControlVolume1D{$A,$B}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
    )
end

"""
    mutable struct ControlVolume1D1F{A,B} <: AbstractControlVolume1D
        w::A
        prim::A
        sw::A

        f::B
        sf::B
    end

1D control volume with 1 distribution function

"""
mutable struct ControlVolume1D1F{A,B} <: AbstractControlVolume1D
    w::A
    prim::A
    sw::A

    f::B
    sf::B
end

function ControlVolume1D1F(W::T, PRIM::T, F) where T
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = zero(W)

    f = deepcopy(F)
    sf = zero(f)

    return ControlVolume1D1F{typeof(w),typeof(f)}(w, prim, sw, f, sf)
end

function Base.show(io::IO, ctr::ControlVolume1D1F{A,B}) where {A,B}
    print(
        io,
        "ControlVolume1D1F{$A,$B}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
        "pdf vars: f\n",
        "pdf slopes: sf\n",
    )
end

"""
    mutable struct ControlVolume1D2F{A,B} <: AbstractControlVolume1D
        w::A
        prim::A
        sw::A

        h::B
        b::B
        sh::B
        sb::B
    end

1D control volume with 2 distribution functions

"""
mutable struct ControlVolume1D2F{A,B} <: AbstractControlVolume1D
    w::A
    prim::A
    sw::A

    h::B
    b::B
    sh::B
    sb::B
end

function ControlVolume1D2F(
    W::T1,
    PRIM::T1,
    H::T2,
    B::T2,
) where {T1<:AbstractArray,T2<:AbstractArray}

    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = zero(W)

    h = deepcopy(H)
    b = deepcopy(B)
    sh = zero(h)
    sb = zero(b)

    return ControlVolume1D2F{typeof(w),typeof(h)}(w, prim, sw, h, b, sh, sb)

end

function Base.show(io::IO, ctr::ControlVolume1D2F{A,B}) where {A,B}
    print(
        io,
        "ControlVolume1D2F{$A,$B}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
        "pdf vars: h, b\n",
        "pdf slopes: sh, sb\n",
    )
end

"""
    mutable struct ControlVolume1D3F{A,B,C,D,E} <: AbstractControlVolume1D
        w::A
        prim::A
        sw::A

        h0::B
        h1::B
        h2::B
        sh0::B
        sh1::B
        sh2::B

        E::C
        B::C
        ϕ::D
        ψ::D
        lorenz::E
    end

1D control volume with 3 distribution functions

"""
mutable struct ControlVolume1D3F{A,B,C,D,E} <: AbstractControlVolume1D
    w::A
    prim::A
    sw::A

    h0::B
    h1::B
    h2::B
    sh0::B
    sh1::B
    sh2::B

    E::C
    B::C
    ϕ::D
    ψ::D
    lorenz::E
end

# deterministic
function ControlVolume1D3F(
    W::AbstractArray{<:Real,2},
    PRIM::AbstractArray{<:Real,2},
    H0::AbstractArray{<:AbstractFloat,3},
    H1::AbstractArray{<:AbstractFloat,3},
    H2::AbstractArray{<:AbstractFloat,3},
    E0::AbstractArray{<:AbstractFloat,1},
    B0::AbstractArray{<:AbstractFloat,1},
    L::AbstractArray{<:AbstractFloat,2},
)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = zero(W)

    h0 = deepcopy(H0)
    h1 = deepcopy(H1)
    h2 = deepcopy(H2)
    sh0 = zero(H0)
    sh1 = zero(H1)
    sh2 = zero(H2)

    E = deepcopy(E0)
    B = deepcopy(B0)
    ϕ = 0.0
    ψ = 0.0
    lorenz = deepcopy(L)

    return ControlVolume1D3F{typeof(w),typeof(h0),typeof(E),typeof(ϕ),typeof(lorenz)}(
        w,
        prim,
        sw,
        h0,
        h1,
        h2,
        sh0,
        sh1,
        sh2,
        E,
        B,
        ϕ,
        ψ,
        lorenz,
    )
end

# Rykov
function ControlVolume1D3F(
    W::AbstractArray{<:Real,1},
    PRIM::AbstractArray{<:Real,1},
    H0::AbstractArray{<:AbstractFloat,1},
    H1::AbstractArray{<:AbstractFloat,1},
    H2::AbstractArray{<:AbstractFloat,1},
)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = zero(W)

    h0 = deepcopy(H0)
    h1 = deepcopy(H1)
    h2 = deepcopy(H2)
    sh0 = zero(H0)
    sh1 = zero(H1)
    sh2 = zero(H2)

    E = nothing
    B = nothing
    ϕ = nothing
    ψ = nothing
    lorenz = nothing

    return ControlVolume1D3F{typeof(w),typeof(h0),typeof(E),typeof(ϕ),typeof(lorenz)}(
        w,
        prim,
        sw,
        h0,
        h1,
        h2,
        sh0,
        sh1,
        sh2,
        E,
        B,
        ϕ,
        ψ,
        lorenz,
    )
end

# stochastic
function ControlVolume1D3F(
    W::AbstractArray{<:Real,3},
    PRIM::AbstractArray{<:Real,3},
    H0::AbstractArray{<:AbstractFloat,4},
    H1::AbstractArray{<:AbstractFloat,4},
    H2::AbstractArray{<:AbstractFloat,4},
    E0::AbstractArray{<:AbstractFloat,2},
    B0::AbstractArray{<:AbstractFloat,2},
    L::AbstractArray{<:AbstractFloat,3},
)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = zero(W)

    h0 = deepcopy(H0)
    h1 = deepcopy(H1)
    h2 = deepcopy(H2)
    sh0 = zero(H0)
    sh1 = zero(H1)
    sh2 = zero(H2)

    E = deepcopy(E0)
    B = deepcopy(B0)
    ϕ = zero(E[1, :]) # here is difference
    ψ = zero(B[1, :])
    lorenz = deepcopy(L)

    return ControlVolume1D3F{typeof(w),typeof(h0),typeof(E),typeof(ϕ),typeof(lorenz)}(
        w,
        prim,
        sw,
        h0,
        h1,
        h2,
        sh0,
        sh1,
        sh2,
        E,
        B,
        ϕ,
        ψ,
        lorenz,
    )
end

function Base.show(io::IO, ctr::ControlVolume1D3F{A,B,C,D,E}) where {A,B,C,D,E}
    print(
        io,
        "ControlVolume1D3F{$A,$B,$C,$D,$E}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
        "pdf vars: h0, h1, h2\n",
        "pdf slopes: sh0, sh1, sh2\n",
        "electric field: $(ctr.E)\n",
        "magnetic field: $(ctr.B)\n",
        "Lorenz force: $(ctr.lorenz)\n",
    )
end

"""
    mutable struct ControlVolume1D4F{A,B,C,D,E} <: AbstractControlVolume1D
        w::A
        prim::A
        sw::A

        h0::B
        h1::B
        h2::B
        h3::B
        sh0::B
        sh1::B
        sh2::B
        sh3::B

        E::C
        B::C
        ϕ::D
        ψ::D
        lorenz::E
    end

1D control volume with 4 distribution functions

"""
mutable struct ControlVolume1D4F{A,B,C,D,E} <: AbstractControlVolume1D
    w::A
    prim::A
    sw::A

    h0::B
    h1::B
    h2::B
    h3::B
    sh0::B
    sh1::B
    sh2::B
    sh3::B

    E::C
    B::C
    ϕ::D
    ψ::D
    lorenz::E
end

#--- deterministic ---#
function ControlVolume1D4F(
    W::AbstractArray{<:Real,2},
    PRIM::AbstractArray{<:Real,2},
    H0::AbstractArray{<:AbstractFloat,2},
    H1::AbstractArray{Float64,2},
    H2::AbstractArray{Float64,2},
    H3::AbstractArray{Float64,2},
    E0::AbstractArray{Float64,1},
    B0::AbstractArray{Float64,1},
    L::AbstractArray{Float64,2},
)

    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = zero(W)

    h0 = deepcopy(H0)
    h1 = deepcopy(H1)
    h2 = deepcopy(H2)
    h3 = deepcopy(H3)
    sh0 = deepcopy(H0)
    sh1 = deepcopy(H1)
    sh2 = deepcopy(H2)
    sh3 = deepcopy(H3)

    E = deepcopy(E0)
    B = deepcopy(B0)
    ϕ = 0.0
    ψ = 0.0
    lorenz = deepcopy(L)

    return ControlVolume1D4F{typeof(w),typeof(h0),typeof(E),typeof(ϕ),typeof(lorenz)}(
        w,
        prim,
        sw,
        h0,
        h1,
        h2,
        h3,
        sh0,
        sh1,
        sh2,
        sh3,
        E,
        B,
        ϕ,
        ψ,
        lorenz,
    )

end

#--- uncertainty quantification ---#
function ControlVolume1D4F(
    W::AbstractArray{<:Real,3},
    PRIM::AbstractArray{<:Real,3},
    H0::AbstractArray{<:AbstractFloat,3},
    H1::AbstractArray{Float64,3},
    H2::AbstractArray{Float64,3},
    H3::AbstractArray{Float64,3},
    E0::AbstractArray{Float64,2},
    B0::AbstractArray{Float64,2},
    L::AbstractArray{Float64,3},
)

    w = deepcopy(W)
    prim = deepcopy(PRIM)
    sw = zero(W)

    h0 = deepcopy(H0)
    h1 = deepcopy(H1)
    h2 = deepcopy(H2)
    h3 = deepcopy(H3)
    sh0 = zero(H0)
    sh1 = zero(H1)
    sh2 = zero(H2)
    sh3 = zero(H3)

    E = deepcopy(E0)
    B = deepcopy(B0)
    ϕ = zero(E[1, :])
    ψ = zero(B[1, :])
    lorenz = deepcopy(L)

    return ControlVolume1D4F{typeof(w),typeof(h0),typeof(E),typeof(ϕ),typeof(lorenz)}(
        w,
        prim,
        sw,
        h0,
        h1,
        h2,
        h3,
        sh0,
        sh1,
        sh2,
        sh3,
        E,
        B,
        ϕ,
        ψ,
        lorenz,
    )

end

function Base.show(io::IO, ctr::ControlVolume1D4F{A,B,C,D,E}) where {A,B,C,D,E}
    print(
        io,
        "ControlVolume1D4F{$A,$B,$C,$D,$E}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
        "pdf vars: h0, h1, h2, h3\n",
        "pdf slopes: sh0, sh1, sh2, sh3\n",
        "electric field: $(ctr.E)\n",
        "magnetic field: $(ctr.B)\n",
        "Lorenz force: $(ctr.lorenz)\n",
    )
end

# ------------------------------------------------------------
# 2D
# ------------------------------------------------------------

"""
    mutable struct ControlVolume2D{A,B} <: AbstractControlVolume2D
        w::A
        prim::A
        sw::B
    end

2D control volume with no distribution function

"""
mutable struct ControlVolume2D{A,B} <: AbstractControlVolume2D
    w::A
    prim::A
    sw::B
end

function ControlVolume2D(W, PRIM)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    #sw = zeros(eltype(W), (axes(W)..., Base.OneTo(2)))
    sw = slope_array(W)

    return ControlVolume2D{typeof(w),typeof(sw)}(w, prim, sw)
end

function Base.show(io::IO, ctr::ControlVolume2D{A,B}) where {A,B}
    print(
        io,
        "ControlVolume2D{$A,$B}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
    )
end

"""
    mutable struct ControlVolume2D1F{A,B,C,D} <: AbstractControlVolume2D
        w::A
        prim::A
        sw::B

        f::C
        sf::D
    end

2D control volume with 1 distribution function

"""
mutable struct ControlVolume2D1F{A,B,C,D} <: AbstractControlVolume2D
    w::A
    prim::A
    sw::B

    f::C
    sf::D
end

function ControlVolume2D1F(W, PRIM, F::AbstractArray)
    w = deepcopy(W)
    prim = deepcopy(PRIM)
    #sw = zeros(eltype(W), (axes(W)..., Base.OneTo(2)))
    sw = slope_array(W)

    f = deepcopy(F)
    #sf = zeros(eltype(F), (axes(F)..., Base.OneTo(2)))
    sf = slope_array(F)

    return ControlVolume2D1F{typeof(w),typeof(sw),typeof(f),typeof(sf)}(w, prim, sw, f, sf)
end

function Base.show(io::IO, ctr::ControlVolume2D1F{A,B,C,D}) where {A,B,C,D}
    print(
        io,
        "ControlVolume2D1F{$A,$B,$C,$D}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
        "pdf vars: f\n",
        "pdf slopes: sf\n",
    )
end

"""
    mutable struct ControlVolume2D2F{A,B,C,D} <: AbstractControlVolume2D
        w::A
        prim::A
        sw::B

        h::C
        b::C
        sh::D
        sb::D
    end

2D control volume with 2 distribution functions

"""
mutable struct ControlVolume2D2F{A,B,C,D} <: AbstractControlVolume2D
    w::A
    prim::A
    sw::B

    h::C
    b::C
    sh::D
    sb::D
end

function ControlVolume2D2F(
    W::AbstractArray,
    PRIM::AbstractArray,
    H::AbstractArray,
    B::AbstractArray,
)

    w = deepcopy(W)
    prim = deepcopy(PRIM)
    #sw = zeros(eltype(W), (axes(W)..., Base.OneTo(2)))
    sw = slope_array(W)

    h = deepcopy(H)
    b = deepcopy(B)
    #sh = zeros(eltype(H), (axes(H)..., Base.OneTo(2)))
    #sb = zeros(eltype(B), (axes(B)..., Base.OneTo(2)))
    sh = slope_array(H)
    sb = slope_array(B)

    return ControlVolume2D2F{typeof(w),typeof(sw),typeof(h),typeof(sh)}(
        w,
        prim,
        sw,
        h,
        b,
        sh,
        sb,
    )

end

function Base.show(io::IO, ctr::ControlVolume2D2F{A,B,C,D}) where {A,B,C,D}
    print(
        io,
        "ControlVolume2D2F{$A,$B,$C,$D}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
        "pdf vars: h, b\n",
        "pdf slopes: sh, sb\n",
    )
end

"""
    mutable struct ControlVolume2D3F{T2,T3,T4,T5,T6,T7,T8} <: AbstractControlVolume2D
        w::T2
        prim::T2
        sw::T3

        h0::T4
        h1::T4
        h2::T4
        sh0::T5
        sh1::T5
        sh2::T5

        E::T6
        B::T6
        ϕ::T7
        ψ::T7
        lorenz::T8
    end

2D control volume with 3 distribution functions

"""
mutable struct ControlVolume2D3F{T2,T3,T4,T5,T6,T7,T8} <: AbstractControlVolume2D
    w::T2
    prim::T2
    sw::T3

    h0::T4
    h1::T4
    h2::T4
    sh0::T5
    sh1::T5
    sh2::T5

    E::T6
    B::T6
    ϕ::T7
    ψ::T7
    lorenz::T8
end

#--- deterministic & stochastic ---#
function ControlVolume2D3F(
    W::AbstractArray,
    PRIM::AbstractArray,
    H0::AbstractArray,
    H1::AbstractArray,
    H2::AbstractArray,
    E0::AbstractArray,
    B0::AbstractArray,
    L::AbstractArray,
)

    w = deepcopy(W)
    prim = deepcopy(PRIM)
    #sw = zeros(eltype(W), (axes(W)..., Base.OneTo(2))) # 2D
    sw = slope_array(W)

    h0 = deepcopy(H0)
    h1 = deepcopy(H1)
    h2 = deepcopy(H2)
    #sh0 = zeros(eltype(H0), (axes(H0)..., Base.OneTo(2)))
    #sh1 = zeros(eltype(H1), (axes(H1)..., Base.OneTo(2)))
    #sh2 = zeros(eltype(H2), (axes(H2)..., Base.OneTo(2)))
    sh0 = slope_array(H0)
    sh1 = slope_array(H1)
    sh2 = slope_array(H2)

    E = deepcopy(E0)
    B = deepcopy(B0)
    ϕ = 0.0
    ψ = 0.0
    lorenz = deepcopy(L)

    return ControlVolume2D3F{
        typeof(w),
        typeof(sw),
        typeof(h0),
        typeof(sh2),
        typeof(E),
        typeof(ϕ),
        typeof(lorenz),
    }(
        w,
        prim,
        sw,
        h0,
        h1,
        h2,
        sh0,
        sh1,
        sh2,
        E,
        B,
        ϕ,
        ψ,
        lorenz,
    )

end

#--- Rykov ---#
function ControlVolume2D3F(
    W::AbstractArray{<:Real,1},
    PRIM::AbstractArray{<:Real,1},
    H0::AbstractArray{<:AbstractFloat,2},
    H1::AbstractArray{<:AbstractFloat,2},
    H2::AbstractArray{<:AbstractFloat,2},
)

    w = deepcopy(W)
    prim = deepcopy(PRIM)
    #sw = zeros(eltype(W), (axes(W)..., Base.OneTo(2))) # 2D
    sw = slope_array(W)

    h0 = deepcopy(H0)
    h1 = deepcopy(H1)
    h2 = deepcopy(H2)
    #sh0 = zeros(eltype(H0), (axes(H0)..., Base.OneTo(2)))
    #sh1 = zeros(eltype(H1), (axes(H1)..., Base.OneTo(2)))
    #sh2 = zeros(eltype(H2), (axes(H2)..., Base.OneTo(2)))
    sh0 = slope_array(H0)
    sh1 = slope_array(H1)
    sh2 = slope_array(H2)

    E = nothing
    B = nothing
    ϕ = nothing
    ψ = nothing
    lorenz = nothing

    return ControlVolume2D3F{
        typeof(w),
        typeof(sw),
        typeof(h0),
        typeof(sh2),
        typeof(E),
        typeof(ϕ),
        typeof(lorenz),
    }(
        w,
        prim,
        sw,
        h0,
        h1,
        h2,
        sh0,
        sh1,
        sh2,
        E,
        B,
        ϕ,
        ψ,
        lorenz,
    )

end

function Base.show(
    io::IO,
    ctr::ControlVolume2D3F{T1,T2,T3,T4,T5,T6,T7},
) where {T1,T2,T3,T4,T5,T6,T7}
    print(
        io,
        "ControlVolume2D3F{$T1,$T2,$T3,$T4,$T5,$T6,$T7}\n",
        "conservative vars: $(ctr.w)\n",
        "primitive vars: $(ctr.prim)\n",
        "conservative slopes: $(ctr.sw)\n",
        "pdf vars: h0, h1, h2\n",
        "pdf slopes: sh0, sh1, sh2\n",
        "electric field: $(ctr.E)\n",
        "magnetic field: $(ctr.B)\n",
        "Lorenz force: $(ctr.lorenz)\n",
    )
end

# ------------------------------------------------------------
# Unstructured cell
# ------------------------------------------------------------

"""
    mutable struct ControlVolumeUS{E,F,A,B,C} <: AbstractUnstructControlVolume
        n::E
        x::F
        dx::F

        w::A
        prim::B
        sw::C
    end

Unstructured control volume with no distribution function

"""
mutable struct ControlVolumeUS{E,F,A,B,C} <: AbstractUnstructControlVolume
    n::E
    x::F
    dx::F

    w::A
    prim::B
    sw::C
end

function ControlVolumeUS(N, X::T, DX::T, W, PRIM) where {T<:Union{Real,AbstractVector}}
    n = deepcopy(N)
    x = deepcopy(X)
    dx = deepcopy(DX)

    w = deepcopy(W)
    prim = deepcopy(PRIM)
    #sw = zeros(eltype(W), axes(W)..., length(N[1]))
    sw = ifelse(
        length(N[1]) == 2,
        slope_array(W; reduction = true),
        slope_array(W; reduction = false),
    )

    return ControlVolumeUS{typeof(n),typeof(x),typeof(w),typeof(prim),typeof(sw)}(
        n,
        x,
        dx,
        w,
        prim,
        sw,
    )
end

"""
    mutable struct ControlVolumeUS1F{E,F,A,B,C,D} <: AbstractUnstructControlVolume
        n::E
        x::F
        dx::F

        w::A
        prim::A
        sw::B

        f::C
        sf::D
    end

Unstructured control volume with 1 distribution function

"""
mutable struct ControlVolumeUS1F{E,F,A,B,C,D} <: AbstractUnstructControlVolume
    n::E
    x::F
    dx::F

    w::A
    prim::A
    sw::B

    f::C
    sf::D
end

function ControlVolumeUS1F(N, X, DX, W, PRIM, F::AbstractArray{T}) where T
    n = deepcopy(N)
    x = deepcopy(X)
    dx = deepcopy(DX)

    w = deepcopy(W)
    prim = deepcopy(PRIM)
    #sw = zeros(eltype(W), axes(W)..., length(N[1]))
    sw = ifelse(
        length(N[1]) == 2,
        slope_array(W; reduction = true),
        slope_array(W; reduction = false),
    )

    f = deepcopy(F)
    #sf = zeros(eltype(F), axes(F)..., length(N[1]))
    sf = ifelse(
        length(N[1]) == 2,
        slope_array(F; reduction = true),
        slope_array(F; reduction = false),
    )

    return ControlVolumeUS1F{typeof(n),typeof(x),typeof(w),typeof(sw),typeof(f),typeof(sf)}(
        n,
        x,
        dx,
        w,
        prim,
        sw,
        f,
        sf,
    )
end

"""
    mutable struct ControlVolumeUS2F{E,F,A,B,C,D} <: AbstractUnstructControlVolume
        n::E
        x::F
        dx::F

        w::A
        prim::A
        sw::B

        h::C
        b::C
        sh::D
        sb::D
    end

Unstructured control volume with 2 distribution functions

"""
mutable struct ControlVolumeUS2F{E,F,A,B,C,D} <: AbstractUnstructControlVolume
    n::E
    x::F
    dx::F

    w::A
    prim::A
    sw::B

    h::C
    b::C
    sh::D
    sb::D
end

function ControlVolumeUS2F(
    N,
    X,
    DX,
    W::T1,
    PRIM::T1,
    H::T2,
    B::T2,
) where {T1<:AbstractArray,T2<:AbstractArray}

    n = deepcopy(N)
    x = deepcopy(X)
    dx = deepcopy(DX)

    w = deepcopy(W)
    prim = deepcopy(PRIM)
    #sw = zeros(eltype(W), axes(W)..., length(N[1]))
    sw = ifelse(
        length(N[1]) == 2,
        slope_array(W; reduction = true),
        slope_array(W; reduction = false),
    )

    h = deepcopy(H)
    b = deepcopy(B)
    #sh = zeros(eltype(H), axes(H)..., length(N[1]))
    #sb = zeros(eltype(B), axes(B)..., length(N[1]))
    sh = ifelse(
        length(N[1]) == 2,
        slope_array(H; reduction = true),
        slope_array(H; reduction = false),
    )
    sb = ifelse(
        length(N[1]) == 2,
        slope_array(B; reduction = true),
        slope_array(B; reduction = false),
    )

    return ControlVolumeUS2F{typeof(n),typeof(x),typeof(w),typeof(sw),typeof(h),typeof(sh)}(
        n,
        x,
        dx,
        w,
        prim,
        sw,
        h,
        b,
        sh,
        sb,
    )

end


text:
try
    rm("test1.h5")
catch
end
try
    rm("test2.h5")
catch
end
function write_in_layers_repack(filename = "test1.h5", dset_name = "test")
    file_list = Glob.glob("*.unw")
    stack_size = (720, 720, length(file_list))

    h5open(filename, "cw") do f
        d_create(
            f,
            dset_name,
            datatype(Float32),
            dataspace(stack_size),
            # "chunk", (1, 1, size(stack_in, 3)),  # Seems better to repack?
        )
    end

    println("Reading layers and writing")
    arr_buf = zeros(Float32, (720, 720))
    h5open(filename, "cw") do hf
        dset = hf[dset_name]
        for (idx, f) in enumerate(file_list)
            arr_buf .= Sario.load(f, do_permute = true)
            dset[:, :, idx] = arr_buf
        end
    end
    InsarTimeseries.repack(filename, dset_name)

    return
end

function write_in_pixel_chunks(filename = "test2.h5", dset_name = "test")
    file_list = Glob.glob("*.unw")
    stack_size = (720, 720, length(file_list))
    m = Array{Float32,3}(undef, stack_size)
    println("Reading all layers into memory")
    for (idx, f) in enumerate(file_list)
        m[:, :, idx] .= Sario.load(f)
    end

    h5open(filename, "cw") do f
        d_create(
            f,
            dset_name,
            datatype(Float32),
            dataspace(stack_size),
            "chunk",
            (1, 1, size(m, 3)),  # Seems better to repack?
        )
    end

    println("writing pixelwise")
    h5open(filename, "r+") do f
        # NOTE: this seems to take 10x slower than repack
        # for col in 1:size(m, 1)
        #     for row in 1:size(m, 2)
        #         f[dset_name][row, col, :] = m[row, col, :]
        #     end
        # end
        f[dset_name][:, :, :] = m

    end
    return
end

@time write_in_layers_repack();
@time write_in_pixel_chunks();

ext: jl


text:
################################################################################
#
#  NfOrd/FactorBaseBound.jl: Bounds for the factor base for Buchmann
#                               algorithm
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2016 Tommy Hofmann
#
################################################################################

################################################################################
#
#  Apres Belabas, Diaz y Diaz, Friedmann
#
################################################################################

# The algorithm is described in
# Belabas, Diaz y Diaz, Friedmann: "Small generators for the ideal class group"

@doc Markdown.doc"""
    factor_base_bound_bdf(O::NfOrd) -> Int

Use the algorithm of Belabas, Diaz y Diaz and Friedmann to find $B$ such that
under GRH the ideal class group of $\mathcal O$ is generated by the prime
ideals of norm bounded by $B$.
"""
function factor_base_bound_bdf(O::NfOrd)
  return _factor_base_bound_bdf(O, 100.0, 100.0)
end

function _factorbase_bound_bdf_right_side(O::NfOrd, x0::Float64, D::Dict{Int, Array{Tuple{Int, Int}, 1}})
  K = nf(O)
  d = degree(K)
  r, s = signature(O)

  curval = Int(ceil(x0))

  R = ArbField(64, false)

  summ = R(0)

  #curval = 100
  logcurval = log(R(curval))

  # small helper function (is this fast?)
  function comp_summand(p::fmpz, m::Int)
    logp = log(R(p))

    pm2 = R(p)^(R(FlintZZ(m)//FlintZZ(2)))

    secondterm = R(1) - m*logp//logcurval

    return logp//pm2 * secondterm
  end

  function comp_summand(p::Int, m::Int)
    return comp_summand(fmpz(p), m)
  end

  p = 2

  while p < curval

    max_exp = _max_power_in(p, curval)

    #println("maximal power is $max_exp")

    lP = _prime_dec_cache(O, p, D)

    for P in lP
      Pnorm = fmpz(p)^P[1]
      if Pnorm < curval
        max_exp = _max_power_in(Pnorm, curval)

        for m in 1:max_exp
          summand = comp_summand(Pnorm, m)
          summ = summ + summand
        end
      end
    end
    p = next_prime(p)
  end

  y = 2*summ - (R(d)*(const_pi(R)^2//R(2)) + r*4*const_catalan(R))//logcurval
  return y::arb
end

function _factor_base_bound_bdf(O::NfOrd, x0::Float64 = 50.0, ste::Float64 = 20.0)
  K = nf(O)
  d = degree(K)
  r, s = signature(O)

  R = ArbField(64, false)

  summ = R(0)

  dec_cache = Dict{Int, Array{Tuple{Int, Int}, 1}}()

  D = log(R(abs(discriminant(O)))) -
      R(degree(K))*(const_euler(R) +
      log(8*const_pi(R))) -
      r * const_pi(R)//R(2)

  y = _factorbase_bound_bdf_right_side(O, x0, dec_cache)

  x1 = 2*x0

  while y < D
    x0 = x1
    x1 = 2*x0
    y = _factorbase_bound_bdf_right_side(O, x1, dec_cache)
  end

  dista = abs(x0-x1)

  while !( y > D && dista < ste)
    if y < D
      x1 = x0 + 3*dista/2
    else
      x1 = x0 - dista/2
    end

    dista = abs(x1-x0)

    x0 = x1
    y = _factorbase_bound_bdf_right_side(O, x0, dec_cache)
  end
  @assert ceil(x0) < typemax(Int)
  return Int(ceil(x0))
end

function _prime_dec_cache(O::NfOrd, p::Int, D::Dict{Int, Array{Tuple{Int, Int}, 1}})
  if haskey(D, p)
    return D[p]
  else
    D[p] = prime_decomposition_type(O, p)
    return D[p]
  end
end

################################################################################
#
#  Apres Bach
#
################################################################################

# The theorem is described in
# Bach: "Explicit bounds for primality testing and related problems"

@doc Markdown.doc"""
    factor_base_bound_bach(O::NfOrd) -> Int

Use the theorem of Bach to find $B$ such that
under GRH the ideal class group of $\mathcal O$ is generated by the prime
ideals of norm bounded by $B$.
"""
function factor_base_bound_bach(O::NfOrd)
  p = 64
  R = ArbField(p, false)
  if degree(O)==2
    r = ceil(6*log(R(abs(discriminant(O))))^2)
  else
    r = ceil(12*log(R(abs(discriminant(O))))^2)
  end
  (b, x) = unique_integer(r)
  @assert b
  @assert x < typemax(Int)
  return Int(x)::Int
end

################################################################################
#
#  Toplevel function
#
################################################################################

@doc Markdown.doc"""
    factor_base_bound_grh(O::NfOrd) -> Int

Returns an integer $B$, such that
under GRH the ideal class group of $\mathcal O$ is generated by the prime
ideals of norm bounded by $B$.
"""
function factor_base_bound_grh(O::NfOrd)
  return min(factor_base_bound_bdf(O), factor_base_bound_bach(O))
end

@doc Markdown.doc"""
    factor_base_bound_minkowski(O::NfOrd) -> Int

Returns an integer $B$, such that
the ideal class group of $\mathcal O$ is generated by the prime
ideals of norm bounded by $B$.
"""
function factor_base_bound_minkowski(O::NfOrd)
  r1, r2 = signature(O)
  n = degree(O)
  b = BigInt(round(floor((4/pi)^r2*factorial(BigInt(n))/BigInt(n)^n*sqrt(1.0*abs(discriminant(O))))))
  return fmpz(b)
end


text:
using SpecialFunctions
using Test
x = collect(0:0.05:50)
"""
Computes the loss of discharge in the Pitztal based on an exponential function.

$(SIGNATURES)
returns how much water is lossed to other catchment in the Pitztal
"""
function loss(Discharge, loss_parameter)
      loss = loss_parameter .* Discharge.^2
      loss[loss.>12.1] .= 12.1
      return loss
end
plot()

"""
Computes how high the discharge would be without loss.

$(SIGNATURES)
Returns the real discharge in Pitztal before loss
"""
function realDischarge(Discharge_after_loss, loss_parameter)
      #loss = loss_parameter .* (Discharge_after_loss + loss).^2
      loss_all = Float64[]
      for Q in Discharge_after_loss
            p = - 1/ loss_parameter .+ 2 .* Q
            q = Q .^2
            if ((p./2).^2 .- q) < 0
                  loss = 12.1
            else
                  loss = - p./2 .+ sqrt.((p./2).^2 .- q)
            end
            append!(loss_all, loss)
      end
      #loss[loss.>12.1] .= 12.1
      return loss_all + Discharge_after_loss
end

@testset "loss" begin
    for i in 1:10
          loss_parameter = rand(0.01:0.0001:0.08)[1]
          Discharge = rand(0:0.1:50)[]
          Discharge2 = rand(0:0.1:50)[]
          Discharge_after_loss = [Discharge, Discharge2] .- loss([Discharge, Discharge2], loss_parameter)
          println(Discharge, " ",Discharge2, " ",loss_parameter)
          @test  Discharge - 0.00000001 <= realDischarge(Discharge_after_loss, loss_parameter)[1] <= Discharge + 0.00000001
           @test  2 - 0.00000001 <= realDischarge(Discharge_after_loss, loss_parameter)[2] <= Discharge2 + 0.00000001
    end
end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
# Plot of various Mathieu functions, inspired by https://dlmf.nist.gov/28.3
using Mathieu, PyPlot

z = LinRange(0,2π,101)
q = LinRange(0,10,101)

## ce₀(z,q)
f = figure(num=1, figsize=(10,8)); f.clf();
plot_surface(q,z,mapreduce(q->Mathieu.ce(0,q,z),hcat,q),cmap="coolwarm")
title("ce₀(z,q)"); xlabel("q"); ylabel("z"); zlim(0,1.5);
gca(projection="3d").view_init(azim=-30); display(f);

## ce₁(z,q)
f = figure(num=2, figsize=(10,8)); f.clf();
plot_surface(q,z,mapreduce(q->Mathieu.ce(1,q,z),hcat,q),cmap="coolwarm")
title("ce₁(z,q)"); xlabel("q"); ylabel("z"); zlim(-1.2,1.2);
gca(projection="3d").view_init(azim=-30); display(f);

## se₁(z,q)
f = figure(num=3, figsize=(10,8)); f.clf();
plot_surface(q,z,mapreduce(q->Mathieu.se(1,q,z),hcat,q),cmap="coolwarm")
title("se₁(z,q)"); xlabel("q"); ylabel("z"); zlim(-1.5,1.5);
gca(projection="3d").view_init(azim=-30); display(f);

## se₂(z,q)
f = figure(num=4, figsize=(10,8)); f.clf();
plot_surface(q,z,mapreduce(q->Mathieu.se(2,q,z),hcat,q),cmap="coolwarm")
title("se₂(z,q)"); xlabel("q"); ylabel("z"); zlim(-1.2,1.2);
gca(projection="3d").view_init(azim=-30); display(f);

## ce₂(z,q)
f = figure(num=5, figsize=(10,8)); f.clf();
plot_surface(q,z,mapreduce(q->Mathieu.ce(2,q,z),hcat,q),cmap="coolwarm")
title("ce₂(z,q)"); xlabel("q"); ylabel("z"); zlim(-1,1.1);
gca(projection="3d").view_init(azim=-30); display(f);


text:
function matgen(n)
  tmp = 1.0 / n / n
  [ tmp * (i - j) * (i + j - 2) for i=1:n, j=1:n ]
end

function main(n)
  t = time()
  n = round(Int, n / 2 * 2)
  a = matgen(n)
  b = matgen(n)
  c = a * b
  v = round(Int, n/2) + 1
  println(c[v, v])
  println(time() - t)
end

function test()
  n = 100
  if length(ARGS) >= 1
    n = parse(Int, ARGS[1])
  end

  println("JIT warming up")
  main(200)

  println("bench")
  x = @timed main(n)
  println("Elapsed: $(x[2]), Allocated: $(x[3]), GC Time: $(x[4])")
end

test()

max_stars_repo_name: wookay/benchmarks

max_issues_repo_name: wookay/benchmarks

max_forks_repo_name: wookay/benchmarks


text:
# Various term manipulation utilities

"""
    dequantify(term::Term, domain::Domain, state::State)

Replaces universally or existentially quantified expressions with their
corresponding conjuctions or disjunctions over the object types they quantify.
"""
function dequantify(term::Term, domain::Domain, state::State)
    if term.name in (:forall, :exists)
        typeconds, query = flatten_conjs(term.args[1]), term.args[2]
        query = dequantify(query, domain, state)
        ground_terms = [query]
        for cond in typeconds
            type, var = cond.name, cond.args[1]
            ground_terms = map(ground_terms) do gt
                return [substitute(gt, Subst(var => obj))
                        for obj in get_objects(domain, state, type)]
            end
            ground_terms = reduce(vcat, ground_terms)
        end
        op = term.name == :forall ? :and : :or
        return Compound(op, ground_terms)
    elseif term.name in (:and, :or, :imply, :not, :when)
        args = Term[dequantify(arg, domain, state) for arg in term.args]
        return Compound(term.name, args)
    else
        return term
    end
end

"""
    conds, effects = flatten_conditions(term::Term)

Flattens a (potentially nested) conditional effect `term` into a list of
condition lists (`conds`) and a list of effect lists (`effects`).
"""
function flatten_conditions(term::Term)
    if term.name == :when # Conditional case
        cond, effect = term.args[1], term.args[2]
        subconds, subeffects = flatten_conditions(effect)
        for sc in subconds prepend!(sc, flatten_conjs(cond)) end
        return subconds, subeffects
    elseif term.name == :and # Conjunctive case
        conds, effects = [Term[]], [Term[]]
        for effect in term.args
            subconds, subeffects = flatten_conditions(effect)
            for (c, e) in zip(subconds, subeffects)
                if isempty(c)
                    append!(effects[1], e)
                else
                    push!(conds, c)
                    push!(effects, e)
                end
            end
        end
        if isempty(effects[1])
            popfirst!(conds)
            popfirst!(effects)
        end
        return conds, effects
    else # Base case
        cond, effect = Term[], Term[term]
        return [cond], [effect]
    end
end

"""
    actions = flatten_conditions(action::GroundAction)

Flattens ground actions with conditional effects into multiple ground actions.
"""
function flatten_conditions(action::GroundAction)
    if (action.effect isa GenericDiff) return [action] end
    actions = GroundAction[]
    for (conds, diff) in zip(action.effect.conds, action.effect.diffs)
        preconds = [action.preconds; conds]
        push!(actions, GroundAction(action.name, action.term, preconds, diff))
    end
    return actions
end

"""
    clauses = to_cnf_clauses(term)

Convert a `term` to CNF and return a list of the clauses as terms. Clauses with
multiple literals have an `or` around them, but single-literal clauses do not.
"""
function to_cnf_clauses(term::Term)
    term = to_cnf(term)
    clauses = map!(c -> (length(c.args) == 1) ? c.args[1] : c,
                   similar(term.args), term.args)
    return clauses
end
to_cnf_clauses(terms::AbstractVector{<:Term}) =
    isempty(terms) ? Term[] : reduce(vcat, (to_cnf_clauses(t) for t in terms))

"""
    clauses = to_dnf_clauses(term)

Convert a `term` to DNF and return a list of the clauses as terms. Clauses with
multiple literals have an `and` around them, but single-literal clauses do not.
"""
function to_dnf_clauses(term::Term)
    term = to_dnf(term)
    clauses = map!(c -> (length(c.args) == 1) ? c.args[1] : c,
                   similar(term.args), term.args)
    return clauses
end
to_dnf_clauses(terms::AbstractVector{<:Term}) =
    isempty(terms) ? Term[] : reduce(vcat, (to_dnf_clauses(t) for t in terms))


text:
struct BathymetryModel{Td<:AbstractMvNormal, Tω<:AbstractFisheryDomain} <: Any
    D::Td
    Ω::Tω
end
## Might be worth making a type system for spatial kernels...
function BathymetryModel(Ω::GriddedFisheryDomain{Tf, Ti},
                         μ::AbstractVecOrMat{Tf},
                         kern::K) where {Tf<:Real, Ti, K<:AbstractCovarianceKernel}
    length(μ) == length(Ω) || throw(DimensionMismatch("Need mean for every grid cell."))
    D = MvNormal(vec(μ), cov(kern, Ω))
    BathymetryModel(D, Ω)
end

"""
    Bathymetry
        bathymetry::Matrix
        Ω::AbstractFisheryDomain

A realization of a `BathymetryModel`. Typically created using the
`rand(BM::BathyModel)` method.
"""
struct Bathymetry{T<:Real, Tω<:AbstractFisheryDomain} <: Any
    bathymetry::Matrix{T}
    Ω::Tω

    function Bathymetry(bathymetry::Matrix{T}, Ω::Tω) where {T<:Real, Tω<:AbstractFisheryDomain}
        all(size(Ω) .== size(bathymetry)) ||
            throw(DimensionMismatch("Bathymetry must conform to domain grid"))
        new{T, Tω}(bathymetry, Ω)
    end
end

function rand(rng::Random.AbstractRNG, BM::BathymetryModel)
    bathy = rand(BM.D)
    bathy .-= minimum(bathy)
    Bathymetry(reshape(bathy, size(BM.Ω)...), BM.Ω)
end
rand(BM::BathymetryModel) = rand(Random.GLOBAL_RNG, BM)

getindex(B::Bathymetry, i, j) = B.bathymetry[i, j]
getindex(B::Bathymetry, i) = B.bathymetry[i]


text:
# Parameter estimation
using Test
using CircStats

@testset "Distance" begin
    @test cdist(0, 0) == 0
    @test cdist(0, 1) == 1
    @test cdist(-1, 0) == 1
    @test cdist(0, -1) == -1
    @test cdist(2.1, 2.0) ≈ -0.1
    @test cdist(340.1, 0.1, true) ≈ 20.0
    let (a, b) = rand(2)
        @test cdist(a, b, false) == cdist(a, b)
        @test cdist(a, b) ≈ deg2rad(cdist(rad2deg(a), rad2deg(b), true))
        @test cdist(a, b, true) ≈ rad2deg(cdist(deg2rad(a), deg2rad(b)))
    end
end

@testset "Parameter est" begin
    # Example 1.1 of Mardia & Jupp (2000)
    let data = [43, 45, 52, 61, 75, 88, 88, 279, 357], μ = 51, R̄ = 0.711, med = 52
        @test cmean(data, true) ≈ μ atol=0.1
        @test cresultant(data, true) ≈ R̄ atol=0.001
        @test cmedian(data, true) ≈ med
        @test cvariance(data, true) ≈ 1 - R̄ atol=0.001
    end
end


text:
@testset "minimal_is_empty_test" begin
    @test isempty(emptyinterval())
    @test !isempty(interval(-Inf,+Inf))
    @test !isempty(interval(1.0,2.0))
    @test !isempty(interval(-1.0,2.0))
    @test !isempty(interval(-3.0,-2.0))
    @test !isempty(interval(-Inf,2.0))
    @test !isempty(interval(-Inf,0.0))
    @test !isempty(interval(-Inf,-0.0))
    @test !isempty(interval(0.0,Inf))
    @test !isempty(interval(-0.0,Inf))
    @test !isempty(interval(-0.0,0.0))
    @test !isempty(interval(0.0,-0.0))
    @test !isempty(interval(0.0,0.0))
    @test !isempty(interval(-0.0,-0.0))
end

@testset "minimal_is_entire_test" begin
    @test !isentire(emptyinterval())
    @test isentire(interval(-Inf,+Inf))
    @test !isentire(interval(1.0,2.0))
    @test !isentire(interval(-1.0,2.0))
    @test !isentire(interval(-3.0,-2.0))
    @test !isentire(interval(-Inf,2.0))
    @test !isentire(interval(-Inf,0.0))
    @test !isentire(interval(-Inf,-0.0))
    @test !isentire(interval(0.0,Inf))
    @test !isentire(interval(-0.0,Inf))
    @test !isentire(interval(-0.0,0.0))
    @test !isentire(interval(0.0,-0.0))
    @test !isentire(interval(0.0,0.0))
    @test !isentire(interval(-0.0,-0.0))
end

@testset "minimal_equal_test" begin
    @test ==(interval(1.0,2.0), interval(1.0,2.0))
    @test !==(interval(1.0,2.1), interval(1.0,2.0))
    @test ==(emptyinterval(), emptyinterval())
    @test !==(emptyinterval(), interval(1.0,2.0))
    @test ==(interval(-Inf,+Inf), interval(-Inf,+Inf))
    @test !==(interval(1.0,2.4), interval(-Inf,+Inf))
    @test ==(interval(1.0,Inf), interval(1.0,Inf))
    @test !==(interval(1.0,2.4), interval(1.0,Inf))
    @test ==(interval(-Inf,2.0), interval(-Inf,2.0))
    @test !==(interval(-Inf,2.4), interval(-Inf,2.0))
    @test ==(interval(-2.0,0.0), interval(-2.0,0.0))
    @test ==(interval(-0.0,2.0), interval(0.0,2.0))
    @test ==(interval(-0.0,-0.0), interval(0.0,0.0))
    @test ==(interval(-0.0,0.0), interval(0.0,0.0))
    @test ==(interval(0.0,-0.0), interval(0.0,0.0))
end

@testset "minimal_subset_test" begin
    @test issubset(emptyinterval(), emptyinterval())
    @test issubset(emptyinterval(), interval(0.0,4.0))
    @test issubset(emptyinterval(), interval(-0.0,4.0))
    @test issubset(emptyinterval(), interval(-0.1,1.0))
    @test issubset(emptyinterval(), interval(-0.1,0.0))
    @test issubset(emptyinterval(), interval(-0.1,-0.0))
    @test issubset(emptyinterval(), interval(-Inf,+Inf))
    @test !issubset(interval(0.0,4.0), emptyinterval())
    @test !issubset(interval(-0.0,4.0), emptyinterval())
    @test !issubset(interval(-0.1,1.0), emptyinterval())
    @test !issubset(interval(-Inf,+Inf), emptyinterval())
    @test issubset(interval(0.0,4.0), interval(-Inf,+Inf))
    @test issubset(interval(-0.0,4.0), interval(-Inf,+Inf))
    @test issubset(interval(-0.1,1.0), interval(-Inf,+Inf))
    @test issubset(interval(-Inf,+Inf), interval(-Inf,+Inf))
    @test issubset(interval(1.0,2.0), interval(1.0,2.0))
    @test issubset(interval(1.0,2.0), interval(0.0,4.0))
    @test issubset(interval(1.0,2.0), interval(-0.0,4.0))
    @test issubset(interval(0.1,0.2), interval(0.0,4.0))
    @test issubset(interval(0.1,0.2), interval(-0.0,4.0))
    @test issubset(interval(-0.1,-0.1), interval(-4.0, 3.4))
    @test issubset(interval(0.0,0.0), interval(-0.0,-0.0))
    @test issubset(interval(-0.0,-0.0), interval(0.0,0.0))
    @test issubset(interval(-0.0,0.0), interval(0.0,0.0))
    @test issubset(interval(-0.0,0.0), interval(0.0,-0.0))
    @test issubset(interval(0.0,-0.0), interval(0.0,0.0))
    @test issubset(interval(0.0,-0.0), interval(-0.0,0.0))
end

@testset "minimal_less_test" begin
    @test ≤(emptyinterval(), emptyinterval())
    @test !≤(interval(1.0,2.0), emptyinterval())
    @test !≤(emptyinterval(), interval(1.0,2.0))
    @test ≤(interval(-Inf,+Inf), interval(-Inf,+Inf))
    @test !≤(interval(1.0,2.0), interval(-Inf,+Inf))
    @test !≤(interval(0.0,2.0), interval(-Inf,+Inf))
    @test !≤(interval(-0.0,2.0), interval(-Inf,+Inf))
    @test !≤(interval(-Inf,+Inf), interval(1.0,2.0))
    @test !≤(interval(-Inf,+Inf), interval(0.0,2.0))
    @test !≤(interval(-Inf,+Inf), interval(-0.0,2.0))
    @test ≤(interval(0.0,2.0), interval(0.0,2.0))
    @test ≤(interval(0.0,2.0), interval(-0.0,2.0))
    @test ≤(interval(0.0,2.0), interval(1.0,2.0))
    @test ≤(interval(-0.0,2.0), interval(1.0,2.0))
    @test ≤(interval(1.0,2.0), interval(1.0,2.0))
    @test ≤(interval(1.0,2.0), interval(3.0,4.0))
    @test ≤(interval(1.0,3.5), interval(3.0,4.0))
    @test ≤(interval(1.0,4.0), interval(3.0,4.0))
    @test ≤(interval(-2.0,-1.0), interval(-2.0,-1.0))
    @test ≤(interval(-3.0,-1.5), interval(-2.0,-1.0))
    @test ≤(interval(0.0,0.0), interval(-0.0,-0.0))
    @test ≤(interval(-0.0,-0.0), interval(0.0,0.0))
    @test ≤(interval(-0.0,0.0), interval(0.0,0.0))
    @test ≤(interval(-0.0,0.0), interval(0.0,-0.0))
    @test ≤(interval(0.0,-0.0), interval(0.0,0.0))
    @test ≤(interval(0.0,-0.0), interval(-0.0,0.0))
end

@testset "minimal_precedes_test" begin
    @test precedes(emptyinterval(), interval(3.0,4.0))
    @test precedes(interval(3.0,4.0), emptyinterval())
    @test precedes(emptyinterval(), emptyinterval())
    @test !precedes(interval(1.0,2.0), interval(-Inf,+Inf))
    @test !precedes(interval(0.0,2.0), interval(-Inf,+Inf))
    @test !precedes(interval(-0.0,2.0), interval(-Inf,+Inf))
    @test !precedes(interval(-Inf,+Inf), interval(1.0,2.0))
    @test !precedes(interval(-Inf,+Inf), interval(-Inf,+Inf))
    @test precedes(interval(1.0,2.0), interval(3.0,4.0))
    @test precedes(interval(1.0,3.0), interval(3.0,4.0))
    @test precedes(interval(-3.0, -1.0), interval(-1.0,0.0))
    @test precedes(interval(-3.0, -1.0), interval(-1.0,-0.0))
    @test !precedes(interval(1.0,3.5), interval(3.0,4.0))
    @test !precedes(interval(1.0,4.0), interval(3.0,4.0))
    @test !precedes(interval(-3.0, -0.1), interval(-1.0,0.0))
    @test precedes(interval(0.0,0.0), interval(-0.0,-0.0))
    @test precedes(interval(-0.0,-0.0), interval(0.0,0.0))
    @test precedes(interval(-0.0,0.0), interval(0.0,0.0))
    @test precedes(interval(-0.0,0.0), interval(0.0,-0.0))
    @test precedes(interval(0.0,-0.0), interval(0.0,0.0))
    @test precedes(interval(0.0,-0.0), interval(-0.0,0.0))
end

@testset "minimal_interior_test" begin
    @test isinterior(emptyinterval(), emptyinterval())
    @test isinterior(emptyinterval(), interval(0.0,4.0))
    @test !isinterior(interval(0.0,4.0), emptyinterval())
    @test isinterior(interval(-Inf,+Inf), interval(-Inf,+Inf))
    @test isinterior(interval(0.0,4.0), interval(-Inf,+Inf))
    @test isinterior(emptyinterval(), interval(-Inf,+Inf))
    @test !isinterior(interval(-Inf,+Inf), interval(0.0,4.0))
    @test !isinterior(interval(0.0,4.0), interval(0.0,4.0))
    @test isinterior(interval(1.0,2.0), interval(0.0,4.0))
    @test !isinterior(interval(-2.0,2.0), interval(-2.0,4.0))
    @test isinterior(interval(-0.0,-0.0), interval(-2.0,4.0))
    @test isinterior(interval(0.0,0.0), interval(-2.0,4.0))
    @test !isinterior(interval(0.0,0.0), interval(-0.0,-0.0))
    @test !isinterior(interval(0.0,4.4), interval(0.0,4.0))
    @test !isinterior(interval(-1.0,-1.0), interval(0.0,4.0))
    @test !isinterior(interval(2.0,2.0), interval(-2.0,-1.0))
end

@testset "minimal_strictly_less_test" begin
    @test <(emptyinterval(), emptyinterval())
    @test !<(interval(1.0,2.0), emptyinterval())
    @test !<(emptyinterval(), interval(1.0,2.0))
    @test <(interval(-Inf,+Inf), interval(-Inf,+Inf))
    @test !<(interval(1.0,2.0), interval(-Inf,+Inf))
    @test !<(interval(-Inf,+Inf), interval(1.0,2.0))
    @test !<(interval(1.0,2.0), interval(1.0,2.0))
    @test <(interval(1.0,2.0), interval(3.0,4.0))
    @test <(interval(1.0,3.5), interval(3.0,4.0))
    @test !<(interval(1.0,4.0), interval(3.0,4.0))
    @test !<(interval(0.0,4.0), interval(0.0,4.0))
    @test !<(interval(-0.0,4.0), interval(0.0,4.0))
    @test !<(interval(-2.0,-1.0), interval(-2.0,-1.0))
    @test <(interval(-3.0,-1.5), interval(-2.0,-1.0))
end

@testset "minimal_strictly_precedes_test" begin
    @test strictprecedes(emptyinterval(), interval(3.0,4.0))
    @test strictprecedes(interval(3.0,4.0), emptyinterval())
    @test strictprecedes(emptyinterval(), emptyinterval())
    @test !strictprecedes(interval(1.0,2.0), interval(-Inf,+Inf))
    @test !strictprecedes(interval(-Inf,+Inf), interval(1.0,2.0))
    @test !strictprecedes(interval(-Inf,+Inf), interval(-Inf,+Inf))
    @test strictprecedes(interval(1.0,2.0), interval(3.0,4.0))
    @test !strictprecedes(interval(1.0,3.0), interval(3.0,4.0))
    @test !strictprecedes(interval(-3.0,-1.0), interval(-1.0,0.0))
    @test !strictprecedes(interval(-3.0,-0.0), interval(0.0,1.0))
    @test !strictprecedes(interval(-3.0,0.0), interval(-0.0,1.0))
    @test !strictprecedes(interval(1.0,3.5), interval(3.0,4.0))
    @test !strictprecedes(interval(1.0,4.0), interval(3.0,4.0))
    @test !strictprecedes(interval(-3.0,-0.1), interval(-1.0,0.0))
end

@testset "minimal_disjoint_test" begin
    @test isdisjoint(emptyinterval(), interval(3.0,4.0))
    @test isdisjoint(interval(3.0,4.0), emptyinterval())
    @test isdisjoint(emptyinterval(), emptyinterval())
    @test isdisjoint(interval(3.0,4.0), interval(1.0,2.0))
    @test !isdisjoint(interval(0.0,0.0), interval(-0.0,-0.0))
    @test !isdisjoint(interval(0.0,-0.0), interval(-0.0,0.0))
    @test !isdisjoint(interval(3.0,4.0), interval(1.0,7.0))
    @test !isdisjoint(interval(3.0,4.0), interval(-Inf,+Inf))
    @test !isdisjoint(interval(-Inf,+Inf), interval(1.0,7.0))
    @test !isdisjoint(interval(-Inf,+Inf), interval(-Inf,+Inf))
end


text:
include("../../../lib/Newton.jl")

module BayesianNewtonSLAB
using ForwardDiff
using ..NewtonMethod

eucLen    = x -> (sum(x.^2.0))^(0.5)
normalize = x -> x / eucLen(x)

function repeat_fill!(to::AbstractArray, fr::AbstractArray)

    len_fr = length(fr)
    println("Repeat:", len_fr, ", ", length(to))
    println(typeof(fr), ", ", typeof(to))
    for i = 1 : length(to)
        to[i] = fr[mod(i-1, len_fr)+1]
    end
    println("Repeat done")
end

function fit(;
    N        :: Integer,
    period   :: Integer,
    beg_t    :: Integer,
    Δt       :: T,
    init_h   :: Array{T},
    init_Q   :: Array{T},
    θ        :: Array{T},
    F        :: Array{T},
    max      :: Integer,
    η        :: T,
    σ_ϵ      :: T,
    σ_Q      :: T,
    σ_h      :: T,
    h_rng    :: Array{T},
    verbose  :: Bool = false,
) where T <: AbstractFloat

    if mod(length(θ), period) != 0 || mod(N, period) !=0
        throw(ArgumentError("Data length should be multiple of [period]"))
    end

    reduced_years = Int(N / period)

    rng1 = collect(beg_t:beg_t+N-1)
    rng2 = rng1 .+ 1

    σ²_ϵ = σ_ϵ ^ 2
    σ²_Q = σ_Q ^ 2
    σ²_h = σ_h ^ 2

    # Extract fixed data
    _F_ph = (F[rng1] + F[rng2]) / 2.0

    _θ    = θ[rng1]
    _θ_p1 = θ[rng2]

    _∂θ∂t_ph = (_θ_p1 - _θ) / Δt
    _θ_ph    = (_θ_p1 + _θ) / 2.0

    x_mem = zeros(T, 2*period)
    x_mem[ 1       :   period] = init_h
    x_mem[period+1 : 2*period] = init_Q

    #=
    calLogPost = function(x)

        local L = 0.0

        h    = repeat(x[1:period],          outer=(reduced_years,))
        Q_ph = repeat(x[period+1:2*period], outer=(reduced_years,))
        h_p1 = circshift(h, -1)

        h_p1 = circshift(h, -1)
        h_ph = (h_p1 + h) / 2.0

        ϵ =  (
            h_ph .* _∂θ∂t_ph -  _F_ph - Q_ph
        )

        L += - ϵ' * ϵ / σ²_ϵ

        # Add prior of h
        _h = h[1:period]
        L += - sum( ((_h .< h_rng[1]) .* (_h .- h_rng[1])).^2 ) / σ²_h
        L += - sum( ((_h .> h_rng[2]) .* (_h .- h_rng[2])).^2 ) / σ²_h

        # Add prior of Q
        L += - sum( (Q_ph[1:period]).^2 ) / σ²_Q

        return L
    end
    =#

    calLogPost = function(x)

        local L = 0.0

        h_ph = repeat(x[1:period],          outer=(reduced_years,))
        Q_ph = repeat(x[period+1:2*period], outer=(reduced_years,))

        ϵ =  (
            h_ph .* _∂θ∂t_ph -  _F_ph - Q_ph
        )

        L += - ϵ' * ϵ / σ²_ϵ

        # Add prior of h
        _h = h_ph[1:period]
        L += - sum( ((_h .< h_rng[1]) .* (_h .- h_rng[1])).^2 ) / σ²_h
        L += - sum( ((_h .> h_rng[2]) .* (_h .- h_rng[2])).^2 ) / σ²_h

        # Add prior of Q
        L += - sum( (Q_ph[1:period]).^2 ) / σ²_Q

        return L
    end

    f_and_∇f = function(x)
        f  = ForwardDiff.gradient(calLogPost, x)
        ∇f = ForwardDiff.hessian( calLogPost, x)

        return f, ∇f
    end

    x_mem[:] = NewtonMethod.fit(;
        f_and_∇f = f_and_∇f,
        η        = η,
        x0       = x_mem,
        max      = max,
        verbose  = verbose
    )

    return x_mem

end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
Base.show(io::IO, d::AbstractDocument) = print(io, "A $(typeof(d))")
Base.show(io::IO, crps::Corpus) = print(io, "A Corpus")
Base.show(io::IO, dtm::DocumentTermMatrix) = print(io, "A DocumentTermMatrix")

function Base.summary(d::AbstractDocument)
    o = ""
    o *= "A $(typeof(d))\n"
    o *= " * Language: $(language(d))\n"
    o *= " * Title: $(title(d))\n"
    o *= " * Author: $(author(d))\n"
    o *= " * Timestamp: $(timestamp(d))\n"

    if typeof(d) <: Union{TokenDocument, NGramDocument}
        o *= " * Snippet: ***SAMPLE TEXT NOT AVAILABLE***"
    else
        sample_text = replace(text(d)[1:min(50, length(text(d)))], r"\s+" => " ")
        o *= " * Snippet: $(sample_text)"
    end
    return o
end

function Base.summary(crps::Corpus)
    n = length(crps.documents)
    n_s = sum(map(d -> typeof(d) == StringDocument{String}, crps.documents))
    n_f = sum(map(d -> typeof(d) == FileDocument, crps.documents))
    n_t = sum(map(d -> typeof(d) == TokenDocument{String}, crps.documents))
    n_ng = sum(map(d -> typeof(d) == NGramDocument{String}, crps.documents))
    o = ""
    o *= "A Corpus with $n documents:\n"
    o *= " * $n_s StringDocument's\n"
    o *= " * $n_f FileDocument's\n"
    o *= " * $n_t TokenDocument's\n"
    o *= " * $n_ng NGramDocument's\n\n"
    o *= "Corpus's lexicon contains $(lexicon_size(crps)) tokens\n"
    o *= "Corpus's index contains $(index_size(crps)) tokens"
    return o
end

function Base.summary(dtm::DocumentTermMatrix)
    n, p = size(dtm.dtm)
    o = "A $n X $p DocumentTermMatrix"
    return o
end

Base.show(io::IO, ::MIME"text/plain", d::AbstractDocument) = print(io, summary(d))
Base.show(io::IO, ::MIME"text/plain", crps::Corpus) = print(io, summary(crps))
Base.show(io::IO, ::MIME"text/plain", dtm::DocumentTermMatrix) = print(io, summary(dtm))


text:
using AmbientForcing
using NetworkDynamics, Distributions
using Graphs
using OrdinaryDiffEq

# This is an example from NetworkDynamics
# https://github.com/PIK-ICoN/NetworkDynamics.jl/blob/main/examples/kuramoto_plasticity.jl
# This only works with NetworkDynamics v0.5.0 or newer

# Further examples can be found in the test folder

const N_plastic = 10 # number of nodes
k = 4  # average degree
g = barabasi_albert(N_plastic, k)

#=
  Berner, Rico, Eckehard Schöll, and Serhiy Yanchuk.
  "Multiclusters in Networks of Adaptively Coupled Phase Oscillators."
  SIAM Journal on Applied Dynamical Systems 18.4 (2019): 2227-2266.
=#

# The coupling function is modeled by a differential algebraic equation with mass matrix
# 0 * de[1] = e[2] * sin(v_s[1] - v_d[1] + α) / N - e[1] is equivalent to
# e[1] = e[2] * sin(v_s[1] - v_d[1] + α) / N
function kuramoto_plastic_edge!(de, e, v_s, v_d, p, t)
    de[1] =  e[2] * sin(v_s[1] - v_d[1] + α) / N_plastic - e[1]
    de[2] = - ϵ * (sin(v_s[1] - v_d[1] + β) + e[2])

    nothing
end

function kuramoto_plastic_vertex!(dv, v, edges, p, t)
    dv .= 0
    for e in edges
        dv .-= e[1]
    end
end

# Global parameters need to be const for type stability
const ϵ = 0.1
const α = .2π
const β = -.95π

# NetworkDynamics Setup
plasticvertex = ODEVertex(f = kuramoto_plastic_vertex!, dim =1)
mass_matrix_plasticedge = zeros(2,2)
mass_matrix_plasticedge[2,2] = 1. # First variables is set to 0

plasticedge = ODEEdge(f = kuramoto_plastic_edge!, dim=2, sym=[:e, :de], coupling=:undirected,mass_matrix = mass_matrix_plasticedge);
kuramoto_plastic! = network_dynamics(plasticvertex, plasticedge, g)

# Ambient Forcing starts here

# Using a random initial condition x0 violates the constraints!
# The constraints are fulfilled when g(x) ≈ 0.
x0_plastic = rand(106)
g_nd = constraint_equations(kuramoto_plastic!)
sum(g_nd(x0_plastic))

# Using zeros as the initial conditions for the ambient forcing algo
x0_plastic = zeros(106)
g_nd = constraint_equations(kuramoto_plastic!)
g_nd(x0_plastic)

# Perturbing all variables at once
Frand = random_force(kuramoto_plastic!, [0.0, 1.0], Uniform)
afoprob = ambient_forcing_problem(kuramoto_plastic!, x0_plastic, 2.0, Frand)
z_new = ambient_forcing(afoprob, x0_plastic, 2.0, Frand)

# As we can see the constraints are not violated!
sum(g_nd(z_new))

# Perturbing only the variables e_22 and de_22
idx = idx_exclusive(kuramoto_plastic!, ["e_22", "de_22"])
Frand = random_force(kuramoto_plastic!, [0.0, 1.0], Uniform, idx)
z_new = ambient_forcing(afoprob, x0_plastic, 2.0, Frand)

# Still the constraints are fulfilled!
sum(g_nd(z_new))


text:
struct RotationMatrix{T}
    TT::T
    TTᵀ::T
end

# ==================================================================================================================
function solveStokes(U, P, gr, Ucartesian,Upolar, g, ρ,
    η, 𝓒, τ, ε, εII,
    F, DoF_U,
    coordinates, RotationMatrices,
    PhaseID, UBC, SF_Stress,
    KKidx, GGidx, MMidx, Δt,
    to;
    solver = :pardiso)

    ∂Ωu = UBC.Ω
    ufix = UBC.vfix
    ifree = UBC.ifree

    # Picard variables
    max_it = 1
    tol = 1e-3

    @timeit to "Stokes" begin

        U0 = deepcopy(U)

        for it in 1:max_it
            @show mean(η.ip)
            KK,GG,MM,Rhs =
                assembly_stokes_cylindric(gr, coordinates, g, ρ, η, 𝓒, εII, PhaseID, KKidx, GGidx, MMidx)

            KK,GG,Rhs = _prepare_matrices(KK, GG, Rhs, RotationMatrices)
            U, Rhs = _apply_bcs(U,KK,Rhs,∂Ωu,ufix)

            if solver == :pardiso
                @timeit to "PCG solver" U, P = StokesPcCG_pardiso(U, P, KK, MM, GG, Rhs, ifree)

            elseif solver == :suitesparse
                @timeit to "PCG solver" U, P = StokesPcCG(U, P, KK, MM, GG, Rhs, ifree)

            end

            U, Ucart, Upolar, Ucartesian =
                updatevelocity2(U, Ucartesian, Upolar, ρ, RotationMatrices.TT, coordinates, gr)

            stress!(F, ε, εII, τ, η, Ucart, gr.e2n, gr.nel, DoF_U, coordinates, SF_Stress, Δt)

            err = mynorm(U .- U0)/mynorm(U0)
            println("Velocity error = $err")
            println("Mean velocity = $(mean(U))")

            if err ≤ tol
                println("Picard iterations converged after $it iterations")
                break
            end

            U0 .= deepcopy(U)
        end
    end

    return Ucartesian, Upolar, U, P, ε, to
end
# ==================================================================================================================

# ==================================================================================================================
function solveStokes(U, P, gr, Ucartesian,Upolar, g, ρ, η, 𝓒, εII,
    coordinates, RotationMatrices,
    PhaseID, UBC,
    KKidx,GGidx,MMidx,
    to;
    solver = :pardiso)

    ∂Ωu = UBC.Ω
    ufix = UBC.vfix
    ifree = UBC.ifree

    @timeit to "Stokes" begin
        KK,GG,MM,Rhs =
           assembly_stokes_cylindric(gr, coordinates, g, ρ, η, 𝓒, εII, PhaseID, KKidx, GGidx, MMidx)

        KK,GG,Rhs = _prepare_matrices(KK, GG, Rhs, RotationMatrices)
        U,Rhs = _apply_bcs(U,KK,Rhs,∂Ωu,ufix)

        if solver == :pardiso
            @timeit to "PCG solver" U, P = StokesPcCG_pardiso(U, P, KK, MM, GG, Rhs, ifree)

        elseif solver == :suitesparse
            @timeit to "PCG solver" U, P = StokesPcCG(U, P, KK, MM, GG, Rhs, ifree)
        end

        U, Ucart, Upolar, Ucartesian =
            updatevelocity2(U, Ucartesian, Upolar, ρ, RotationMatrices.TT, coordinates, gr)
    end

    return Ucartesian, Upolar, U, Ucart, P, to
end
# ==================================================================================================================

@inline initvelocity(nn) = [Point2D{Cartesian}(0.0,0.0) for _ = 1:nn], [Point2D{Polar}(0.0,0.0) for _ = 1:nn]

function remove_wind(r, Uθ::Vector{Point2D{Polar}})
    rmax = maximum(r)
    id = findall(x->x==rmax,r)
    mean([@inbounds Uθ[i].x for i in id])
end

function cartesianvelocity!(U, Ucart, θ, r)
    for (c, i) in enumerate(1:2:(length(U)-1))
        Ucart[i] = U[i+1]*sin(θ[c]) + r[c]*U[i]*cos(θ[c])
        Ucart[i+1] = U[i+1]*cos(θ[c]) - r[c]*U[i]*sin(θ[c])
    end
end

@inline function updatevelocity(U, Ucartesian, Upolar, θ, r, TT)
    nU = length(U)
    Ucart = similar(U)
    # -- Fill tangential and radial velocities
    @inbounds for (c,i) in enumerate(1:2:nU-1)
        Upolar[c].x = U[i]
        Upolar[c].z = U[i+1]
    end

    # -- Remove net rotation
    U_wind = remove_wind(r, Upolar)
    invr = @. (r-1.22)
    @views (@. U[1:2:nU-1] -= U_wind*invr)
    extrema(@. U[1:2:nU-1] - U_wind)
    # a, b = extrema(@views U[1:2:nU-1])
    # c = 0.5*(a+b)
    # @views (@. U[1:2:nU-1] -= c*invr)

    # -- Convert to Cartesian
    Ucart = TT*U

    # -- Fill cartesian velocities
    @inbounds for (c, i) in enumerate(1:2:nU-1)
        Upolar[c].x = U[i]
        # Ucart[i] = Upolar[c].z*sin(θ[c]) + r[c]*Upolar[c].x*cos(θ[c])
        # Ucart[i+1] = Upolar[c].z*cos(θ[c]) - r[c]*Upolar[c].x*sin(θ[c])
        Ucartesian[c].x = Ucart[i]
        Ucartesian[c].z = Ucart[i+1]
    end
    return U, Ucart, Upolar, Ucartesian
end

modulus(v::Vector{T}) where T = √(dot(v, v))

modulus(v::NTuple{N, T}) where {N,T} = √(dot(v, v))

modulus(v::NTuple{2, T}) where {T} = √(v[1]*v[1] + v[2]*v[2])

function modulus(v::Vector{T}, u::Vector{T}) where T
    m = similar(v)
    Threads.@threads for i in eachindex(v)
        @inbounds m[i] = modulus(
            (u[i], v[i])
        )
    end
    m
end

function cross_product(v::Vector{T}, u::Vector{T}) where T
    m = similar(v)
    Threads.@threads for i in eachindex(v)
        @inbounds m[i] = cross_product(
            (v[i], u[i]),
            (v[i], u[i]),
        )
    end
    m
end

"""
    2D cross-product : ||u||⋅||v||⋅sin(θ)
"""
function cross_product(x, v)
    mx, mv = modulus(x), modulus(v) # vector modulus
    θ = acos(dot(x, v)/(mx*mv)) # angle between vectors
    mx*mv*sin(θ) # 2D cross product
end

@inline angular_velocity(x, v) = modulus(v)*sin(x[1])/x[2] # 2D cross product

function v_perp(gr, Ux, Uz)
    v = similar(Ux)
    for i in eachindex(Ux)
        v[i] = cross_product(
            (gr.x[i], gr.z[i]),
            (Ux[i], Uz[i]),
        )
    end
end

function v_perp(gr, U::Vector{Point2D{T}}) where {T}
    v = Vector{Float64}(undef, length(U))
    Threads.@threads for i in eachindex(U)
        @inbounds v[i] = angular_velocity(
            (gr.θ[i], gr.r[i]),
            (U[i].x, U[i].z),
        )
    end
    v
end

function update_velocity_structures(U, Upolar, Ucartesian, TT)
    nU = length(U)
    # Cartesian velocity
    Ucart = TT*U
    Threads.@threads for i in 1:2:nU-1
        @inbounds Upolar[div(i,2)+1].x = U[i]
        @inbounds Upolar[div(i,2)+1].z = U[i+1]
        @inbounds Ucartesian[div(i,2)+1].x = Ucart[i]
        @inbounds Ucartesian[div(i,2)+1].z = Ucart[i+1]
    end
    return Upolar, Ucartesian, Ucart
end

## Nullspace removal
@inline function updatevelocity2(U, Ucartesian, Upolar, ρ, TT, coordinates, gr)

    θ, r = coordinates.θ, coordinates.r

    # Net rotation
    Uθ = @views U[1:2:end]
    angular_momentum = volume_integral(Uθ.*ρ, gr.e2n, θ, r)
    moment = volume_integral(gr.r.*ρ,  gr.e2n, θ, r)
    ω = angular_momentum/moment
    @views @. U[1:2:end] -= ω*gr.r

    # Update velocity
    Upolar, Ucartesian, Ucart = update_velocity_structures(U, Upolar, Ucartesian, TT)

    return U, Ucart, Upolar, Ucartesian
end

@inline function _apply_bcs(U,KK,Rhs,∂Ω,vfix)
    U[∂Ω] = vfix # write prescribed temperature values
    Rhs -= KK[:,∂Ω] * vfix
    return U,Rhs
end

@inline function _prepare_matrices(KK,GG,Fb,RotationMatrices)
    KK = SparseMatrixCSC(Symmetric(KK,:L))
    dropzeros!(KK)
    dropzeros!(GG)
    KK = RotationMatrices.TTᵀ*KK*RotationMatrices.TT
    GG = RotationMatrices.TTᵀ*GG
    Fb = RotationMatrices.TTᵀ*Fb
    KK = SparseMatrixCSC(Symmetric(KK,:L))
    return KK, GG, Fb
end

@inline function rotation_matrix(theta)
    npt   = length(theta)
    # -- memory ALLOCATIONS
    nnz_T = 4 * npt
    KKi_T = Vector{Int64}(undef,nnz_T)
    KKj_T = similar(KKi_T)
    KK1_T = fill(1.0,nnz_T)

    # Diagonal slots: cos(theta). This will use up the first 2*npt slots
    # KKi_T[1:2*npt]      = 1:2*npt
    # KKj_T[1:2*npt]      = 1:2*npt
    @inbounds for i in 1:2*npt
        KKi_T[i] = i
        KKj_T[i] = i
    end

    idx     = 1:2:2*npt-1
    idx2    = 3*npt+1:4*npt

    @inbounds @fastmath for (ii, i) in enumerate(2*npt+1:3*npt)
        cost           = cos(theta[ii])
        sint           = sin(theta[ii])
        i0             = idx[ii]
        i2             = idx2[ii]
        # Diagonal slots: cos(theta). This will use up the first 2*npt slots
        KK1_T[i0]      = cost   # element (1,1) in the rotation matrix [TT]
        KK1_T[i0+1]    = cost   # element (2,2) in the rotation matrix [TT]
        # Off-diagonal slots: sin(theta) and -sin(theta). This will use up the last 2*npt slots
        KKi_T[i]       = i0
        KKj_T[i]       = i0+1
        KK1_T[i]       = sint  # element (1,2) in the rotation matrix [TT] (ORIGINAL: flipped sign -> transposed rotation)
        # KK1_T[i]       = -sint  # element (1,2) in the rotation matrix [TT]
        KKi_T[i2]      = i0+1
        KKj_T[i2]      = i0
        KK1_T[i2]      = -sint # element (2,1) in the rotation matrix [TT] (ORIGINAL:flipped sign -> transposed rotation)
        # KK1_T[i2]      = sint # element (2,1) in the rotation matrix [TT]
    end

    TT = sparse(KKi_T,KKj_T,KK1_T)
    transTT = SparseMatrixCSC(TT')

    return RotationMatrix(TT, transTT)

end

function assembly_stokes_cylindric(gr, coordinates, g, ρ, η, 𝓒, εII, PhaseID, KKidx, GGidx, MMidx)

    EL2NOD, EL2NODP, theta, r = gr.e2n, gr.e2n_p1, coordinates.θ, coordinates.r

    # ============================================ MODEL AND BLOCKING PARAMETERS
    ndim            = 2
    nnodel          = size(EL2NOD,1)
    nel             = size(EL2NOD,2)
    nelblk          = min(nel, 1200)
    nblk            = ceil(Int,nel/nelblk)
    il              = one(nelblk); iu = nelblk
    # =========== PREPARE INTEGRATION POINTS & DERIVATIVES wrt LOCAL COORDINATES
    nip             = 7
    # nUdof           = ndim*maximum(EL2NOD)
    # nPdof           = maximum(EL2NODP)
    nUdofel         = ndim * nnodel
    nPdofel         = size(EL2NODP,1)
    ni, nn, nnP, nn3 = Val(nip), Val(nnodel), Val(3), Val(3)
    N,dNds,~,w_ip   = _get_SF(ni,nn)
    NP,~,~,~        = _get_SF(ni,nnP)
    N3,~,dN3ds,_    = _get_SF(ni,nn3)
    # ============================================================= ALLOCATIONS
    detJ_PL         = Vector{Float64}(undef,nelblk)
    R_31            = similar(detJ_PL) # =  detJa_PL*dxi_dth
    R_21            = similar(detJ_PL) # = -detJa_PL*deta_dth
    Th_31           = similar(detJ_PL) # = -detJa_PL*dxi_dr
    Th_21           = similar(detJ_PL) # =  detJa_PL*deta_dr
    # NxN             = Matrix{Float64}(undef,3,3)
    # ang_dist        = Matrix{Float64}(undef,3,nelblk)
    # VCOORD_th       = similar(ang_dist)
    # VCOORD_r        = similar(ang_dist)
    dNdx            = Matrix{Float64}(undef,nelblk,size(dNds[1],2))
    dNdy            = similar(dNdx)
    ω               = Vector{Float64}(undef,Int(nelblk))
    invJx_double    = Matrix{Float64}(undef,Int(nelblk),ndim)    # storage for x-components of Jacobi matrix
    invJz_double    = similar(invJx_double)                     # storage for z-components of Jacobi matrix
    # ==================================== STORAGE FOR DATA OF ELEMENTS IN BLOCK
    K_blk           = fill(0.0, nelblk, Int(nUdofel*(nUdofel+1)/2))
        # symmetric stiffness matrix, dim=vel dofs, but only upper triangle
    M_blk           = fill(0.0, nelblk, Int(nPdofel*(nPdofel+1)/2))
        # symmetric mass matrix, dim=pressure nodes, but only upper triangle
    G_blk           = fill(0.0, nelblk, nPdofel*nUdofel)
        # asymmetric gradient matrix, vel dofs x pressure nodes
    Fb_blk          = fill(0.0, nelblk, nUdofel)
        # storage for entries in bouyancy force vector

    # ====================================  STORAGE FOR DATA OF ALL ELEMENT MATRICES/VECTORS
    K_all           = Matrix{Float64}(undef, Int(nUdofel*(nUdofel+1)/2),nel)
    M_all           = Matrix{Float64}(undef, Int(nPdofel*(nPdofel+1)/2),nel)
    G_all           = Matrix{Float64}(undef, nUdofel*nPdofel,nel)
    Fb_all          = Matrix{Float64}(undef, nUdofel,nel)

    #=========================================================================
    BLOCK LOOP - MATRIX COMPUTATION
    =========================================================================#
    @inbounds for ib in 1:nblk
        #===========================================================================
        CALCULATE JACOBIAN, ITS DETERMINANT AND INVERSE

        NOTE: For triangular elements with non-curved edges the Jacobian is
              the same for all integration points (i.e. calculated once
              before the integration loop). Further, linear 3-node shape are
              sufficient to calculate the Jacobian.
        ===========================================================================#
        VCOORD_th = view(theta,:,il:iu)
        VCOORD_r = view(r,:,il:iu)

        J_th = gemmt(VCOORD_th', dN3ds')
        J_r = gemmt(VCOORD_r', dN3ds')

        _fill_R_J!(J_th,J_r,VCOORD_th,VCOORD_r,
            R_31,R_21,Th_31,Th_21,detJ_PL)

        # ------------------------ NUMERICAL INTEGRATION LOOP (GAUSS QUADRATURE)
        _ip_loop!(K_blk, M_blk, G_blk, Fb_blk,
            g, η, 𝓒, ρ, εII,
            EL2NOD, PhaseID, il:iu,
            dNds, w_ip, nip, nnodel, nelblk, nPdofel,nUdofel,
            dNdx,dNdy,N,N3,NP,
            ω,invJx_double, invJz_double, detJ_PL,
            R_21,R_31,Th_21,Th_31,
            VCOORD_th,VCOORD_r)

        # --------------------------------------- WRITE DATA INTO GLOBAL STORAGE
        @views begin
            K_all[:,il:iu] .= K_blk'
            G_all[:,il:iu] .= G_blk'
            M_all[:,il:iu] .= M_blk'
            Fb_all[:,il:iu] .= Fb_blk'
        end

        # -------------------------------------------------- RESET BLOCK MATRICES
        fill!(K_blk, 0.0)
        fill!(G_blk, 0.0)
        fill!(M_blk, 0.0)
        fill!(Fb_blk, 0.0)

        # -------------------------------- READJUST START, END AND SIZE OF BLOCK
        il += nelblk
        if ib == nblk-1
           # ------------ Account for different number of elements in last block
           nelblk       = nel - il + 1  # new block size
           # --------------------------------- Reallocate at blocks at the edge
           detJ_PL      = Vector{Float64}(undef,nelblk)
           R_31         = similar(detJ_PL)
           R_21         = similar(detJ_PL)
           Th_31        = similar(detJ_PL)
           Th_21        = similar(detJ_PL)
           dNdx         = Matrix{Float64}(undef,nelblk,size(dNds[1],2))
           dNdy         = similar(dNdx)
           K_blk        = fill(0.0, nelblk, Int(nUdofel*(nUdofel+1)/2))
           M_blk        = fill(0.0, nelblk, Int(nPdofel*(nPdofel+1)/2))
           G_blk        = fill(0.0, nelblk, nPdofel*nUdofel)
           Fb_blk       = fill(0.0, nelblk, nUdofel)
           ω            = Vector{Float64}(undef,Int(nelblk))
           invJx_double = Matrix{Float64}(undef,Int(nelblk),ndim)
           invJz_double = similar(invJx_double)
        end
        iu += nelblk

    end # end block loop

    # #===========================================================================
    # ASSEMBLY OF GLOBAL SPARSE MATRICES AND RHS-VECTOR
    # ===========================================================================#
    KK, GG, MM, Fb = _create_stokes_matrices(EL2NOD,nUdofel,ndim,
        KKidx, GGidx, MMidx,
        K_all, G_all, M_all, Fb_all)

    return KK,GG,MM,Fb

end # END OF ASSEMBLY FUNCTION

#===============================================================================
                        BEGINING OF ASSEMBLYE SUBFUNCTIONS
===============================================================================#
@inline function _create_stokes_matrices(EL2NOD,EL2NODP,nUdofel,nPdofel,ndim,
    K_all,G_all,M_all,Fb_all)

    nel     = size(EL2NOD,2)
    indx_j  = Array{Int64}(undef,nUdofel,nUdofel)
    indx_i  = Array{Int64}(undef,nUdofel,nUdofel)
    dummy   = 1:nUdofel
    @inbounds for i in 1:nUdofel
        indx_j[i,:] = dummy
        indx_i[:,i] = dummy'
    end

    # indx_i      = indx_j';
    indx_i      = tril(indx_i); indxx_i = vec(indx_i); filter!(x->x>0,indxx_i)
    indx_j      = tril(indx_j); indxx_j = vec(indx_j); filter!(x->x>0,indxx_j)

    #==========================================================================
    ASSEMBLY OF GLOBAL STIFFNESS MATRIX
    ==========================================================================#
    EL2DOF                    = fill(0,nUdofel, nel)
    EL2DOF[1:ndim:nUdofel,:] .= @. ndim*(EL2NOD-1)+1
    EL2DOF[2:ndim:nUdofel,:] .= @. ndim*(EL2NOD-1)+2
    K_i        = deepcopy(EL2DOF[vec(indxx_i),:])
    K_j        = deepcopy(EL2DOF[vec(indxx_j),:])
    indx       = K_i .< K_j
    tmp        = deepcopy(K_j[indx])
    K_j[indx]  = deepcopy(K_i[indx])
    K_i[indx]  = tmp
    #  convert triplet data to sparse global stiffness matrix (assembly)
    KK         = sparse(vec(K_i) , vec(K_j), vec(K_all))

    #==========================================================================
    ASSEMBLY OF GLOBAL GRADIENT MATRIX
    ==========================================================================#
    G_i    = repeat(EL2DOF,nPdofel,1) # global velocity dofs
    indx_j = repeat(collect(1:nPdofel)',nUdofel,1)
    G_j    = EL2NODP[indx_j,:]        # global pressure nodes
    # convert triplet data to sparse global gradient matrix (assembly)
    GG     = sparse(vec(G_i) , vec(G_j) ,vec(G_all))

    #==========================================================================
    ASSEMBLY OF GLOBAL FORCE VECTOR
    ==========================================================================#
    Fb     = accumarray(vec(EL2DOF),vec(Fb_all));

    #==========================================================================
    ASSEMBLY OF GLOBAL (1/VISCOSITY)-SCALED MASS MATRIX
    ==========================================================================#
    indx_j  = Array{Int64}(undef,nPdofel,nPdofel)
    indx_i  = Array{Int64}(undef,nPdofel,nPdofel)
    dummy   = collect(1:nPdofel)
    @inbounds for i in 1:nPdofel
        indx_j[i,:] = dummy
        indx_i[:,i] = dummy'
    end
    indx_i    = tril(indx_i); indxx_i = vec(indx_i); filter!(x->x>0,indxx_i)
    indx_j    = tril(indx_j); indxx_j = vec(indx_j); filter!(x->x>0,indxx_j)
    M_i       = deepcopy(EL2NODP[indxx_i,:]); MM_i = vec(M_i);
    M_j       = deepcopy(EL2NODP[indxx_j,:]); MM_j = vec(M_j);
    indx      = MM_i .< MM_j
    tmp       = deepcopy(MM_j[indx])
    M_j[indx] = deepcopy(MM_i[indx])
    M_i[indx] = tmp;
    # convert triplet data to sparse global matrix
    MM        = sparse(vec(M_i) , vec(M_j) , vec(M_all))

    return KK,GG,MM,Fb

end

@inline function _create_stokes_matrices(EL2NOD,nUdofel,ndim,
    KKidx,GGidx,MMidx,
    K_all,G_all,M_all,Fb_all)

    nel     = size(EL2NOD,2)

    #==========================================================================
    ASSEMBLY OF GLOBAL STIFFNESS MATRIX
    ==========================================================================#
    #  convert triplet data to sparse global stiffness matrix (assembly)
    KK     = sparse(KKidx._i ,KKidx._j, vec(K_all))

    #==========================================================================
    ASSEMBLY OF GLOBAL GRADIENT MATRIX
    ==========================================================================#
    # convert triplet data to sparse global gradient matrix (assembly)
    GG     = sparse(GGidx._i ,GGidx._j ,vec(G_all))

    #==========================================================================
    ASSEMBLY OF GLOBAL FORCE VECTOR
    ==========================================================================#
    EL2DOF = Array{Int64}(undef,nUdofel, nel)
    @views EL2DOF[1:ndim:nUdofel,:] .= @. ndim*(EL2NOD-1)+1
    @views EL2DOF[2:ndim:nUdofel,:] .= @. ndim*(EL2NOD-1)+2
    Fb     = accumarray(vec(EL2DOF),vec(Fb_all));

    #==========================================================================
    ASSEMBLY OF GLOBAL (1/VISCOSITY)-SCALED MASS MATRIX
    ==========================================================================#
    # convert triplet data to sparse global matrix
    MM     = sparse(MMidx._i ,MMidx._j , vec(M_all))

    return KK,GG,MM,Fb

end

@inline function _fill_R_J!(J_th,J_r,VCOORD_th,VCOORD_r,
    R_31,R_21,Th_31,Th_21,detJ_PL)

    VCOORD_r1 = view(VCOORD_r,1,:)
    VCOORD_r2 = view(VCOORD_r,2,:)
    VCOORD_r3 = view(VCOORD_r,3,:)
    VCOORD_th1 = view(VCOORD_th,1,:)
    VCOORD_th2 = view(VCOORD_th,2,:)
    VCOORD_th3 = view(VCOORD_th,3,:)

    @turbo for i in 1:length(R_31)
        detJ_PL[i] = J_th[i,1]*J_r[i,2] - J_th[i,2]*J_r[i,1]
        R_31[i]    = VCOORD_r3[i]  - VCOORD_r1[i]  # =  detJa_PL*dxi_dth
        R_21[i]    = VCOORD_r2[i]  - VCOORD_r1[i]  # = -detJa_PL*deta_dth
        Th_31[i]   = VCOORD_th3[i] - VCOORD_th1[i] # = -detJa_PL*dxi_dr
        Th_21[i]   = VCOORD_th2[i] - VCOORD_th1[i] # =  detJa_PL*deta_dr
    end

    return R_31,R_21,Th_31,Th_21,detJ_PL
end

@inline function _ip_loop!(K_blk, M_blk, G_blk, Fb_blk,
    g, η, 𝓒, ρ, εII,
    EL2NOD, PhaseID,els,
    dNds, w_ip, nip, nnodel, nelblk, nPdofel,nUdofel,
    dNdx, dNdy, N, N3, NP,
    ω, invJx_double, invJz_double, detJ_PL,
    R_21, R_31, Th_21,Th_31,
    VCOORD_th, VCOORD_r)

    sin_ip  = similar(detJ_PL)
    cos_ip  = similar(detJ_PL)
    th_ip   = similar(detJ_PL)
    r_ip    = similar(detJ_PL)
    NP_blk  = Matrix{Float64}(undef, length(detJ_PL), 3)

    @inbounds for ip in 1:nip

        N_ip = N3[ip]
        NP_blk .= NP[ip].*ones(nelblk)

        #===========================================================================================================
        CALCULATE 2nd JACOBIAN (FROM CARTESIAN TO POLAR COORDINATES --> curved edges), ITS DETERMINANT AND INVERSE

        NOTE: For triangular elements with curved edges the Jacobian needs to be computed at each integration
               point (inside the integration loop).
        ===========================================================================================================#
        gemmt!(th_ip,  VCOORD_th', N_ip')
        gemmt!(r_ip,  VCOORD_r', N_ip')

        #=======================================================================
        PROPERTIES OF ELEMENTS AT ip-TH EVALUATION POINT
        =======================================================================#
        # interpolate density onto integration point
        Dens_blk    = _element_density(ρ, EL2NOD, PhaseID, els, NP[ip])
        # interpolate viscosity onto integration point
        Visc_blk    = _element_viscosity(η, EL2NOD, PhaseID, els, NP[ip], εII, r_ip, ip)
        # Gravitational force at ip-th integration point
        Fg_blk      = g * Dens_blk

        _derivative_weights!(dNds[ip],ω,dNdx,dNdy,w_ip[ip],th_ip,r_ip,sin_ip,cos_ip,
            R_21,R_31,Th_21,Th_31,detJ_PL,invJx_double,invJz_double)

        _fill_Kblk!(K_blk, 𝓒, Visc_blk, Val(η),
            ω, dNdx, dNdy, nnodel, els, ip)

        _fill_Gblk!(G_blk,NP_blk,
            ω, dNdx, dNdy, nPdofel,nUdofel,nnodel)

        _fill_Mblk!(M_blk,Visc_blk,NP_blk,
            ω, nPdofel)

        _fill_Fbblk!(Fb_blk, Fg_blk, N[ip],
            ω, sin_ip,cos_ip,nUdofel)

    end # end integration point loop
end # ENF OF IP_LOOP FUNCTION

@inline function _derivative_weights!(dNds,ω,dNdx,dNdy,w_ip,th_ip,r_ip,sin_ip,cos_ip,
    R_21,R_31,Th_21,Th_31,detJ_PL,invJx_double,invJz_double)

    # --- Start loop (faster and more memory friendly than single - line declarations)
    for i in 1:length(th_ip)
        # -------------------------------- i-th extractions for some extra speed
        detJ      = detJ_PL[i]
        invdetJ   = 1/detJ
        sin_ip0   = sin(th_ip[i])
        cos_ip0   = cos(th_ip[i])
        sin_ip[i] = sin_ip0
        cos_ip[i] = cos_ip0
        R_31_i    = R_31[i]
        R_21_i    = R_21[i]
        Th_31_i   = Th_31[i]
        Th_21_i   = Th_21[i]
        r_i       = r_ip[i]
        invr_i    = 1/r_i
        # --------------------------------------- Inverse of the double Jacobian
        invJx_double[i,1] = ( R_31_i*cos_ip0 - Th_31_i*sin_ip0*r_i )* invr_i * invdetJ
        invJx_double[i,2] = (-R_21_i*cos_ip0 + Th_21_i*sin_ip0*r_i )* invr_i * invdetJ
        invJz_double[i,1] = (-R_31_i*sin_ip0 - Th_31_i*cos_ip0*r_i )* invr_i * invdetJ
        invJz_double[i,2] = ( R_21_i*sin_ip0 + Th_21_i*cos_ip0*r_i )* invr_i * invdetJ
        # ---------------------------- Numerical integration of element matrices
        ω[i]  = r_i * detJ * w_ip
    end

    # --------------------------------------- Derivatives wrt global coordinates
    mul!(dNdx, invJx_double, dNds)
    mul!(dNdy, invJz_double, dNds)

end

@inline function _fill_Kblk!(K_blk, 𝓒, Visc_blk, ::Val{Isotropic}, ω, dNdx, dNdy, nnodel,els,ip)
    C1 = 4/3; # Used instead of 'D' matrix in standard assembly
    C2 = 2/3; # see Zienkiewicz book, Vol 2, 4th edition, p 519

    # C1 = 2; # incompressible
    # C2 = 0; #
    indx = 1
    nn = size(K_blk,1)
    for i in 1:nnodel
        # |-> ugly as fuck, but faster and more memory friendly than MATLAB style semivectorised original version
        for j in i:nnodel
            # x-velocity equation (1st, 3rd, 5th,... rows of stiffness matrices)
            @inbounds @simd for k in 1:nn
                # x-velocity (1st, 3th, 5th,... columns)
                K_blk[k, indx] +=  Visc_blk[k]*ω[k] *
                    ( C1 * dNdx[k,i]*dNdx[k,j] + dNdy[k,i]*dNdy[k,j])
            end
            indx   += 1
            @inbounds @simd for k in 1:nn
                # y-velocity equation (2nd, 4th, 6th,... rows of stiffness matrices)
                K_blk[k, indx] +=  Visc_blk[k]*ω[k] *
                    ( -C2 * dNdx[k,i]*dNdy[k,j] + dNdy[k,i]*dNdx[k,j])
            end
            indx   += 1
        end

        for j in i:nnodel
            if j > i
                # x-velocity equation (3rd, 5th, 7th... rows of stiffness matrices)
                @inbounds @simd for k in 1:nn
                    # x-velocity (1st, 3th, 5th,... columns)
                    K_blk[k, indx] +=  Visc_blk[k]*ω[k] *
                        (-C2 * dNdy[k,i]*dNdx[k,j] + dNdx[k,i]*dNdy[k,j])
                end
                indx   += 1
            end
            @inbounds @simd for k in 1:nn
                # y-velocity equation (2nd, 4th, 6th,... rows of stiffness matrices)
                K_blk[k, indx] += Visc_blk[k]*ω[k] *
                    ( C1 * dNdy[k,i]*dNdy[k,j] + dNdx[k,i]*dNdx[k,j])
            end
            indx += 1
        end

    end

end

@inline function _fill_Kblk!(K_blk, 𝓒, Visc_blk, ::Val{Anisotropic}, ω, dNdx, dNdy, nnodel,els,ip)

    # Unpack viscous tensor
    η11 = view(𝓒.η11,els,ip).*Visc_blk
    η22 = view(𝓒.η33,els,ip).*Visc_blk
    η33 = view(𝓒.η55,els,ip).*Visc_blk
    η12 = view(𝓒.η13,els,ip).*Visc_blk
    η13 = view(𝓒.η15,els,ip).*Visc_blk.*1
    η23 = view(𝓒.η35,els,ip).*Visc_blk.*1

    # C = 0.5; # Used instead of 'D' matrix in standard assembly
    C = 1 # Used instead of 'D' matrix in standard assembly
    indx = 1
    nn = size(K_blk,1)
    for i in 1:nnodel
        # |-> ugly as fuck, but faster and more memory friendly than semivectorised original version
        for j in i:nnodel
            # x-velocity equation (1st, 3rd, 5th,... rows of stiffness matrices)
            @inbounds @simd for k in 1:nn
                # x-velocity (1st, 3th, 5th,... columns)
                K_blk[k,indx] += ω[k] *
                    (dNdx[k,i]*(η11[k]*dNdx[k,j] + C*η13[k]*dNdy[k,j]) +
                     dNdy[k,i]*(η13[k]*dNdx[k,j] + C*η33[k]*dNdy[k,j]))
            end
            indx += 1
            @inbounds @simd for k in 1:nn
                # y-velocity equation (2nd, 4th, 6th,... rows of stiffness matrices)
                K_blk[k,indx] += ω[k] *
                    (dNdx[k,i]*(η12[k]*dNdy[k,j] + C*η13[k]*dNdx[k,j]) +
                     dNdy[k,i]*(η23[k]*dNdy[k,j] + C*η33[k]*dNdx[k,j]))
            end
            indx += 1
        end

        for j in i:nnodel
            if j > i
                # x-velocity equation (3rd, 5th, 7th... rows of stiffness matrices)
                @inbounds @simd for k in 1:nn
                    # x-velocity (1st, 3th, 5th,... columns)
                    K_blk[k,indx] += ω[k] *
                        (dNdx[k,i]*(η13[k]*dNdx[k,j] + C*η33[k]*dNdy[k,j]) +
                         dNdy[k,i]*(η12[k]*dNdx[k,j] + C*η23[k]*dNdy[k,j]))
                end
                indx   += 1
            end
            @inbounds @simd for k in 1:nn
                # y-velocity equation (2nd, 4th, 6th,... rows of stiffness matrices)
                K_blk[k,indx] += ω[k] *
                    (dNdx[k,i]*(η23[k]*dNdy[k,j] + C*η33[k]*dNdx[k,j]) +
                     dNdy[k,i]*(η22[k]*dNdy[k,j] + C*η23[k]*dNdx[k,j]))
            end
            indx += 1
        end

    end

end

# @inline function _fill_Kblk!(K_blk, 𝓒, Visc_blk, ::Val{Anisotropic}, ω, dNdx, dNdy, nnodel,els,ip)

#     # Unpack viscous tensor
#     η11 = view(𝓒.η11,els,ip).*Visc_blk
#     η22 = view(𝓒.η33,els,ip).*Visc_blk
#     η33 = view(𝓒.η55,els,ip).*Visc_blk
#     η12 = view(𝓒.η13,els,ip).*Visc_blk
#     η13 = view(𝓒.η15,els,ip).*Visc_blk.*1
#     η23 = view(𝓒.η35,els,ip).*Visc_blk.*1

#     C1 =  2/3; # Used instead of 'D' matrix in standard assembly
#     C2 = -1/3; # Used instead of 'D' matrix in standard assembly
#     C3 =  0.5; # Used instead of 'D' matrix in standard assembly
#     indx = 1
#     nn = size(K_blk,1)
#     for i in 1:nnodel
#         # |-> ugly as fuck, but faster and more memory friendly than semivectorised original version
#         for j in i:nnodel
#             # x-velocity equation (1st, 3rd, 5th,... rows of stiffness matrices)
#             @inbounds @simd for k in 1:nn
#                 # x-velocity (1st, 3th, 5th,... columns)
#                 K_blk[k,indx] += ω[k] *
#                     dNdx[k,i]*( η11[k]*C1*dNdx[k,j] + η12[k]*C2*dNdx[k,j] + η13[k]*(C3*dNdx[k,j])) +
#                     dNdy[k,i]*( η13[k]*C1*dNdx[k,j] + η23[k]*C2*dNdx[k,j] + η33[k]*(C3*dNdx[k,j]))
#             end
#             indx   += 1
#             @inbounds @simd for k in 1:nn
#                 # y-velocity equation (2nd, 4th, 6th,... rows of stiffness matrices)
#                 K_blk[k,indx] += ω[k] *

#                     (dNdx[k,i]*(2*η12[k]*dNdy[k,j] + η13[k]*dNdx[k,j])/2 +
#                      dNdy[k,i]*(2*η23[k]*dNdy[k,j] + η33[k]*dNdx[k,j])/2)

#                     dNdx[k,i]*( η11[k]*C2*dNdy[k,j] + η12[k]*C1*dNdy[k,j] + η13[k]*(C3*dNdy[k,j])) +
#                     dNdy[k,i]*( η13[k]*C2*dNdy[k,j] + η23[k]*C1*dNdy[k,j] + η33[k]*(C3*dNdy[k,j]))
#             end
#             indx   += 1
#         end

#         for j in i:nnodel
#             if j > i
#                 # x-velocity equation (3rd, 5th, 7th... rows of stiffness matrices)
#                 @inbounds @simd for k in 1:nn
#                     # x-velocity (1st, 3th, 5th,... columns)
#                     K_blk[k,indx] += ω[k] *
#                         dNdx[k,i]*( η13[k]*C1*dNdx[k,j] + η23[k]*C2*dNdx[k,j] + η33[k]*(C3*dNdx[k,j])) +
#                         dNdy[k,i]*( η12[k]*C1*dNdx[k,j] + η22[k]*C2*dNdx[k,j] + η23[k]*(C3*dNdx[k,j]))
#                 end
#                 indx   += 1
#             end
#             @inbounds @simd for k in 1:nn
#                 # y-velocity equation (2nd, 4th, 6th,... rows of stiffness matrices)
#                 K_blk[k,indx] += ω[k] *
#                     dNdx[k,i]*( η13[k]*C2*dNdy[k,j] + η23[k]*C1*dNdy[k,j] + η33[k]*(C3*dNdy[k,j])) +
#                     dNdy[k,i]*( η12[k]*C2*dNdy[k,j] + η22[k]*C1*dNdy[k,j] + η23[k]*(C3*dNdy[k,j]))
#             end
#             indx += 1
#         end

#     end

# end

@inline function _fill_Gblk!(G_blk,NP_blk, ω,dNdx,dNdy, nPdofel,nUdofel,nnodel)
    N = ones(1,nnodel)
    @inbounds  for i in 1:nPdofel
        tmp1 =  ω .* view(NP_blk,:,i)
        # tmp2 = repeat(tmp1,1,nnodel)
        tmp2 = tmp1.*N
        ii = (i-1)*nUdofel .+ (1:2:nUdofel)
        @views G_blk[:,ii] .+= tmp2 .* dNdx
        ii = ii .+ 1
        @views G_blk[:,ii] .+= tmp2 .* dNdy
    end

end

@inline function _fill_Mblk!(M_blk,Visc_blk,NP_blk, ω, nPdofel)
    # ASSEMBLY OF (1/VISCOSITY)-SCALED MASS MATRICES FOR ALL ELEMENTS IN BLOCK
    weight3 = ω./Visc_blk;
    # "weight" divided by viscosity at integration point in all elements
    indx = 1;
    for i in 1:nPdofel
        @inbounds @simd for j in i:nPdofel
            M_blk[:,indx] .+= weight3.*view(NP_blk,:,i).*view(NP_blk,:,j)
            indx += 1
        end
    end
end

@inline function _fill_Fbblk!(Fb_blk,Fg_blk,NU,
    ω, sin_ip,cos_ip,nUdofel)
    # ASSEMBLY OF BUOYANCY FORCE VECTORS FOR ALL ELEMENTS IN BLOCK
    cte = ω .* Fg_blk * NU
    @views begin
        Fb_blk[:,1:2:nUdofel] .+= sin_ip .* cte
        Fb_blk[:,2:2:nUdofel] .+= cos_ip .* cte
    end
end

@inline function _element_viscosity(η, EL2NOD, PhaseID, els, Nip, εII, r_ip, ip)
    idx = view(EL2NOD, 1:3, els)
    Visc_blk = vec(view(η.val,idx)'*Nip')

    # Make sure a column-vector is returned
    return  Visc_blk
end

@inline function _element_viscosity(η::TemperatureDependantPlastic, e2n, PhaseID, els, Nip, εII, r_ip, ip)
    # constant yield stress
    VonMises = η.τ_VonMises
    # slope of yield surface
    slope = 20*VonMises
    # allocate viscosity at the integration point
    η_eff = Vector{Float64}(undef, length(els))

    @inbounds for i in 1:length(els)

        # indices of the element vertices
        idx = @views e2n[1:3, els[i]]
        # interpolate temperature to the ip i.e. T_ip = T_nop ⋅ N
        ηT = η.node[idx[1]]*Nip[1] + η.node[idx[2]]*Nip[2] + η.node[idx[3]]*Nip[3]
        # 2nd invariant of the strain rate at integration point
        εII_ip = εII[els[i], ip]

        # should be false only in the 1st time step
        if εII_ip != 0.0
            # yield stress
            # τy = min(VonMises, (2.22-r_ip[i])*slope)
            τy = (2.22-r_ip[i])*slope
            # 'plastic viscosity'
            ηy = 0.5*τy/εII[els[i], ip]
            # effective viscosity (cap it: 1e-3 ≤ η ≤ 1e3)
            η_eff[i] = η.ip[els[i], ip] = max(min(1/(1/ηT + 1/ηy), 1e3), 1e-3)
            # η_eff[i] = η.ip[els[i], ip] = 1/(1/ηT + 1/ηy)

        else
            η_eff[i] = η.ip[els[i], ip] = max(min(ηT, 1e3), 1e-3)
            # η_eff[i] = η.ip[els[i], ip] = ηT

        end

    end

    return  η_eff
end

_element_viscosity(η::Isoviscous{T}, EL2NOD, PhaseID, els, Nip, εII, r_ip, ip) where{T} = η.val

@inline function _element_density(ρ,EL2NOD,PhaseID,els,Nip)
    idx = view(EL2NOD,1:3,els)
    Dens_blk = vec(view(ρ,idx)'*Nip')

    return  Dens_blk
end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
using PlotlyJS

function tri_edges(verts_matr::AbstractMatrix, tri_matr::AbstractMatrix;
                   linecolor="rgb(70,70,70)", linewidth=1)
    #def triangulation edges from vertices and triangles, of size (m, 3), resp (n,3)
    if minimum(tri_matr) == 0
        tri_verts = [hcat(verts_matr[t, :]) for t in eachrow(tri_matr .+ 1)]
    else
        tri_verts = [hcat(verts_matr[t, :]) for t in eachrow(tri_matr)]
    end
    Xe = []
    Ye = []
    Ze = []
    for T in tri_verts
        append!(Xe, T[:, 1])
        append!(Xe, [T[1,1], NaN])
        append!(Ye, T[:, 2])
        append!(Ye, [T[1,2], NaN])
        append!(Ze, T[:, 3])
        append!(Ze, [T[1,3], NaN])
    end
    return PlotlyJS.scatter3d(
                   x=Xe,
                   y=Ye,
                   z=Ze,
                   mode="lines",
                   name=" ",
                   line=attr(color=linecolor, width=1.5))
end;

function black_bg(; bg_color = "rgb(20, 20, 20)",
                    grid_color="rgb(150, 150, 150)",
                    zeroline=false)
return  Dict{Any, Any}( #black background
        :showbackground => true,
        :backgroundcolor => bg_color,
        :gridcolor => grid_color,
        :zeroline => zeroline)
end;


text:
function KL(p::AbstractMvNormal, q::AbstractMvNormal)
    # This is the generic implementation for AbstractMvNormal, you might need to specialize for your type
    length(p) == length(q) ||
        throw(DimensionMismatch("Distributions p and q have different dimensions $(length(p)) and $(length(q))"))
    Σp, Σq = cov.((p, q))
    Δμ = mean(p) - mean(q)
    0.5 * (tr(Σq \ Σp) + dot(Δμ, Σq \ Δμ) - length(p) + logdet(Σq) - logdet(Σp))
end

function KL(p::MvNormal, q::MvNormal)
    # We use p.Σ and q.Σ to take the advantage that they are defined as PDMats objects
    length(p) == length(q) ||
        throw(DimensionMismatch("Distributions p and q have different dimensions $(length(p)) and $(length(q))"))
    0.5 * (tr(q.Σ \ p.Σ) + invquad(q.Σ, mean(p) - mean(q)) - length(p) + logdet(q.Σ) - logdet(p.Σ))
end