text:
import Base.getindex
import ControlSystems.numeric_type

"""
    G = LQG(sys::AbstractStateSpace, Q1, Q2, R1, R2; qQ=0, qR=0, M = I, N = I)

Return an LQG object that describes the closed control loop around the process `sys=ss(A,B,C,D)`
where the controller is of LQG-type. The controller is specified by weight matrices `Q1,Q2`
that penalizes state deviations and control signal variance respectively, and covariance
matrices `R1,R2` which specify state drift and measurement covariance respectively.

`qQ` and `qR` can be set to incorporate loop transfer recovery, i.e.,
```julia
L = lqr(A, B, Q1+qQ*C'C, Q2)
K = kalman(A, C, R1+qR*B*B', R2)
```
Increasing `qQ` will add more cost in output direction, e.g., encouraging the use of cheap control, while
increasing `qR` adds fictious dynamics noise, makes the observer faster in the direction we control.

`M` is a matrix that defines the controlled variables `z`, i.e., the variables for which you provide reference signals. If no `M` is provided, the default is to consider all state variables of the system as controlled. The definitions of `z` and `y` are given below
```
y = C*x
z = M*x
```
`size(M, 1)` determines the size of the `Q1` matrix you need to supply.

`N` is a matrix that defines how the dynamics noise `v` enters the system, i.e. If no `N` is provided, the default is to consider all state variables being affected by independent noise components. The definition of `v` is given below
```
x′ = A*x + B*u + N*v
```
`size(N, 2)` determines the size of the `R1` matrix you need to supply.

# Example

```julia
s = tf("s")
P = [1/(s+1) 2/(s+2); 1/(s+3) 1/(s-1)]
sys = ss(P)
eye(n) = Matrix{Float64}(I,n,n) # For convinience

qQ = 1
qR = 1
Q1 = 10eye(4)
Q2 = 1eye(2)
R1 = 1eye(6)
R2 = 1eye(2)

G = LQG(sys, Q1, Q2, R1, R2, qQ=qQ, qR=qR)

Gcl = G.cl
T = G.T
S = G.S
sigmaplot([S,T],exp10.(range(-3, stop=3, length=1000)))
stepplot(Gcl)
```

"""
struct LQGProblem
    sys::ExtendedStateSpace
    Q1::AbstractMatrix
    Q2::AbstractMatrix
    R1::AbstractMatrix
    R2::AbstractMatrix
    qQ::Real
    qR::Real
end

ControlSystems.isdiscrete(l::LQGProblem) = ControlSystems.isdiscrete(l.sys)
ControlSystems.iscontinuous(l::LQGProblem) = ControlSystems.iscontinuous(l.sys)

function LQGProblem(
    sys::AbstractStateSpace,
    Q1::AbstractVecOrMat,
    Q2::AbstractVecOrMat,
    R1::AbstractVecOrMat,
    R2::AbstractVecOrMat;
    qQ = 0,
    qR = 0,
    kwargs...,
)
    sys isa ExtendedStateSpace || (sys = ExtendedStateSpace(sys))
    @unpack B1, C1, C2 = sys
    Q1 = Q1 isa AbstractVector ? diagm(Q1) : Q1
    Q2 = Q2 isa AbstractVector ? diagm(Q2) : Q2
    R1 = R1 isa AbstractVector ? diagm(R1) : R1
    R2 = R2 isa AbstractVector ? diagm(R2) : R2
    size(Q1, 1) == size(C1,1) || throw(ArgumentError("The size of Q1 is determined by C1, not by the state."))
    size(R2, 1) == size(C2,1) || throw(ArgumentError("The size of R2 is determined by C2, not by the state."))
    size(R1, 1) == size(B1,2) || throw(ArgumentError("The size of R1 is determined by B1, not by the state."))
    LQGProblem(sys, Q1, Q2, R1, R2, qQ, qR; kwargs...)
end

"""
    LQGProblem(P::ExtendedStateSpace)

If only an `ExtendedStateSpace` system is provided, the system `P` is assumed to correspond to the H₂ optimal control problem with
```
C1'C1    = Q1
D12'D12  = Q2

B1*B1'   = R1
D21*D21' = R2
```
and an `LQGProblem` with the above covariance matrices is returned. The system
description in the returned LQGProblem will have `B1 = C1 = I`.
See Ch. 13 in Robust and optimal control for reference.
"""
function LQGProblem(P::ExtendedStateSpace)
    @unpack A, B1, B2, C1, C2, D11, D12, D21, D22 = P
    # all(iszero, D22) || throw(ArgumentError("Non-zero D22 not handled."))
    all(iszero, D11) || throw(ArgumentError("Non-zero D11 not handled."))
    all(iszero, D12'C1) || throw(ArgumentError("D12'C1 should be 0"))
    all(iszero, D21*B1') || throw(ArgumentError("D21*B1' should be 0"))
    Q1 = C1'C1
    Q2 = D12'D12
    R1 = B1*B1'
    R2 = D21*D21'
    B1 = I(P.nx)
    C1 = I(P.nx)

    P = ss(A, B1, B2, C1, C2; D22, Ts = P.timeevol)
    # P = ss(A, B1, B2, C1, C2; Ts = P.timeevol)
    LQGProblem(P, Q1, Q2, R1, R2)
end

function ControlSystems.kalman(l::LQGProblem)
    @unpack A, C2, B1, R1, qR, B2, R2 = l
    fun = isdiscrete(l) ? dkalman : kalman
    K = fun(A, C2, Hermitian(B1*R1*B1' + qR * B2 * B2'), R2)
end

function ControlSystems.lqr(l::LQGProblem)
    @unpack A, B2, C1, Q1, qQ, C2, Q2 = l
    fun = isdiscrete(l) ? dlqr : lqr
    L = fun(A, B2, Hermitian(C1'Q1*C1 + qQ * C2'C2), Q2)
end

"""
    ControlSystems.lqr(P::AbstractStateSpace, Q1::AbstractMatrix, Q2::AbstractMatrix, Q3::AbstractMatrix)

Calculate the feedback gain of the LQR cost function augmented with control differences
```math
x'Q1*x + u'Q2*u + Δu'Q3*Δu
Δu = u(k) - u(k-1)
```
"""
function ControlSystems.lqr(P::AbstractStateSpace, Q1::AbstractMatrix, Q2::AbstractMatrix, Q3::AbstractMatrix)
    Pd = add_input_differentiator(P)
    S = zeros(Pd.nx, P.nu)
    S[P.nx+1:end, :] = -Q3
	X, _, L = MatrixEquations.ared(Pd.A, Pd.B, Q2+Q3, cat(Q1, Q3, dims=(1,2)), S) # ME has cost matrices reversed
    L[:, 1:P.nx]
end

"""
    ControlSystems.dare(P::AbstractStateSpace, Q1::AbstractMatrix, Q2::AbstractMatrix, Q3::AbstractMatrix)

Solve the discrete-time algebraic Riccati equation for a LQR cost augmented with control differences
```math
x'Q1*x + u'Q2*u + Δu'Q3*Δu
Δu = u(k) - u(k-1)
```
"""
function ControlSystems.dare(P::AbstractStateSpace, Q1::AbstractMatrix, Q2::AbstractMatrix, Q3::AbstractMatrix)
    # The reference cited in MatrixEquations.ared, W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations
    # defines the cost function as x'Q1x + u'Q2u + 2x'Su.
    # The Δu term expands to u+'Q2u + u'Q2u - 2u Q3 u+, so the factor 2 is already accounted for
    Pd = add_input_differentiator(P)
    S = zeros(Pd.nx, P.nu)
    S[P.nx+1:end, :] = -Q3
	X, _, L = MatrixEquations.ared(Pd.A, Pd.B, Q2+Q3, cat(Q1, Q3, dims=(1,2)), S) # ME has cost matrices reversed
    X[1:P.nx, 1:P.nx]
end

function static_gain_compensation(l::LQGProblem, L = lqr(l))
    @unpack A, C1, B1, B2, D11, D12  = l
    pinv(D12 - (C1 - D12*L) * inv(A - B2*L) * B2)
end

function static_gain_compensation(A, B, C, D, L)
    pinv(D - (C - D*L) * inv(A - B*L) * B) # if D is nonzero
    # pinv(C * ((B * L - A) \ B))
end

"""
    extended_controller(K::AbstractStateSpace)

Takes a controller and returns an `ExtendedStateSpace` version which has augmented input `[r; y]` and output `y` (`z` output is 0-dim).
"""
function extended_controller(K::AbstractStateSpace)
    nx,nu,ny = K.nx, K.nu, K.ny
    A,B,C,D = ssdata(K)
    @error("This has not been verified")

    B1 = zeros(nx, nx) # dynamics not affected by r
    B2 = B # input y
    D21 = C#   K*r
    C2 = -C # - K*y
    C1 = zeros(0, nx)
    ss(A, B1, B2, C1, C2; D21, Ts = K.timeevol)
end

"""
    extended_controller(l::LQGProblem, L = lqr(l), K = kalman(l))

Returns an expression for the controller that is obtained when state-feedback `u = -L(xᵣ-x̂)` is combined with a Kalman filter with gain `K` that produces state estimates x̂. The controller is an instance of `ExtendedStateSpace` where `C2 = -L, D21 = L` and `B2 = K`.

The returned system has *inputs* `[xᵣ; y]` and outputs the control signal `u`. If a reference model `R` is used to generate state references `xᵣ`, the controller from `e = ry - y -> u` is given by
```julia
Ce = extended_controller(l)
Ce = named_ss(Ce; x = :xC, y = :u, u = [R.y; :y^l.ny]) # Name the inputs of Ce the same as the outputs of `R`.
connect([R, Ce]; u1 = R.y, y1 = R.y, w1 = [:ry^l.ny, :y^l.ny], z1=[:u])
```

Since the negative part of the feedback is built into the returned system, we have
```
C = observer_controller(l)
Ce = extended_controller(l)
system_mapping(Ce) == -C
````
"""
function extended_controller(l::LQGProblem, L = lqr(l), K = kalman(l))
    A,B,C,D = ssdata(system_mapping(l))
    Ac = A - B*L - K*C + K*D*L # 8.26b
    nx = l.nx
    B1 = zeros(nx, nx) # dynamics not affected by r
    B2 = K # input y
    D21 = L #   L*xᵣ # should be D21?
    C2 = -L # - L*x̂
    C1 = zeros(0, nx)
    ss(Ac, B1, B2, C1, C2; D21, Ts = l.timeevol)
end

"""
    observer_controller(l::LQGProblem, L = lqr(l), K = kalman(l))

Returns an expression for the feedback controller `u = Cy` that is obtained when state-feedback `u = -Lx̂` is combined with a Kalman filter with gain `K` that produces state estimates x̂.

Note: the transfer function returned is only a representation of the controller in the simple setting described above, e.g., it is not valid if the actual input contains anything that is not produced by a pure feedback from observed states. To obtain a controller that takes references into account, see `extended_controller`.
"""
function ControlSystems.observer_controller(l::LQGProblem, L::AbstractMatrix = lqr(l), K::AbstractMatrix = kalman(l))
    A,B,C,D = ssdata(system_mapping(l))
    Ac = A - B*L - K*C + K*D*L # 8.26b
    Bc = K
    Cc = L
    Dc = 0
    ss(Ac, Bc, Cc, Dc, l.timeevol)
end

function ff_controller(l::LQGProblem, L = lqr(l), K = kalman(l))
    Ae,Be,Ce,De = ssdata(system_mapping(l))
    Ac = Ae - Be*L - K*Ce + K*De*L # 8.26c
    Bc = Be*static_gain_compensation(l, L)
    Cc = L
    Dc = 0
    return 1 - ss(Ac, Bc, Cc, Dc, l.timeevol)
end

"""
    closedloop(l::LQGProblem, L = lqr(l), K = kalman(l))

Closed-loop system as defined in Glad and Ljung eq. 8.28

The return value will be the closed loop from reference only, other disturbance signals (B1) are ignored. See `feedback` for a more advanced option.
"""
function closedloop(l::LQGProblem, L = lqr(l), K = kalman(l))
    # todo: reimplement as lft
    P = system_mapping(l)
    @unpack A, B2, C2, C1 = l
    n = P.nx
    # Lr = pinv(C1 * ((P.B * L[:, 1:n] - P.A) \ P.B))
    Lr = static_gain_compensation(l, L[:, 1:n])
    # Lr = (D - (C - D*L) * inv(A - B*L) * B)
    if any(!isfinite, Lr) || all(iszero, Lr)
        @warn "Could not compensate for static gain automatically." Lr
        Lr = 1
    end
    Acl = [A-B2*L B2*L; zero(A) A-K*C2] # 8.28
    BLr = B2 * Lr # QUESTION: should be B1 here? Glad Ljung has B2
    Bcl = [BLr; zero(BLr)]
    Ccl = [C1 zero(C1)]
    syscl = ss(Acl, Bcl, Ccl, 0, l.timeevol)
end

# function closedloop(l::LQGProblem)
#     K = observer_controller(l)
#     Ke = extended_controller(K)
#     lft(l.sys, Ke)
# end

system_mapping(l::LQGProblem) = system_mapping(l.sys)

performance_mapping(l::LQGProblem) = performance_mapping(l.sys)

function Base.getproperty(G::LQGProblem, s::Symbol)
    if s ∈ fieldnames(LQGProblem)
        return getfield(G, s)
    end
    return getproperty(G.sys, s)
    error("No property named $s")
end

sensdoc = """
```
         ▲
         │e₁
         │  ┌─────┐
d₁────+──┴──►  P  ├─────┬──►e₄
      │-    └─────┘     │
      │                 │
      │     ┌─────┐     │
 e₂◄──┴─────┤  C  ◄──┬──+───d₂
            └─────┘  │
                     │e₃
                     ▼
```
- Input  sensitivity is the transfer function from d₁ to e₁,                 (I + CP)⁻¹
- Output sensitivity is the transfer function from d₂ to e₂,                 (I + PC)⁻¹
- Input  complementary sensitivity is the transfer function from d₁ to e₂, CP(I + CP)⁻¹
- Output complementary sensitivity is the transfer function from d₂ to e₄, PC(I + PC)⁻¹
"""

"""
See [`output_sensitivity`](@ref)
$sensdoc
"""
function sensitivity(args...)# Sensitivity function
    return output_sensitivity(args...)
end

"""
See [`output_comp_sensitivity`](@ref)
$sensdoc
"""
function comp_sensitivity(args...) # Complementary sensitivity function
    return output_comp_sensitivity(args...)
end

G_PS(P, C) = P*input_sensitivity(P, C)
function G_PS(l::LQGProblem) # Load disturbance to output
    return system_mapping(l) * input_sensitivity(l)
end

G_CS(P, C) = C*output_sensitivity(P, C)
function G_CS(l::LQGProblem) # Noise to control signal
    return observer_controller(l) * output_sensitivity(l)
end

loopgain(P,C) = P*C
function loopgain(l::LQGProblem)
    return system_mapping(l)*observer_controller(l)
end

function returndifference(l::LQGProblem)
    PC = loopgain(L)
    p = size(l.C2, 1)
    return ss(Matrix{numeric_type(PC)}(I, p, p), l.timeevol) + PC
end

function stabilityrobustness(l::LQGProblem)
    PC = loopgain(L)
    p = size(l.C2, 1)
    return ss(Matrix{numeric_type(PC)}(I, p, p), l.timeevol) + inv(PC)
end

"""
    input_sensitivity(P, C)
    input_sensitivity(l::LQGProblem)

Transfer function from load disturbance to total plant input.
- "Input" signifies that the transfer function is from the input of the plant.
$sensdoc
"""
function input_sensitivity(P,C)
    T = feedback(C * P)
    ss(I(noutputs(T)), P.timeevol) - T
end
input_sensitivity(l::LQGProblem) = input_sensitivity(system_mapping(l), observer_controller(l))

"""
    output_sensitivity(P, C)
    output_sensitivity(l::LQGProblem)

Transfer function from measurement noise / reference to control signal.
- "output" signifies that the transfer function is from the output of the plant.
$sensdoc
"""
function output_sensitivity(P,C)
    PC = P*C
    S = feedback(ss(Matrix{numeric_type(PC)}(I, ninputs(PC), ninputs(PC)), P.timeevol), PC)
    S.C .*= -1
    S.B .*= -1
    S
end
output_sensitivity(l::LQGProblem) = output_sensitivity(system_mapping(l), observer_controller(l))

"""
    input_comp_sensitivity(P,C)
    input_comp_sensitivity(l::LQGProblem)

Transfer function from load disturbance to control signal.
- "Input" signifies that the transfer function is from the input of the plant.
- "Complimentary" signifies that the transfer function is to an output (in this case controller output)
$sensdoc
"""
function input_comp_sensitivity(P,C)
    T = feedback(C * P)
end
input_comp_sensitivity(l::LQGProblem) = input_comp_sensitivity(system_mapping(l), observer_controller(l))

"""
    output_comp_sensitivity(P,C)
    output_comp_sensitivity(l::LQGProblem)

Transfer function from measurement noise / reference to plant output.
- "output" signifies that the transfer function is from the output of the plant.
- "Complimentary" signifies that the transfer function is to an output (in this case plant output)
$sensdoc
"""
function output_comp_sensitivity(P,C)
    S = output_sensitivity(P,C)
    ss(I(noutputs(S)), P.timeevol) - S
end
output_comp_sensitivity(l::LQGProblem) = output_comp_sensitivity(system_mapping(l), observer_controller(l))

plot(G::LQGProblem) = gangoffourplot(G)

function gangoffour(l::LQGProblem)
    sensitivity(l), G_PS(l), G_CS(l), comp_sensitivity(l)
end

function ControlSystems.gangoffourplot(l::LQGProblem, args...; sigma = true, kwargs...)
    S,D,N,T = gangoffour(l)
    bp = (args...; kwargs...) -> sigma ? sigmaplot(args...; kwargs...) : bodeplot(args...; plotphase=false, kwargs...)
    f1 = bp(S, args...; show=false, title="S = 1/(1+PC)", kwargs...)
    Plots.hline!([1], l=(:black, :dash), primary=false)
    try
        mag, freq = hinfnorm(S)
        isfinite(mag) && isfinite(freq) && (freq>0) && Plots.scatter!([freq], [mag], label="Mₛ = $(round(mag, digits=2))")
    catch
    end
    f2 = bodeplot(D, args...; show=false, title="D = P/(1+PC)", plotphase=false, kwargs...)
    Plots.hline!([1], l=(:black, :dash), primary=false)
    f3 = bodeplot(N, args...; show=false, title="N = C/(1+PC)", plotphase=false, kwargs...)
    f4 = bp(T, args...; show=false, title="T = PC/(1+PC)", kwargs...)
    Plots.hline!([1], l=(:black, :dash), primary=false)
    try
        mag, freq = hinfnorm(T)
        isfinite(mag) && isfinite(freq) && (freq>0) && Plots.scatter!([freq], [mag], label="Mₜ = $(round(mag, digits=2))")
    catch
    end
    Plots.plot(f1,f2,f3,f4)
end

"""
    gangoffourplot2(P, C, args...; sigma = true, kwargs...)

Plot the gang-of-four. If the closed-loop is MIMO, output sensitivity functions will be shown.

$sensdoc
"""
function gangoffourplot2(P, C, args...; sigma = true, kwargs...)
    S,D,N,T = gangoffour2(P,C)
    bp = (args...; kwargs...) -> sigma ? sigmaplot(args...; kwargs...) : bodeplot(args...; plotphase=false, kwargs...)
    f1 = bp(S, args...; show=false, title="S = 1/(1+PC)", kwargs...)
    # Plots.hline!([1], l=(:black, :dash), primary=false)
    Plots.hline!([1.0 1.1 1.2], l=(:dash, [:green :orange :red]), sp=1, lab=["1.0" "1.1" "1.2"], ylims=(-3,1.8))
    f2 = bodeplot(D, args...; show=false, title="D = P/(1+PC)", plotphase=false, kwargs...)
    Plots.hline!([1], l=(:black, :dash), primary=false)
    f3 = bodeplot(N, args...; show=false, title="N = C/(1+PC)", plotphase=false, kwargs...)
    f4 = bp(T, args...; show=false, title="T = PC/(1+PC)", ylims=(-3,1.8), kwargs...)
    Plots.hline!([1], l=(:black, :dash), primary=false)
    Plots.plot(f1,f2,f3,f4, ticks=:default, ylabel="", legend=:bottomright)
end

function gangofsevenplot(P, C, F, args...; sigma = true, ylabel="", layout=4, kwargs...)
    S,D,CS,T = gangoffour2(P,C)
    RY = T*F
    RU = CS*F
    RE = S*F
    bp = (args...; kwargs...) -> sigma ? sigmaplot!(args...; kwargs...) : bodeplot!(args...; plotphase=false, kwargs...)
    Plots.plot(; layout, ticks=:default, xscale=:log10, link=:both)
    bp(S, args...; show=false, title="S = 1/(1+PC)", sp=1, kwargs...)
    # Plots.hline!([1], l=(:black, :dash), primary=false, sp=1)
    Plots.hline!([1.0 1.1 1.2], l=(:dash, [:green :orange :red]), sp=1, lab=["1.0" "1.1" "1.2"], ylims=(-3,1.8))
    bodeplot!(D, args...; show=false, title="PS = P/(1+PC)", plotphase=false, sp=2, kwargs...)
    Plots.hline!([1], l=(:black, :dash), primary=false, sp=2)
    bodeplot!(CS, args...; show=false, title="CS = C/(1+PC)", plotphase=false, sp=3, kwargs...)
    Plots.hline!([1], l=(:black, :dash), primary=false, sp=3)
    bp(T, args...; show=false, title="T = PC/(1+PC)", sp=4, kwargs...)
    Plots.hline!([1], l=(:black, :dash), primary=false, sp=4)

    bp(RE, args...; show=false, title="S = 1/(1+PC)", lab="SF = r->e", l=(:dash, :red), sp=1, kwargs...)
    bp(RY, args...; show=false, title="T = PC/(1+PC)", lab="TF = r->y", l=(:dash, :green), sp=4, kwargs...)
    bp(RU, args...; show=false, title="CS = C/(1+PC)", lab="CSF = r->u", l=(:dash, :blue), sp=3, kwargs...)

    # Plots.plot(f1,f2,f3,f4,f7,f6,f5; ticks=:default, layout)
end

"""
    S, PS, CS, T = gangoffour2(P::LTISystem, C::LTISystem)

Return the gang-of-four. If the closed-loop is MIMO, output sensitivity functions will be shown.

$sensdoc
"""
function gangoffour2(P::LTISystem, C::LTISystem)
    S = output_sensitivity(P, C)
    PS = G_PS(P, C)
    CS = G_CS(P, C)
    T = output_comp_sensitivity(P, C)
    return S, PS, CS, T
end


text:
function FDiv3Vec!(F,cCG,v1CG,v2CG,v3CG,CG,Param)
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;
nz=Param.Grid.nz;
dXdxIF = Param.cache.dXdxIF;
dXdxIC = Param.cache.dXdxIC;
# Contravariant components

vCon = Param.CacheC1
DvCon = Param.CacheC2

vCon .= (v1CG.*Param.dXdxIC11 .+ v2CG.*Param.dXdxIC12) .* cCG;
mul!(reshape(DvCon,OP,OP*NF*nz),CG.DS,reshape(vCon,OP,OP*nz*NF))
F .= F .- DvCon

vCon .= (v1CG.*Param.dXdxIC21 .+ v2CG.*Param.dXdxIC22) .* cCG;
mul!(reshape(PermutedDimsArray(DvCon,(2,1,3,4)),OP,OP*NF*nz),CG.DS,reshape(PermutedDimsArray(vCon,(2,1,3,4)),OP,OP*nz*NF))
F .= F .- DvCon

@views vCon[:,:,:,1:end-1] .= 0.5.*((v1CG[:,:,:,1:end-1] .+ v1CG[:,:,:,2:end]).* Param.dXdxIF31 .+
      (v2CG[:,:,:,1:end-1] .+ v2CG[:,:,:,2:end]).* Param.dXdxIF32) .+
       v3CG[:,:,:,2:end-1] .* Param.dXdxIF33;
@views vCon[:,:,:,1:end-1] .= 0.5.*(cCG[:,:,:,1:end-1] .+ cCG[:,:,:,2:end]).*vCon[:,:,:,1:end-1];
if nz>1
  @views DvCon[:,:,:,1] .= 0.5.*vCon[:,:,:,1];
  @views DvCon[:,:,:,2:end-1] .= 0.5.*(vCon[:,:,:,2:end-1] .- vCon[:,:,:,1:end-2]);
  @views DvCon[:,:,:,end] .= -0.5.*vCon[:,:,:,end-1]; # 0.5 Metric
else
  @views DvCon[:,:,:,1] .= 0
end
F .= F .- DvCon

end

function FDiv3Vec(cCG,v1CG,v2CG,v3CG,CG,Param)
# noch uberarbeiten
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;
nz=Param.Grid.nz;
# Contravariant components
v1Con=v1CG.*
    Param.dXdxIC[:,:,:,:,1,1] +
   v2CG.*
    Param.dXdxIC[:,:,:,:,1,2];
v2Con=v1CG.*
    Param.dXdxIC[:,:,:,:,2,1] +
   v2CG.*
    Param.dXdxIC[:,:,:,:,2,2];

v3Con=0.5*((v1CG[:,:,:,1:end-1]+v1CG[:,:,:,2:end]).*
        Param.dXdxIF[:,:,:,2:nz,3,1]+
      (v2CG[:,:,:,1:end-1]+v2CG[:,:,:,2:end]).*
        Param.dXdxIF[:,:,:,2:nz,3,2]) +
       v3CG[:,:,:,2:nz].*
        Param.dXdxIF[:,:,:,2:nz,3,3];

Dv1Con=reshape(
  CG.DS*reshape(v1Con.*cCG,OP,OP*NF*nz),
   OP,OP,NF,nz);
Dv2Con=permute(
  reshape(
  CG.DS*reshape(
  permute(v2Con.*cCG,
  [2,1,3,4]),
  OP,OP*NF*nz),
  OP,OP,NF,nz),
  [2,1,3,4]);
v3Con=0.5*(cCG[:,:,:,1:nz-1]+cCG[:,:,:,2:nz]).*v3Con;
Dv3Con=zeros(OP,OP,NF,nz);
if nz>1
  Dv3Con[:,:,:,1]=0.5*v3Con[:,:,:,1];
  Dv3Con[:,:,:,2:nz-1]=0.5*(v3Con[:,:,:,2:nz-1]-v3Con[:,:,:,1:nz-2]);
  Dv3Con[:,:,:,nz]=-0.5*v3Con[:,:,:,nz-1]; # 0.5 Metric
end
divCG=Dv1Con+Dv2Con+Dv3Con;

return divCG
end

ext: jl


text:
using BinaryBuilder, Pkg
using BinaryBuilder: download_source, sourcify, init_jll_package

# Read in input `.json` file
json = String(read(ARGS[1]))
buff = IOBuffer(strip(json))
objs = []
while !eof(buff)
    push!(objs, BinaryBuilder.JSON.parse(buff))
end

# Merge the multiple outputs into one
merged = BinaryBuilder.merge_json_objects(objs)
BinaryBuilder.cleanup_merged_object!(merged)

# Download all sources
download_source.(merged["sources"]; verbose=true)

# Also initialize JLL package directories
src_name = merged["name"]
code_dir = joinpath(Pkg.devdir(), "$(src_name)_jll")
deploy_repo = "JuliaBinaryWrappers/$(src_name)_jll.jl"
# Always start from a clean slate
rm(code_dir; recursive=true, force=true)

# Then export platforms to file
open(ARGS[2], "w") do io
    println(io, join(triplet.(merged["platforms"]), " "))
end


text:
using Base.Threads

@views function Markers2Cells(p,phase,xc,yc,dx,dy,ncx,ncy)
    weight      = zeros(Float64, (ncx, ncy))
    phase_th    = [similar(phase) for _ = 1:nthreads()] # per thread
    weight_th   = [similar(weight) for _ = 1:nthreads()] # per thread
    @threads for tid=1:nthreads()
        fill!(phase_th[tid] , 0)
        fill!(weight_th[tid], 0)
    end
    chunks = Iterators.partition(1:p.nmark, p.nmark ÷ nthreads())
    @sync for chunk in chunks
        Threads.@spawn begin
            tid = threadid()
            # fill!(phase_th[tid], 0)  # DON'T
            # fill!(weight_th[tid], 0)
            for k in chunk
                if p.phase[k]>=0
                # Get the indices:
                i = p.cellx[k]
                j = p.celly[k]
                # Relative distances
                dxm = 2.0 * abs(xc[i] - p.x[k])
                dym = 2.0 * abs(yc[j] - p.y[k])
                # Increment cell counts
                area = (1.0 - dxm / dx) * (1.0 - dym / dy)
                phase_th[tid][i, j] += p.phase[k] * area
                weight_th[tid][i, j] += area
                end
            end
        end
    end
    phase  .= reduce(+, phase_th)
    weight .= reduce(+, weight_th)
    phase ./= weight
    return
end

@views function CountMarkersPerCell(p,mpc)#,xc,yc,dx,dy,ncx,ncy)
mpc_th   = [similar(mpc) for _ = 1:nthreads()] # per thread
@threads for tid=1:nthreads()
    fill!(mpc_th[tid], 0)
end
chunks = Iterators.partition(1:p.nmark, p.nmark ÷ nthreads())
    @sync for chunk in chunks
        Threads.@spawn begin
            tid = threadid()
            # fill!(mpc_th[tid], 0)
            for k in chunk
                # if p.phase[k]>=0
                i = p.cellx[k]
                j = p.celly[k]
                 # Get the column:
                # dstx = p.x[k] - xc[1];
                # i    = Int64(round(ceil( (dstx/dx) + 0.5)));
                # # Get the line:
                # dsty = p.y[k] - yc[1];
                # j    = Int64(round(ceil( (dsty/dy) + 0.5)));
                # Increment cell counts
                mpc_th[tid][i, j] += 1
                # end
            end
        end
    end
    fill!(mpc,0)
    mpc .= reduce(+, mpc_th)

    # # Count number of marker per cell
    # @avx for j=1:ncy, i=1:ncx
    #     mpc[i,j] = 0.0
    # end
    # @simd for k=1:p.nmark # @avx ne marche pas ici
    #     if (p.phase[k]>=0)
    #         # Get the column:
    #         dstx = p.x[k] - xc[1];
    #         i    = Int64(round(ceil( (dstx/dx) + 0.5)));
    #         # Get the line:
    #         dsty = p.y[k] - yc[1];
    #         j    = Int64(round(ceil( (dsty/dy) + 0.5)));
    #         # Increment cell count
    #         @inbounds mpc[i,j] += 1.0
    #     end
    # end
    return
end

@views function LocateMarkers(p,dx,dy,xc,yc,xmin,xmax,ymin,ymax)
    # Find marker cell indices
    @threads for k=1:p.nmark
        if (p.x[k]<xmin || p.x[k]>xmax || p.y[k]<ymin || p.y[k]>ymax)
            p.phase[k] = -1
        end
        if p.phase[k]>=0
            dstx         = p.x[k] - xc[1]
            i            = ceil(Int, dstx / dx + 0.5)
            dsty         = p.y[k] - yc[1]
            j            = ceil(Int, dsty / dy + 0.5)
            p.cellx[k]   = i
            p.celly[k]   = j
        end
    end
end

@views function ListMarkers( p, ncx, ncy )
liste = hcat([[Int[] for i in 1:ncx] for j in 1:ncy]...)
for k in 1:p.nmark
    if p.phase[k]>=0
        i = p.cellx[k]
        j = p.celly[k]
        push!(liste[i,j], k)
    end
end
return liste
end


text:
module SFFM
import Base: *, size, show, getindex, +, -, setindex!
import Jacobi, LinearAlgebra, SparseArrays
import Plots, StatsBase, KernelDensity

abstract type AbstractModel end

"""
Construct a SFFM model object.

    Model(
        T::Array{Float64,2},
        C::Array{Float64,1},
        r::NamedTuple{(:r, :R)};
        Bounds::Array{<:Real,2} = [-Inf Inf; -Inf Inf],
    )

# Arguments
- `T::Array{Float64,2}`: generator matrix for the CTMC ``φ(t)``
- `C::Array{Float64,1}`: vector of rates ``d/dt X(t)=C[i]`` for ``i=φ(t)``.
- `r::NamedTuple{(:r, :R)}`: rates for the second fluid.
    - `:r(x::Array{Real})`, a function  which takes arrays of x-values and
        returns a row vector of values for each x-value. i.e. `:r([0;1])`
        returns a `2×NPhases` array where the first row contains all the
        ``rᵢ(0)`` and row 2 all the ``rᵢ(1)`` values.
    - `:R(x::Array{Real})`: has the same structure/behaviour as ``:r`` but
        returns the integral of ``:r``. i.e. `Rᵢ(x)=∫ˣrᵢ(y)dy`.
- `Bounds::Array{<:Real,2}`: contains the bounds for the model. The first row
    are the L and R bounds for ``X(t)`` and the second row the bounds for
    ``Y(t)`` (although the bounds for ``Y(t)`` don't actually do anything yet).

# Outputs
- a model object which is a tuple with fields
    - `:T`: as input
    - `:C`: as input
    - `:r`: a named tuple with fields `(:r, :R, :a)`, `:r` and `:R` are as input
        and `:a = abs.(:r)` returns the absolute values of the rates.
    - `Bounds`: as input
)
"""
struct Model <: AbstractModel
    T::Array{<:Real}
    C::Array{<:Real,1}
    r::NamedTuple{(:r, :R, :a)}
    Bounds::Array{<:Real}
end
function Model(
    T::Array{<:Real},
    C::Array{<:Real,1},
    r::NamedTuple{(:r, :R)};
    Bounds::Array{<:Real,2} = [-Inf Inf; -Inf Inf],
    v::Bool = false,
)
    a(x) = abs.(r.r(x))
    r = (r = r.r, R = r.R, a = a)

    v && println("UPDATE: Model object created with fields ", fieldnames(SFFM.Model))
    return Model(
        T,
        C,
        r,
        Bounds,
    )
end
Model() = Model([0],[0],(r=0, R=0, a=0),[0])

struct AugmentedModel <: AbstractModel
    T::Array{<:Real}
    C::Array{<:Real,1}
    r::NamedTuple{(:r, :R, :a)}
    Bounds::Array{<:Real}
end

function _duplicate_zero_states(T::Array{<:Real,2},C::Array{<:Real,1}, r::NamedTuple{(:r, :R)})
    n_0 = sum(C.==0)

    # assign rates for the augmented model
    C_aug = Array{Float64,1}(undef,length(C)+n_0)
    c_zero = 0
    plus_idx = falses(length(C_aug)) # zero states associated with +
    neg_idx = falses(length(C_aug)) # zero states associated with -
    for i in 1:length(C)
        C_aug[i+c_zero] = C[i]
        if C[i] == 0
            C_aug[i+c_zero+1] = C[i]
            plus_idx[i+c_zero] = true
            neg_idx[i+c_zero+1] = true
            c_zero += 1
        end
    end

    # assign second fluid rates
    function r_aug_inner(x)
        out = zeros(length(C_aug))
        out[(C_aug.!=0).|(plus_idx)] = r.r(x)
        out[neg_idx] = r.r(x)[C.==0]
        return out
    end
    function R_aug(x)
        out = zeros(length(C_aug))
        out[(C_aug.!=0).|(plus_idx)] = r.R(x)
        out[neg_idx] = r.R(x)[C.==0]
        return out
    end
    r_aug = (r = r_aug_inner, R = R_aug)

    # assign augmented generator
    c_zero = 0
    T_aug = zeros(length(C_aug),length(C_aug))
    for i in 1:length(C)
        if C[i] == 0
            # duplicate
            T_aug[i+c_zero,(C_aug.!=0).|(plus_idx)] = T[i,:]
            T_aug[i+c_zero+1,(C_aug.!=0).|(neg_idx)] = T[i,:]
            c_zero += 1
        elseif C[i] < 0
            T_aug[i+c_zero,(C_aug.!=0).|(neg_idx)] = T[i,:]
        elseif C[i] > 0
            T_aug[i+c_zero,(C_aug.!=0).|(plus_idx)] = T[i,:]
        end
    end
    return T_aug, C_aug, r_aug
end

function AugmentedModel(model::Model)
    if (any(model.C.==0))
        T_aug, C_aug = _duplicate_zero_states(model.T,model.C)
    else # no zero states, no augmentation needed
       return model
    end
end

"""

    NPhases(model::Model)

the number of states in the state space
"""
NPhases(model::Model) = length(model.C)

phases(model::Model) = 1:NPhases(model)

"""

    modelDicts(model::Model)

input: a Model object

outputs:
     - SDict: a dictionary with keys `"+","-","0","bullet"`
    and corresponding values `findall(model.C .> 0)`, `findall(model.C .< 0)`,
    `findall(model.C .== 0)`, `findall(model.C .!= 0)`, respectively.

     - TDict: a dictionary of submatrices of `T` with keys
    `"ℓm"` with ``ℓ,m∈{+,-,0,bullet}`` and corresponding values
    `model.T[S[ℓ],S[m]]`.
"""
function modelDicts(model::Model)
    nPhases = NPhases(model)
    SDict = Dict{String,Array}("S" => 1:nPhases)
    SDict["+"] = findall(model.C .> 0)
    SDict["-"] = findall(model.C .< 0)
    SDict["0"] = findall(model.C .== 0)
    SDict["bullet"] = findall(model.C .!= 0)

    TDict = Dict{String,Array}("T" => model.T)
    for ℓ in ["+" "-" "0" "bullet"], m in ["+" "-" "0" "bullet"]
        TDict[ℓ*m] = model.T[SDict[ℓ], SDict[m]]
    end

    return SDict, TDict
end

"""

    TDict(model::Model)

"""
function TDict(model::Model)

end

"""

    Mesh

Abstract type representing a mesh for a numerical scheme.
"""
abstract type Mesh end

include("METools.jl")

const PlusMinusIndex = Union{String,Tuple{String,String}}

const PhaseIndex = Union{Tuple{Union{Int64,Colon},Union{Int64,Colon}},Int64,Colon}

const GeneratorIndex = Union{PlusMinusIndex,Tuple{PlusMinusIndex,PhaseIndex}}

const IndexDict = Dict{Tuple{String,Union{Int64,Colon}},BitArray{1}}

const PartitionedGenerator = Dict{GeneratorIndex,SparseArrays.SparseMatrixCSC{Float64,Int64}}

"""

    MakeFil(
        model::SFFM.Model,
        Nodes::Array{<:Real,1},
        )

Construct dict with entries indexing which cells belong to Fᵢᵐ.
"""
function MakeFil(
    model::SFFM.Model,
    Nodes::Array{<:Real,1},
    )
    meshNIntervals = length(Nodes) - 1
    Δtemp = Nodes[2:end] - Nodes[1:end-1]

    Fil = IndexDict()

    ## Construct the sets Fᵐ = ⋃ᵢ Fᵢᵐ, global index for sets of type m
    idxPlus = model.r.r(Nodes[1:end-1].+Δtemp[:]/2).>0
    idxZero = model.r.r(Nodes[1:end-1].+Δtemp[:]/2).==0
    idxMinus = model.r.r(Nodes[1:end-1].+Δtemp[:]/2).<0
    for i in 1:NPhases(model)
        Fil[("+",i)] = idxPlus[:,i]
        Fil[("0",i)] = idxZero[:,i]
        Fil[("-",i)] = idxMinus[:,i]
        if model.C[i] .<= 0
            Fil[("p+",i)] = [model.r.r(model.Bounds[1,1])[i]].>0
            Fil[("p0",i)] = [model.r.r(model.Bounds[1,1])[i]].==0
            Fil[("p-",i)] = [model.r.r(model.Bounds[1,1])[i]].<0
        end
        if model.C[i] .>= 0
            Fil[("q+",i)] = [model.r.r(model.Bounds[1,end])[i]].>0
            Fil[("q0",i)] = [model.r.r(model.Bounds[1,end])[i]].==0
            Fil[("q-",i)] = [model.r.r(model.Bounds[1,end])[i]].<0
        end
    end
    currKeys = keys(Fil)
    for ℓ in ["+", "-", "0"], i = 1:NPhases(model)
        if ((ℓ,i) ∉ currKeys)
            Fil[(ℓ,i)] = falses(meshNIntervals)
        end
        if (("p"*ℓ,i) ∉ currKeys) && (model.C[i] <= 0)
            Fil["p"*ℓ,i] = falses(1)
        end
        if (("p"*ℓ,i) ∉ currKeys) && (model.C[i] > 0)
            Fil["p"*ℓ,i] = falses(0)
        end
        if (("q"*ℓ,i) ∉ currKeys) && (model.C[i] >= 0)
            Fil["q"*ℓ,i] = falses(1)
        end
        if (("q"*ℓ,i) ∉ currKeys) && (model.C[i] < 0)
            Fil["q"*ℓ,i] = falses(0)
        end
    end
    for ℓ in ["+", "-", "0"]
        Fil[(ℓ,:)] = falses(meshNIntervals * NPhases(model))
        Fil[("p"*ℓ,:)] = trues(0)
        Fil[("q"*ℓ,:)] = trues(0)
        for i = 1:NPhases(model)
            idx = findall(Fil[(ℓ,i)]) .+ (i - 1) * meshNIntervals
            Fil[(ℓ,:)][idx] .= true
            Fil[("p"*ℓ,:)] = [Fil[("p"*ℓ,:)]; Fil[("p"*ℓ,i)]]
            Fil[("q"*ℓ,:)] = [Fil[("q"*ℓ,:)]; Fil[("q"*ℓ,i)]]
        end
    end
    return Fil
end

function MakeDict(
    B::Union{Array{<:Real,2},SparseArrays.SparseMatrixCSC{<:Real,Int64}},
    model::Model,
    mesh::Mesh,
    )

    ## Make a Dictionary so that the blocks of B are easy to access
    N₋ = sum(model.C.<=0)
    N₊ = sum(model.C.<=0)

    BDict = PartitionedGenerator()

    ppositions = cumsum(model.C .<= 0)
    qpositions = cumsum(model.C .>= 0)
    for ℓ in ["+", "-", "0"], m in ["+", "-", "0"]
        for i = 1:NPhases(model), j = 1:NPhases(model)
            FilBases = repeat(mesh.Fil[(ℓ,i)]', NBases(mesh), 1)[:]
            pitemp = falses(N₋)
            qitemp = falses(N₊)
            pjtemp = falses(N₋)
            qjtemp = falses(N₊)
            if model.C[i] <= 0
                if length(pitemp) > 0
                    pitemp[ppositions[i]] = mesh.Fil[("p"*ℓ,i)][1]
                end
            end
            if model.C[j] <= 0
                if length(pjtemp) > 0
                    pjtemp[ppositions[j]] = mesh.Fil[("p"*m,j)][1]
                end
            end
            if model.C[i] >= 0
                if length(qitemp) > 0
                    qitemp[qpositions[i]] = mesh.Fil[("q"*ℓ,i)][1]
                end
            end
            if model.C[j] >= 0
                if length(qjtemp) > 0
                    qjtemp[qpositions[j]] = mesh.Fil[("q"*m,j)][1]
                end
            end
            i_idx = [
                pitemp
                falses((i - 1) * TotalNBases(mesh))
                FilBases
                falses(NPhases(model) * TotalNBases(mesh) - i * TotalNBases(mesh))
                qitemp
            ]
            FjmBases = repeat(mesh.Fil[(m, j)]', NBases(mesh), 1)[:]
            j_idx = [
                pjtemp
                falses((j - 1) * TotalNBases(mesh))
                FjmBases
                falses(NPhases(model) * TotalNBases(mesh) - j * TotalNBases(mesh))
                qjtemp
            ]
            BDict[((ℓ, m),(i, j))] = B[i_idx, j_idx]
        end
        # below we need to use repeat(mesh.Fil[ℓ]', NBases(mesh), 1)[:] to
        # expand the index mesh.Fil[ℓ] from cells to all basis function
        FlBases =
            [mesh.Fil["p"*ℓ,:]; repeat(mesh.Fil[ℓ,:]', NBases(mesh), 1)[:]; mesh.Fil["q"*ℓ,:]]
        FmBases =
            [mesh.Fil["p"*m,:]; repeat(mesh.Fil[m,:]', NBases(mesh), 1)[:]; mesh.Fil["q"*m,:]]
        BDict[((ℓ,m),(:,:))] = B[FlBases, FmBases]
    end
    return BDict
end

abstract type Generator <: AbstractArray{Real,2} end

checksquare(A::AbstractArray{<:Any,2}, str::String) = if !(size(A,1)==size(A,2)); throw(DomainError(str*" must be square")); end
checksquare(A::AbstractArray{<:Any,2}) = checksquare(A, "")

struct LazyGenerator  <: Generator
    blocks::Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2}}
    boundary_flux::NamedTuple{(:upper,:lower),Tuple{NamedTuple{(:in,:out),Tuple{Vector{Float64},Vector{Float64}}},NamedTuple{(:in,:out),Tuple{Vector{Float64},Vector{Float64}}}}}
    T::Array{<:Real,2}
    C::Array{<:Real,1}
    Δ::Array{<:Real,1}
    D::Union{Array{Float64,2},LinearAlgebra.Diagonal{Bool,Array{Bool,1}}}
    pmidx::Union{Array{Bool,2},BitArray{2}}
    Fil::IndexDict
    function LazyGenerator(
        blocks::Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2}},
        boundary_flux::NamedTuple{(:upper,:lower),Tuple{NamedTuple{(:in,:out),Tuple{Vector{Float64},Vector{Float64}}},NamedTuple{(:in,:out),Tuple{Vector{Float64},Vector{Float64}}}}},
        T::Array{<:Real,2},
        C::Array{<:Real,1},
        Δ::Array{<:Real,1},
        D::Union{Array{Float64,2},LinearAlgebra.Diagonal{Bool,Array{Bool,1}}},
        pmidx::Union{Array{Bool,2},BitArray{2}},
        Fil::IndexDict,
    )
        s = size(blocks[1])
        for b in 1:4
            checksquare(blocks[b],"blocks")
            !(s == size(blocks[b])) && throw(DomainError("blocks must be the same size"))
        end
        checksquare(T,"T")
        checksquare(D,"D")
        !(s == size(D)) && throw(DomainError("blocks must be the same size as D"))
        checksquare(pmidx,"pmidx")
        !(size(T) == size(pmidx)) && throw(DomainError("pmidx must be the same size as T"))
        !(length(C) == size(T,1)) && throw(DomainError("C must be the same length as T"))

        return new(blocks,boundary_flux,T,C,Δ,D,pmidx,Fil)
    end
end
function LazyGenerator(
    blocks::Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}},
    boundary_flux::NamedTuple{(:in, :out),Tuple{Array{Float64,1},Array{Float64,1}}},
    T::Array{<:Real,2},
    C::Array{<:Real,1},
    Δ::Array{<:Real,1},
    D::Union{Array{Float64,2},LinearAlgebra.Diagonal{Bool,Array{Bool,1}}},
    pmidx::Union{Array{Bool,2},BitArray{2}},
    Fil::IndexDict,
)
    blocks = (blocks[1],blocks[2],blocks[2],blocks[3])
    boundary_flux = (upper = boundary_flux, lower = boundary_flux)
    return LazyGenerator(
        blocks,
        boundary_flux,
        T,
        C,
        Δ,
        D,
        pmidx,
        Fil,
    )
end

function size(B::SFFM.LazyGenerator)
    sz = size(B.T,1)*size(B.blocks[1],1)*length(B.Δ) + sum(B.C.<=0) + sum(B.C.>=0)
    return (sz,sz)
end
size(B::SFFM.LazyGenerator, n::Int) = size(B)[n]

function *(u::AbstractArray{<:Real,2}, B::LazyGenerator)
    output_type = typeof(u)

    sz_u_1 = size(u,1)
    sz_u_2 = size(u,2)
    sz_B_1 = size(B,1)
    sz_B_2 = size(B,2)
    !(sz_u_2 == sz_B_1) && throw(DomainError("Dimension mismatch, u*B, length(u) must be size(B,1)"))
    N₋ = sum(B.C.<=0)
    N₊ = sum(B.C.>=0)
    if output_type <: SparseArrays.SparseMatrixCSC
        v = SparseArrays.spzeros(sz_u_1,sz_B_2)
    else
        v = zeros(sz_u_1,sz_B_2)
    end
    size_delta = length(B.Δ)
    size_blocks = size(B.blocks[1],1)
    size_T = size(B.T,1)

    # boundaries
    # at lower
    v[:,1:N₋] += u[:,1:N₋]*B.T[B.C.<=0,B.C.<=0]
    # in to lower
    idxdown = N₋ .+ ((1:size_blocks).+size_blocks*size_delta*(findall(B.C .<= 0) .- 1)')[:]
    v[:,1:N₋] += u[:,idxdown]*LinearAlgebra.kron(
        LinearAlgebra.diagm(0 => abs.(B.C[B.C.<=0])),
        B.boundary_flux.lower.in/B.Δ[1],
    )
    # out of lower
    idxup = N₋ .+ (size_blocks*size_delta*(findall(B.C .> 0).-1)' .+ (1:size_blocks))[:]
    v[:,idxup] += u[:,1:N₋]*kron(B.T[B.C.<=0,B.C.>0],B.boundary_flux.lower.out')

    # at upper
    v[:,end-N₊+1:end] += u[:,end-N₊+1:end]*B.T[B.C.>=0,B.C.>=0]
    # in to upper
    idxup = N₋ .+ ((1:size_blocks).+size_blocks*size_delta*(findall(B.C .> 0) .- 1)')[:] .+
        (size_blocks*size_delta - size_blocks)
    v[:,end-N₊+1:end] += u[:,idxup]*LinearAlgebra.kron(
        LinearAlgebra.diagm(0 => B.C[B.C.>0]),
        B.boundary_flux.upper.in/B.Δ[end],
    )
    # out of upper
    idxdown = N₋ .+ (size_blocks*size_delta*(findall(B.C .< 0).-1)' .+ (1:size_blocks))[:] .+
        (size_blocks*size_delta - size_blocks)
    v[:,idxdown] += u[:,end-N₊+1:end]*kron(B.T[B.C.>=0,B.C.<0],B.boundary_flux.upper.out')

    # innards
    for i in 1:size_T, j in 1:size_T
        if i == j
            # mult on diagonal
            for k in 1:size_delta
                k_idx = (i-1)*size_blocks*size_delta .+ (k-1)*size_blocks .+ (1:size_blocks) .+ N₋
                for ℓ in 1:size_delta
                    if (k == ℓ+1) && (B.C[i] > 0)
                        ℓ_idx = k_idx .- size_blocks
                        v[:,k_idx] += B.C[i]*(u[:,ℓ_idx]*B.blocks[4])/B.Δ[ℓ]
                    elseif k == ℓ
                        v[:,k_idx] += (u[:,k_idx]*(abs(B.C[i])*B.blocks[2 + (B.C[i].<0)]/B.Δ[ℓ] + B.T[i,j]*LinearAlgebra.I))
                    elseif (k == ℓ-1) && (B.C[i] < 0)
                        ℓ_idx = k_idx .+ size_blocks
                        v[:,k_idx] += abs(B.C[i])*(u[:,ℓ_idx]*B.blocks[1])/B.Δ[ℓ]
                    end
                end
            end
        elseif B.pmidx[i,j]
            # changes from S₊ to S₋ etc.
            for k in 1:size_delta
                for ℓ in 1:size_delta
                    if k == ℓ
                        i_idx = (i-1)*size_blocks*size_delta .+ (k-1)*size_blocks .+ (1:size_blocks) .+ N₋
                        j_idx = (j-1)*size_blocks*size_delta .+ (k-1)*size_blocks .+ (1:size_blocks) .+ N₋
                        v[:,j_idx] += (u[:,i_idx]*(B.T[i,j]*B.D))
                    end
                end
            end
        else
            i_idx = (i-1)*size_blocks*size_delta .+ (1:size_blocks*size_delta) .+ N₋
            j_idx = (j-1)*size_blocks*size_delta .+ (1:size_blocks*size_delta) .+ N₋
            v[:,j_idx] += (u[:,i_idx]*B.T[i,j])
        end
    end
    return v
end

function *(B::LazyGenerator, u::AbstractArray{<:Real,2})
    output_type = typeof(u)
    sz_u_1 = size(u,1)
    sz_u_2 = size(u,2)
    sz_B_1 = size(B,1)
    sz_B_2 = size(B,2)
    !(sz_u_2 == sz_B_1) && throw(DomainError("Dimension mismatch, u*B, length(u) must be size(B,1)"))
    N₋ = sum(B.C.<=0)
    N₊ = sum(B.C.>=0)
    if output_type <: SparseArrays.SparseMatrixCSC
        v = SparseArrays.spzeros(sz_u_1,sz_B_2)
    else
        v = zeros(sz_u_1,sz_B_2)
    end
    size_delta = length(B.Δ)
    size_blocks = size(B.blocks[1],1)
    size_T = size(B.T,1)

    # boundaries
    # at lower
    v[1:N₋,:] += B.T[B.C.<=0,B.C.<=0]*u[1:N₋,:]
    # in to lower
    idxdown = N₋ .+ ((1:size_blocks).+size_blocks*size_delta*(findall(B.C .<= 0) .- 1)')[:]
    v[idxdown,:] += LinearAlgebra.kron(
        LinearAlgebra.diagm(0 => abs.(B.C[B.C.<=0])),
        B.boundary_flux.lower.in/B.Δ[1],
    )*u[1:N₋,:]
    # out of lower
    idxup = N₋ .+ (size_blocks*size_delta*(findall(B.C .> 0).-1)' .+ (1:size_blocks))[:]
    v[1:N₋,:] += kron(B.T[B.C.<=0,B.C.>0],B.boundary_flux.lower.out')*u[idxup,:]

    # at upper
    v[end-N₊+1:end,:] += B.T[B.C.>=0,B.C.>=0]*u[end-N₊+1:end,:]
    # in to upper
    idxup = N₋ .+ ((1:size_blocks).+size_blocks*size_delta*(findall(B.C .> 0) .- 1)')[:] .+
        (size_blocks*size_delta - size_blocks)
    v[idxup,:] += LinearAlgebra.kron(
        LinearAlgebra.diagm(0 => B.C[B.C.>0]),
        B.boundary_flux.upper.in/B.Δ[end],
    )*u[end-N₊+1:end,:]
    # out of upper
    idxdown = N₋ .+ (size_blocks*size_delta*(findall(B.C .< 0).-1)' .+ (1:size_blocks))[:] .+
        (size_blocks*size_delta - size_blocks)
    v[end-N₊+1:end,:] += kron(B.T[B.C.>=0,B.C.<0],B.boundary_flux.upper.out')*u[idxdown,:]

    # innards
    for i in 1:size_T, j in 1:size_T
        if i == j
            # mult on diagonal
            for k in 1:size_delta
                k_idx = (i-1)*size_blocks*size_delta .+ (k-1)*size_blocks .+ (1:size_blocks) .+ N₋
                for ℓ in 1:size_delta
                    if (k == ℓ+1) && (B.C[i] > 0)
                        ℓ_idx = k_idx .- size_blocks
                        v[ℓ_idx,:] += B.C[i]*(B.blocks[4]*u[k_idx,:])/B.Δ[ℓ]
                    elseif k == ℓ
                        v[k_idx,:] += ((abs(B.C[i])*B.blocks[2 + (B.C[i].<0)]/B.Δ[ℓ] + B.T[i,j]*LinearAlgebra.I)*u[k_idx,:])
                    elseif (k == ℓ-1) && (B.C[i] < 0)
                        ℓ_idx = k_idx .+ size_blocks
                        v[ℓ_idx,:] += abs(B.C[i])*(B.blocks[1]*u[k_idx,:])/B.Δ[ℓ]
                    end
                end
            end
        elseif B.pmidx[i,j]
            # changes from S₊ to S₋ etc.
            for k in 1:size_delta
                for ℓ in 1:size_delta
                    if k == ℓ
                        i_idx = (i-1)*size_blocks*size_delta .+ (k-1)*size_blocks .+ (1:size_blocks) .+ N₋
                        j_idx = (j-1)*size_blocks*size_delta .+ (k-1)*size_blocks .+ (1:size_blocks) .+ N₋
                        v[i_idx,:] += (B.T[i,j]*B.D)*u[j_idx,:]
                    end
                end
            end
        else
            i_idx = (i-1)*size_blocks*size_delta .+ (1:size_blocks*size_delta) .+ N₋
            j_idx = (j-1)*size_blocks*size_delta .+ (1:size_blocks*size_delta) .+ N₋
            v[i_idx,:] += B.T[i,j]*u[j_idx,:]
        end
    end
    return v
end

*(B::LazyGenerator, u::LazyGenerator) = SparseArrays.SparseMatrixCSC(B)*u

function show(io::IO, mime::MIME"text/plain", B::SFFM.LazyGenerator)
    if VERSION >= v"1.6"
        show(io, mime, SparseArrays.SparseMatrixCSC(Matrix(LinearAlgebra.I(size(B,1)))*B))
    else
        show(io, mime, Matrix(SparseArrays.SparseMatrixCSC(Matrix(LinearAlgebra.I(size(B,1)))*B)))
    end
end
# show(B::SFFM.LazyGenerator) = show(stdout, B)

function getindex(B::SFFM.LazyGenerator,row::Int,col::Int)
    checkbounds(B,row,col)

    sz_B_1 = size(B,1)

    N₋ = sum(B.C.<=0)
    N₊ = sum(B.C.>=0)

    v = 0.0

    size_delta = length(B.Δ)
    size_blocks = size(B.blocks[1],1)

    if (row ∈ 1:N₋) && (col ∈ 1:N₋)
        v += B.T[B.C.<=0,B.C.<=0][row,col]
    elseif (row ∉ 1:N₋) && (row ∉ (sz_B_1 .+ 1) .- (1:N₊)) && (col ∈ 1:N₋) # in to lower
        idxdown = N₋ .+ ((1:size_blocks).+size_blocks*size_delta*(findall(B.C .<= 0) .- 1)')[:]
        idx = findfirst(x -> x==row, idxdown)
        if (nothing!=idx)
            v += LinearAlgebra.kron(
                LinearAlgebra.diagm(0 => abs.(B.C[B.C.<=0])),
                B.boundary_flux.lower.in/B.Δ[1],
            )[idx,col]
        end
    elseif (row ∈ 1:N₋) && (col ∉ 1:N₋) && (col ∉ sz_B_1 .+ 1 .- (1:N₊)) # out of lower
        idxup = N₋ .+ (size_blocks*size_delta*(findall(B.C .> 0).-1)' .+ (1:size_blocks))[:]
        idx = findfirst(x -> x==col,idxup)
        (nothing!=idx) && (v += kron(B.T[B.C.<=0,B.C.>0],B.boundary_flux.lower.out')[row,idx])
    elseif (row ∈ (sz_B_1 .+ 1) .- (1:N₊)) && (col ∈ (sz_B_1 .+ 1) .- (1:N₊)) # at upper
        v += B.T[B.C.>=0,B.C.>=0][row-sz_B_1+N₊, col-sz_B_1+N₊]
    elseif (row ∉ (sz_B_1 .+ 1) .- (1:N₊)) && (col ∈ (sz_B_1 .+ 1) .- (1:N₊)) # in to upper
        idxup = N₋ .+ ((1:size_blocks).+size_blocks*size_delta*(findall(B.C .> 0) .- 1)')[:] .+
            (size_blocks*size_delta - size_blocks)
        idx = findfirst(x -> x==row, idxup)
        if (nothing!=idx)
            v += LinearAlgebra.kron(
                LinearAlgebra.diagm(0 => B.C[B.C.>0]),
                B.boundary_flux.upper.in/B.Δ[end],
            )[idx, col-sz_B_1+N₊]
        end
    elseif (row ∈ sz_B_1 .+ 1 .- (1:N₊)) && (col ∉ sz_B_1 .+ 1 .- (1:N₊)) # out of upper
        idxdown = N₋ .+ (size_blocks*size_delta*(findall(B.C .< 0).-1)' .+ (1:size_blocks))[:] .+
            (size_blocks*size_delta - size_blocks)
        idx = findfirst(x -> x==col, idxdown)
        (nothing!=idx) && (v += kron(B.T[B.C.>=0,B.C.<0],B.boundary_flux.upper.out')[row-sz_B_1+N₊, idx])
    else
        # innards
        # find if (T⊗I)[row,col] != 0
        row_shift = row-N₋
        col_shift = col-N₋
        if mod(row_shift,size_delta*size_blocks) == mod(col_shift,size_delta*size_blocks)
            i = (row_shift-1)÷(size_delta*size_blocks) + 1
            j = (col_shift-1)÷(size_delta*size_blocks) + 1
            !(B.pmidx[i,j]) && (v += B.T[i,j])
        end
        # for i in 1:size_T, j in 1:size_T
        # phase block
        i = (row_shift-1)÷(size_delta*size_blocks) + 1
        j = (col_shift-1)÷(size_delta*size_blocks) + 1
        if i == j
            # on diagonal
            # find block position [b_r,b_c] in the remaining tridiagonal matrix
            b_r = (row_shift-1 - (i-1)*(size_delta*size_blocks))÷size_blocks + 1
            b_c = (col_shift-1 - (j-1)*(size_delta*size_blocks))÷size_blocks + 1
            if (b_c == b_r+1) && (B.C[i] > 0) && (b_r+1 <= size_delta)
                # find position [r,c] in remaining block
                r = row_shift - (i-1)*(size_delta*size_blocks) - (b_r-1)*size_blocks
                c = col_shift - (j-1)*(size_delta*size_blocks) - (b_c-1)*size_blocks
                v += B.C[i]*B.blocks[4][r,c]/B.Δ[r]
            elseif b_c == b_r
                r = row_shift - (i-1)*(size_delta*size_blocks) - (b_r-1)*size_blocks
                c = col_shift - (j-1)*(size_delta*size_blocks) - (b_c-1)*size_blocks
                v += abs(B.C[i])*B.blocks[2 + (B.C[i].<0)][r,c]/B.Δ[r]
            elseif (b_c == b_r-1) && (B.C[i] < 0) && (b_r-1 > 0)
                r = row_shift - (i-1)*(size_delta*size_blocks) - (b_r-1)*size_blocks
                c = col_shift - (j-1)*(size_delta*size_blocks) - (b_c-1)*size_blocks
                v += abs(B.C[i])*B.blocks[1][r,c]/B.Δ[r]
            end
        elseif B.pmidx[i,j]
            # changes from S₊ to S₋ etc.
            # find block position [b_r,b_c] in the remaining tridiagonal matrix
            b_r = (row_shift-1 - (i-1)*(size_delta*size_blocks))÷size_blocks + 1
            b_c = (col_shift-1 - (j-1)*(size_delta*size_blocks))÷size_blocks + 1
            if b_r == b_c
                r = row_shift - (i-1)*(size_delta*size_blocks) - (b_r-1)*size_blocks
                c = col_shift - (j-1)*(size_delta*size_blocks) - (b_c-1)*size_blocks
                v += B.T[i,j]*B.D[r,c]
            end
        end
    end
    return v
end

# function getindex_correct(B::SFFM.LazyGenerator,row::Int,col::Int)
#     checkbounds(B,row,col)

#     ei = zeros(1,size(B,1))
#     ei[row] = 1
#     return (ei*B)[col]
# end

function getindex(B::SFFM.LazyGenerator,plus_minus_index::PlusMinusIndex,phase_index::PhaseIndex)
    N₋ = sum(B.C.<=0)
    N₊ = sum(B.C.<=0)

    out = []

    ppositions = cumsum(B.C .<= 0)
    qpositions = cumsum(B.C .>= 0)

    ℓ,m = plus_minus_index
    if typeof(phase_index[1])==Colon
        i_range = 1:length(B.C)
    else
        i_range = phase_index[1]
    end
    if typeof(phase_index[2])==Colon
        j_range = 1:length(B.C)
    else
        j_range = phase_index[2]
    end

    for i = i_range, j = j_range
        FilBases = repeat(B.Fil[ℓ,i]', size(B.D,1), 1)[:]
        pitemp = falses(N₋)
        qitemp = falses(N₊)
        pjtemp = falses(N₋)
        qjtemp = falses(N₊)
        if B.C[i] <= 0
            if length(pitemp) > 0
                pitemp[ppositions[i]] = B.Fil["p"*ℓ,i][1]
            end
        end
        if B.C[j] <= 0
            if length(pjtemp) > 0
                pjtemp[ppositions[j]] = B.Fil["p"*m,j][1]
            end
        end
        if B.C[i] >= 0
            if length(qitemp) > 0
                qitemp[qpositions[i]] = B.Fil["q"*ℓ,i][1]
            end
        end
        if B.C[j] >= 0
            if length(qjtemp) > 0
                qjtemp[qpositions[j]] = B.Fil["q"*m,j][1]
            end
        end
        i_idx = [
            pitemp
            falses((i - 1) * size(B.D,1)*length(B.Δ))
            FilBases
            falses(length(B.C) * size(B.D,1)*length(B.Δ) - i * size(B.D,1)*length(B.Δ))
            qitemp
        ]
        FjmBases = repeat(B.Fil[m,j]', size(B.D,1), 1)[:]
        j_idx = [
            pjtemp
            falses((j - 1) * size(B.D,1)*length(B.Δ))
            FjmBases
            falses(length(B.C) * size(B.D,1)*length(B.Δ) - j * size(B.D,1)*length(B.Δ))
            qjtemp
        ]
        if (typeof(phase_index[1]) != Colon) && (typeof(phase_index[2]) != Colon)
            out = B[i_idx, j_idx]
        end
    end
    # we use repeat(mesh.Fil[ℓ]', NBases(mesh), 1)[:] to
    # expand the index mesh.Fil[ℓ] from cells to all basis function
    if (typeof(phase_index[1]) == Colon) && (typeof(phase_index[2]) == Colon)
        FlBases =
            [B.Fil["p"*ℓ,:]; repeat(B.Fil[ℓ,:]', size(B.D,1), 1)[:]; B.Fil["q"*ℓ,:]]
        FmBases =
            [B.Fil["p"*m,:]; repeat(B.Fil[m,:]', size(B.D,1), 1)[:]; B.Fil["q"*m,:]]
        out = B[FlBases, FmBases]
    end

    return out
end

getindex(B::SFFM.LazyGenerator,idx::Tuple{PlusMinusIndex,PhaseIndex}) = getindex(B,idx[1],idx[2])

# this is a bit slow...
function MakeDict(B::SFFM.LazyGenerator)
    BDict = SFFM.PartitionedGenerator()
    for i in 1:length(B.C), j in 1:length(B.C)
        for k in ["+","-","0"], ℓ in ["+","-","0"]
            plus_minus_index = (k,ℓ)
            phase_index = (i,j)
            BDict[plus_minus_index, phase_index] = B[plus_minus_index, phase_index]
        end
    end
    i = :;
    j = :;
    for k in ["+","-","0"], ℓ in ["+","-","0"]
        plus_minus_index = (k,ℓ)
        phase_index = (i,j)
        BDict[plus_minus_index, phase_index] = B[plus_minus_index, phase_index]
    end
    return BDict
end

struct FullGenerator <: Generator
    BDict::PartitionedGenerator
    B::Union{Array{Float64,Int64}, SparseArrays.SparseMatrixCSC{Float64,Int}}
    Fil::IndexDict
end

FullGenerator(BDict::PartitionedGenerator,B::Union{Array{Float64,Int64}, SparseArrays.SparseMatrixCSC{Float64,Int}}) =
    FullGenerator(BDict,B,IndexDict())

size(B::FullGenerator) = size(B.B)
getindex(B::FullGenerator,i::Int,j::Int) = B.B[i,j]
getindex(B::FullGenerator,plus_minus_index::PlusMinusIndex,phase_index::PhaseIndex) =
    B.BDict[plus_minus_index, phase_index]
getindex(B::FullGenerator,idx::Tuple{PlusMinusIndex,PhaseIndex}) =
    getindex(B,idx[1],idx[2])
+(A::AbstractArray{<:Real,2}, B::FullGenerator) = A+B.B
+(B::FullGenerator, A::AbstractArray{<:Real,2}) = B.B+A
+(A::FullGenerator, B::FullGenerator) = A.B+B.B
-(A::AbstractArray{<:Real,2}, B::FullGenerator) = A-B.B
-(B::FullGenerator, A::AbstractArray{<:Real,2}) = B.B-A
-(A::FullGenerator, B::FullGenerator) = A.B-B.B
*(A::AbstractArray{<:Real,2}, B::FullGenerator) = A*B.B
*(B::FullGenerator, A::AbstractArray{<:Real,2}) = B.B*A
*(A::FullGenerator, B::FullGenerator) = A.B*B.B

function show(io::IO, mime::MIME"text/plain", B::FullGenerator)
    if VERSION >= v"1.6"
        show(io, mime, B.B)
        println("\n with partitions")
        show(io, mime, keys(B.BDict))
    else
        show(io, mime, Matrix(B.B))
        println("\n with partitions")
        show(io, mime, keys(B.BDict))
    end
end
# show(B::FullGenerator) = show(stdout, B)

function MakeLazyGenerator(model::Model, mesh::Mesh; v::Bool=false)
    throw(DomainError("Can construct LazyGenerator for DGMesh, FRAPMesh, only"))
end
function MakeFullGenerator(model::Model, mesh::Mesh; v::Bool=false)
    lazy = MakeLazyGenerator(model,mesh; v=v)
    return materialise(lazy)
end

function MakeDict(
    B::Union{Array{<:Real,2},SparseArrays.SparseMatrixCSC{<:Real,Int64}},
    C::Array{<:Real,1},
    order::Int,
    Fil::IndexDict,
    )

    ## Make a Dictionary so that the blocks of B are easy to access
    N₋ = sum(C.<=0)
    N₊ = sum(C.<=0)

    BDict = PartitionedGenerator()

    n_phases = length(C)
    total_n_bases = Int((size(B,1)-N₋-N₊)/n_phases)

    ppositions = cumsum(C .<= 0)
    qpositions = cumsum(C .>= 0)
    for ℓ in ["+", "-", "0"], m in ["+", "-", "0"]
        for i = 1:n_phases, j = 1:n_phases
            FilBases = repeat(Fil[(ℓ,i)]', order, 1)[:]
            pitemp = falses(N₋)
            qitemp = falses(N₊)
            pjtemp = falses(N₋)
            qjtemp = falses(N₊)
            if C[i] <= 0
                if length(pitemp) > 0
                    pitemp[ppositions[i]] = Fil[("p"*ℓ,i)][1]
                end
            end
            if C[j] <= 0
                if length(pjtemp) > 0
                    pjtemp[ppositions[j]] = Fil[("p"*m,j)][1]
                end
            end
            if C[i] >= 0
                if length(qitemp) > 0
                    qitemp[qpositions[i]] = Fil[("q"*ℓ,i)][1]
                end
            end
            if C[j] >= 0
                if length(qjtemp) > 0
                    qjtemp[qpositions[j]] = Fil[("q"*m,j)][1]
                end
            end
            i_idx = [
                pitemp
                falses((i - 1) * total_n_bases)
                FilBases
                falses(n_phases * total_n_bases - i * total_n_bases)
                qitemp
            ]
            FjmBases = repeat(Fil[(m, j)]', order, 1)[:]
            j_idx = [
                pjtemp
                falses((j - 1) * order)
                FjmBases
                falses(n_phases * order - j * order)
                qjtemp
            ]
            BDict[((ℓ, m),(i, j))] = B[i_idx, j_idx]
        end
        # below we need to use repeat(Fil[ℓ]', order, 1)[:] to
        # expand the index Fil[ℓ] from cells to all basis function
        FlBases =
            [Fil["p"*ℓ,:]; repeat(Fil[ℓ,:]', order, 1)[:]; Fil["q"*ℓ,:]]
        FmBases =
            [Fil["p"*m,:]; repeat(Fil[m,:]', order, 1)[:]; Fil["q"*m,:]]
        BDict[((ℓ,m),(:,:))] = B[FlBases, FmBases]
    end
    return BDict
end

function materialise(lzB::LazyGenerator)
    B = SparseArrays.SparseMatrixCSC{Float64,Int}(LinearAlgebra.I(size(lzB,1)))*lzB
    BDict = MakeDict(B,lzB.C,size(lzB.D,1),lzB.Fil)
    # BDict = MakeDict(lzB) # in the interest of speed use the above
    return FullGenerator(BDict,B,lzB.Fil)
end

function MakeQBDidx(model::Model,mesh::Mesh)
    ## Make QBD index
    N₊ = sum(model.C .>= 0)
    N₋ = sum(model.C .<= 0)

    c = N₋
    QBDidx = zeros(Int, NPhases(model) * TotalNBases(mesh) + N₊ + N₋)
    for k = 1:NIntervals(mesh), i = 1:NPhases(model), n = 1:NBases(mesh)
        c += 1
        QBDidx[c] = (i - 1) * TotalNBases(mesh) + (k - 1) * NBases(mesh) + n + N₋
    end
    QBDidx[1:N₋] = 1:N₋
    QBDidx[(end-N₊+1):end] = (NPhases(model) * TotalNBases(mesh) + N₋) .+ (1:N₊)

    return QBDidx
end

include("polynomials.jl")
include("DGBase.jl")
include("Operators.jl")
include("FVM.jl")
include("FRAPApproximation.jl")
include("Distributions.jl")
include("SimulateSFFM.jl")
include("SFM.jl")
include("Plots.jl")

"""
Construct all the DG operators.

    MakeAll(
        model::SFFM.Model,
        mesh::DGMesh;
        approxType::String = "projection"
    )

# Arguments
- `model`: a model object as output from Model
- `mesh`: a Mesh object
- `approxType::String`: (optional) argument specifying how to approximate R (in
    `MakeR()`)

# Output
- a tuple with keys
    - `Matrices`: see `MakeMatrices`
    - `MatricesR`: see `MakeMatricesR`
    - `B`: see `MakeFullGenerator`
    - `D`: see `MakeD`
    - `DR`: see `MakeDR`
"""
function MakeAll(
    model::SFFM.Model,
    mesh::DGMesh;
    approxType::String = "projection"
)

    # Matrices = MakeMatrices(model, mesh)

    B = MakeFullGenerator(model, mesh)
    R = MakeR(model, mesh, approxType = approxType)
    D = MakeD(mesh, B, R)
    return (
        B = B,
        R = R,
        D = D,
    )
end

end

# for k in keys(lzB.Fil), ℓ in keys(lzB.Fil)
#     plus_minus_index = (k[1],ℓ[1])
#     phase_index = (k[2],ℓ[2])
#     display(plus_minus_index); display( phase_index)
# end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
# ---
# title: 559. Maximum Depth of N-ary Tree
# id: problem559
# author: Indigo
# date: 2021-02-16
# difficulty: Easy
# categories: Tree, Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/maximum-depth-of-n-ary-tree/description/>
# hidden: true
# ---
#
# Given a n-ary tree, find its maximum depth.
#
# The maximum depth is the number of nodes along the longest path from the root
# node down to the farthest leaf node.
#
# _Nary-Tree input serialization  is represented in their level order traversal,
# each group of children is separated by the null value (See examples)._
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2018/10/12/narytreeexample.png)
#
#
#
#     Input: root = [1,null,3,2,4,null,5,6]
#     Output: 3
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2019/11/08/sample_4_964.png)
#
#
#
#     Input: root = [1,null,2,3,4,5,null,null,6,7,null,8,null,9,10,null,null,11,null,12,null,13,null,null,14]
#     Output: 5
#
#
#
#
# **Constraints:**
#
#   * The depth of the n-ary tree is less than or equal to `1000`.
#   * The total number of nodes is between `[0, 104]`.
#
#
## @lc code=start
using LeetCode

function max_depth_559(tree::Vector{Vector{Int}})
    isempty(tree) && return 0
    q = Queue{Tuple{Int, Int}}()
    enqueue!(q, (1, 1))
    res = 1
    while !isempty(q)
        frt = dequeue!(q)
        res = frt[2]
        for son in tree[frt[1]]
            enqueue!(q, (son, res + 1))
        end
    end
    res
end
## @lc code=end


text:
module Avro

include("common.jl")
include("schema.jl")
include("io.jl")
include("generic.jl")
include("data_file.jl")

export DataFile,
       Generic,
       Io,
       Schemas

end


text:
using Gomalish
using CxxWrap
using Test
using Random

for T in [Float64, Float32]
    @testset "dot $T" begin
        stdv1 = StdVector(T[1,2,3])
        stdv2 = StdVector(T[-1,-2,-3])
        mv1 = vector{T}(stdv1)
        mv2 = vector{T}(stdv2)
        print_all(mv1)
        print_all(mv2)
        @test Gomalish.dot(mv1, mv2) ≈ -T(14)
        @test mv1[1] ≈ stdv1[1]
        @test mv1[2] ≈ stdv1[2]
        @test mv1[3] ≈ stdv1[3]
    end
end

for T in [Float64, Float32]
    Random.seed!(12345)
    @testset "cg $T" begin
        coo = COO{T}("matrixfiles/sample.mtx")
        print_all(coo)
        A = CRS{T}(coo)
        jl_x = one(T) .+ rand(T, get_row(A))
        jl_b = one(T) .+ rand(T, get_row(A))
        x = vector{T}(StdVector(jl_x))
        b = vector{T}(StdVector(jl_b))

        solver = CG{CRS{T},T}()
        precond = Jacobi{CRS{T},T}()
        set_create_precond(solver, precond)
        set_apply_precond(solver, precond)
        set_tol(solver, 1.0e-12)
        set_maxiter(solver, get_row(A))
        print_all(A)
        println("---")
        print_all(x)
        println("---")
        print_all(b)

        if solver_check(solve(solver, A, x, b))
            return
        end
        jl_A = T[
            2 -1 0
            -1 2 -1
            0 -1 2
        ]
        println("Result by Gomalish/monolish")
        print_all(x)
        println("Result by Julia")
        jl_x = jl_A \ jl_b
        println(jl_x)
        @test x[1] ≈ jl_x[1]
        @test x[2] ≈ jl_x[2]
        @test x[3] ≈ jl_x[3]
    end
end

for T in [Float64, Float32]
    @testset "cg-impl $T" begin
        DIM = 100
    coo = tridiagonal_toeplitz_matrix(DIM, T(11.0), T(-1.0));
        A = CRS{T}(coo);

        println("===== Matrix informatrion =====")
        println("# of rows : ", get_row(A))
        println("# of cols : ", get_col(A))
        println("# of nnz  : ", get_nnz(A))

        # initial x is rand(0~1)
        x = vector{T}(get_row(A), T(0.0), T(1.0))
        # initial b is {1, 1, 1, ...,1}
        b = vector{T}(get_row(A), T(1.0))
        tol::T = 1.0e-12

        r = vector{T}(get_row(A), T(0.0))
        p = vector{T}(get_row(A), T(0.0))
        q = vector{T}(get_row(A), T(0.0))

        matvec(A, x, q)
        sub(b, q, r)
        Gomalish.copy(r, p)

        for iter in 1:get_row(A)
            matvec(A, p, q)
            tmp = Gomalish.dot(r, r)
            α = tmp / Gomalish.dot(p, q);
            axpy(α, p, x)
            axpy(-α, q, r)
            β = Gomalish.dot(r, r) / tmp
            xpay(β, r, p)
            resid = nrm2(r);
            println(iter + 1)
            @show resid
            (resid < tol) && return
            isnan(resid) && return
        end
    end
end

for T in [Float64, Float32]
    @testset "template solver $T" begin
        Random.seed!(12345)
        coo = COO{T}("matrixfiles/sample.mtx")
        A = CRS{T}(coo)
        jl_x = one(T) .+ rand(T, get_row(A))
        jl_b = one(T) .+ rand(T, get_row(A))
        x = vector{T}(StdVector(jl_x))
        b = vector{T}(StdVector(jl_b))
        solver = CG{CRS{T},T}()
        precond = Jacobi{CRS{T},T}()
        set_create_precond(solver, precond)
        set_apply_precond(solver, precond)
        set_tol(solver, 1.0e-12)
        set_maxiter(solver, get_row(A))
        if solver_check(solve(solver, A, x, b))
        return
        end
        print_all(x)
        jl_A = T[
            2 -1 0
            -1 2 -1
            0 -1 2
        ]
        println("Result by Gomalish/monolish")
        print_all(x)
        println("Result by Julia")
        jl_x = jl_A \ jl_b
        println(jl_x)
        @test x[1] ≈ jl_x[1]
        @test x[2] ≈ jl_x[2]
        @test x[3] ≈ jl_x[3]
    end
end

for T in [Float64, Float32]
    @testset "template solver part2 $T" begin
        Random.seed!(12345)
        coo = COO{T}("matrixfiles/sample.mtx")
        A = Dense{T}(coo)
        jl_x = one(T) .+ rand(T, get_row(A))
        jl_b = one(T) .+ rand(T, get_row(A))
        x = vector{T}(StdVector(jl_x))
        b = vector{T}(StdVector(jl_b))
        solver = BiCGSTAB{Dense{T},T}()
        precond = Gomalish.none{Dense{T},T}()
        set_create_precond(solver, precond)
        set_apply_precond(solver, precond)
        set_tol(solver, 1.0e-12)
        set_maxiter(solver, get_row(A))
        if solver_check(solve(solver, A, x, b))
        return
        end
        print_all(x)
        jl_A = T[
            2 -1 0
            -1 2 -1
            0 -1 2
        ]
        println("Result by Gomalish/monolish")
            print_all(x)
        println("Result by Julia")
        jl_x = jl_A \ jl_b
        println(jl_x)
        @test x[1] ≈ jl_x[1]
        @test x[2] ≈ jl_x[2]
        @test x[3] ≈ jl_x[3]
    end
end

for T in [Float64, Float32]
    @testset "template solver part3 $T" begin
        Random.seed!(12345)
        coo = COO{T}("matrixfiles/sample.mtx")
        A = Dense{T}(coo)
        jl_x = one(T) .+ rand(T, get_row(A))
        jl_b = one(T) .+ rand(T, get_row(A))
        x = vector{T}(StdVector(jl_x))
        b = vector{T}(StdVector(jl_b))
        solver = Gomalish.LU{Dense{T},T}()
        @test Gomalish.name(solver) == "monolish::equation::LU"
        precond = Gomalish.none{Dense{T},T}()
        set_create_precond(solver, precond)
        set_apply_precond(solver, precond)
        set_tol(solver, 1.0e-12)
        set_maxiter(solver, get_row(A))
        if solver_check(solve(solver, A, x, b))
            return
        end
        print_all(x)
        jl_A = T[
            2 -1 0
           -1 2 -1
            0 -1 2
        ]
        println("Result by Gomalish/monolish")
        print_all(x)
        println("Result by Julia")
        jl_x = jl_A \ jl_b
        println(jl_x)
        @test x[1] ≈ jl_x[1]
        @test x[2] ≈ jl_x[2]
        @test x[3] ≈ jl_x[3]
    end
end

max_stars_repo_name: AtelierArith/gomalish

max_issues_repo_name: AtelierArith/gomalish

max_forks_repo_name: AtelierArith/gomalish


text:
module TestTransform
    import ManifoldLearning: fit, transform, reconstruct, DataTransform,
                             ZScoreTransform, UnitRangeTransform
    using Base.Test

    X = rand(5, 8)

    t = fit(ZScoreTransform, X; center=false, scale=false)
    Y = transform(t, X)
    @test isa(t, DataTransform)
    @test isempty(t.mean)
    @test isempty(t.scale)
    @test isequal(X, Y)
    @test_approx_eq transform(t, X[:,1]) Y[:,1]
    @test_approx_eq reconstruct(t, Y[:,1]) X[:,1]
    @test_approx_eq reconstruct(t, Y) X

    t = fit(ZScoreTransform, X; center=false)
    Y = transform(t, X)
    @test isa(t, DataTransform)
    @test isempty(t.mean)
    @test length(t.scale) == 5
    @test_approx_eq Y X ./ std(X, 2)
    @test_approx_eq transform(t, X[:,1]) Y[:,1]
    @test_approx_eq reconstruct(t, Y[:,1]) X[:,1]
    @test_approx_eq reconstruct(t, Y) X

    t = fit(ZScoreTransform, X; scale=false)
    Y = transform(t, X)
    @test isa(t, DataTransform)
    @test length(t.mean) == 5
    @test isempty(t.scale)
    @test_approx_eq Y X .- mean(X, 2)
    @test_approx_eq transform(t, X[:,1]) Y[:,1]
    @test_approx_eq reconstruct(t, Y[:,1]) X[:,1]
    @test_approx_eq reconstruct(t, Y) X

    t = fit(ZScoreTransform, X)
    Y = transform(t, X)
    @test isa(t, DataTransform)
    @test length(t.mean) == 5
    @test length(t.scale) == 5
    @test_approx_eq Y (X .- mean(X, 2)) ./ std(X, 2)
    @test_approx_eq transform(t, X[:,1]) Y[:,1]
    @test_approx_eq reconstruct(t, Y[:,1]) X[:,1]
    @test_approx_eq reconstruct(t, Y) X

    t = fit(UnitRangeTransform, X)
    Y = transform(t, X)
    @test length(t.min) == 5
    @test length(t.scale) == 5
    @test_approx_eq Y (X .- minimum(X, 2)) ./ (maximum(X, 2) .- minimum(X, 2))
    @test_approx_eq transform(t, X[:,1]) Y[:,1]
    @test_approx_eq reconstruct(t, Y[:,1]) X[:,1]
    @test_approx_eq reconstruct(t, Y) X

    t = fit(UnitRangeTransform, X; unit=false)
    Y = transform(t, X)
    @test length(t.min) == 5
    @test length(t.scale) == 5
    @test_approx_eq Y X ./ (maximum(X, 2) .- minimum(X, 2))
    @test_approx_eq transform(t, X[:,1]) Y[:,1]
    @test_approx_eq reconstruct(t, Y[:,1]) X[:,1]
    @test_approx_eq reconstruct(t, Y) X
end


text:
fn = "data/vinnakota_kemp_kushmeric_2006_exp45.cellml"
doc = readxml(fn)
docroot = doc.root
maths = extract_mathml(fn)
@test length(maths) == 24


text:
module SecureSessions

using HttpServer
using MbedTLS

export
# Secure cookies
    create_secure_session_cookie,
    get_session_cookie_data,
    invalidate_cookie!,
# Password hashing
    StoredPassword,
    password_is_valid,
    password_is_valid,
# utils
    username_is_permissible,
    password_is_permissible

include("utils.jl")
include("secure_cookies.jl")
include("password_hash.jl")
include("pbkdf2.jl")

# Globals
session_timeout = 5 * 60 * 1000        # Duration of a session's validity in milliseconds
key_length      = 32                   # Key length for AES 256-bit cipher in CBC mode
block_size      = 16                   # IV  length for AES 256-bit cipher in CBC mode
const_key       = csrng(key_length)    # Symmetric key for encrypting secret_keys (with 256-bit encryption)
const_iv        = csrng(block_size)    # IV for encrypting secret_keys
http_only       = true
encrypted_sessions_only = true
timeout_str     = string(convert(Int64, 0.001 * session_timeout))    # Session timeout in seconds, represented as a string

end # module


text:
using Base.Test
using DataFrames

let
	p = PooledDataArray(DataVector[9,9,8,NA,1,1])
	pcopy = copy(p)
	@assert levels(p) == [1,8,9]
	@assert levels(set_levels(p, ["a", "b", "c"])) == ["a", "b", "c"]
	@assert removeNA(set_levels(p, DataVector["a", "b", NA])) == ["b", "a", "a"]
	@assert removeNA(set_levels(p, DataVector["a", "b", "a"])) == ["a", "a", "b", "a", "a"]
	@assert levels(set_levels(p, DataVector["a", "b", "a"])) == ["a", "b"]
	@assert levels(set_levels(p, [1 => 111])) == [111, 8, 9]
	@assert levels(set_levels(p, [1 => 111, 8 => NA])) == [111, 9]
	@assert levels(PooledDataArray(p, [9,8,1])) == [9,8,1]
	@assert levels(PooledDataArray(p, [9,8])) == [9,8]
	@assert removeNA(PooledDataArray(p, [9,8])) == [9,9,8]
	@assert levels(PooledDataArray(p, levels(p)[[3,2,1]])) == [9,8,1]
	v = [1:6]
	@assert isequal(p, reorder(p))
	@assert levels(reorder(p, v)) == [9,8,1]
	df = @DataFrame(v => v, x => rand(6))
	@assert levels(reorder(minimum, p, df)) == [9,8,1]
	@assert isequal(p, pcopy)

	@assert levels(set_levels!(copy(p), [10,80,90])) == [10, 80, 90]
	@assert levels(set_levels!(copy(p), [1,8,1])) == [1, 8]
	@assert levels(set_levels!(copy(p), DataVector[1,8,NA])) == [1, 8]
	@assert levels(set_levels!(copy(p), [1,8,9, 10])) == [1, 8, 9, 10]
	@assert levels(set_levels!(copy(p), [1 => 111])) == [111, 8, 9]
	@assert levels(set_levels!(copy(p), [1 => 111, 8 => NA])) == [111, 9]

	pp = PooledDataArray(Any[])
	@assert length(pp) == 0
	@assert length(levels(pp)) == 0

	# Test explicitly setting refs type
	testarray = [1,1,2,2,0,0,3,3]
	testdata = DataVector[1,1,2,2,0,0,3,3]
	for t in Any[testarray, testdata]
	    for R in [Uint8, Uint16, Uint32, Uint64]
	        @assert eltype(PooledDataArray(t, R).refs) == R
	        @assert eltype(PooledDataArray(t, [1,2,3], R).refs) == R
	        @assert eltype(PooledDataArray(t, [1,2,3], t .== 0, R).refs) == R
	    end
	end
end


text:
using ITensors

#
# DMRG calculation of the extended Hubbard model
# ground state wavefunction, and spin densities
#

# Usage:
#
# Run:
#
# $ julia exthubbard.jl input.jl
#

# Include the specified input file
# Include the file from the first positional
# argument. If none is specified, use
# input.jl by default.
include(get(ARGS, 1, "input.jl"))

sweeps = Sweeps(nsweep, sweeps_args)

@show sweeps
@show N, Npart

sites = siteinds("Electron", N; conserve_qns=true)

ampo = AutoMPO()
for i in 1:N
  ampo .+= U, "Nupdn", i
end
for b in 1:(N - 1)
  ampo .+= -t1, "Cdagup", b, "Cup", b + 1
  ampo .+= -t1, "Cdagup", b + 1, "Cup", b
  ampo .+= -t1, "Cdagdn", b, "Cdn", b + 1
  ampo .+= -t1, "Cdagdn", b + 1, "Cdn", b
  ampo .+= V1, "Ntot", b, "Ntot", b + 1
end
for b in 1:(N - 2)
  ampo .+= -t2, "Cdagup", b, "Cup", b + 2
  ampo .+= -t2, "Cdagup", b + 2, "Cup", b
  ampo .+= -t2, "Cdagdn", b, "Cdn", b + 2
  ampo .+= -t2, "Cdagdn", b + 2, "Cdn", b
end
H = MPO(ampo, sites)

state = ["Emp" for n in 1:N]
p = Ref(Npart)
for i in N:-1:1
  if p[] > i
    println("Doubly occupying site $i")
    state[i] = "UpDn"
    p[] -= 2
  elseif p[] > 0
    println("Singly occupying site $i")
    state[i] = (isodd(i) ? "Up" : "Dn")
    p[] -= 1
  end
end
# Initialize wavefunction to be bond
# dimension 10 random MPS with number
# of particles the same as `state`
psi0 = randomMPS(sites, state, 10)

# Check total number of particles:
@show flux(psi0)

# Start DMRG calculation:
energy, psi = dmrg(H, psi0, sweeps)

upd = fill(0.0, N)
dnd = fill(0.0, N)
for j in 1:N
  orthogonalize!(psi, j)
  psidag_j = dag(prime(psi[j], "Site"))
  upd[j] = scalar(psidag_j * op(sites, "Nup", j) * psi[j])
  dnd[j] = scalar(psidag_j * op(sites, "Ndn", j) * psi[j])
end

println("Up Density:")
for j in 1:N
  println("$j $(upd[j])")
end
println()

println("Dn Density:")
for j in 1:N
  println("$j $(dnd[j])")
end
println()

println("Total Density:")
for j in 1:N
  println("$j $(upd[j]+dnd[j])")
end
println()

println("\nGround State Energy = $energy")


text:
using GlobalSensitivity, Test

f = X -> X[1] + 2 * X[2] + 3 * X[3] + 4 * X[7] * X[12]
res1 = gsa(f,FractionalFactorial(),num_params = 12,N=10)

@test res1[1][1:3] == [1.0, 2.0, 3.0]
@test all(res1[1][4:end] .== 0.0)


text:
module Test

using ..BasisFunctions, DomainSets, DomainIntegrals, StaticArrays
using FillArrays: Ones, Eye
using Test, LinearAlgebra
BF = BasisFunctions

using BasisFunctions: dimension_tuple

export Delimit, delimit
export affine_point_in_domain, point_outside_domain, random_point_in_domain, fixed_point_in_domain
export random_index
export test_tolerance
include("test_util_functions.jl")

export test_generic_dict_interface
export supports_approximation, supports_interpolation
export suitable_interpolation_grid, suitable_function
include("test_dictionary.jl")

export test_generic_operator_interface
include("test_operator.jl")

export test_orthogonality_orthonormality
include("test_functionality.jl")

export generic_test_discrete_measure, generic_test_measure
function generic_test_discrete_measure(measure)
    io = IOBuffer()
    show(io, measure)
    @test length(take!(io))>0
    @test length(weights(measure)) == length(points(measure))

end

function generic_test_measure(measure)
    io = IOBuffer()
    show(io, measure)
    @test length(take!(io))>0
    support(measure)
    x = rand()
    @test weightfun(measure,x) ≈ weightfunction(measure)(x)
end
end


text:
struct SaturnBarycenter <: abstractBarycenter end

naif_id(::Type{SaturnBarycenter}) = 6
parent(::Type{SaturnBarycenter}) = SSB


text:
# This file was generated, do not modify it. # hide
ẑ₁ = (W, z) |> KNNRegressor(K=5)
ẑ₂ = (W, z) |> RidgeRegressor(lambda=2.5);

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
"""
    Lib

Low-level blosc api - wrappers under corresponding c functions
"""
module Lib
    using Blosc2_jll
    const lib = libblosc2
    include("defines.jl")
    include("blosc1/blosc1.jl")
    include("blosc2/blosc2.jl")
end


text:
using CSV
using DataFrames
using DelimitedFiles
using SLiDE

# QUESTION: Do we have mapping capabilities for all units?
DATA_DIR = joinpath("data", "input")
MAP_DIR = joinpath("data", "coremaps", "parse")

# Read file names from output directory and isolate csv output files.
files = readdir(joinpath(SLiDE.SLIDE_DIR, DATA_DIR))
files = files[occursin.(".csv", files)]

# Read files and save all units present as well as the associated file name.
df_units = DataFrame()
for f in files
    df = read_file(joinpath(DATA_DIR,f))
    df = unique(df[:,:units])
    d = Dict(:file => f, :units => [[col for col in eachcol(df)]...;])
    global df_units = [df_units; DataFrame(d)]
end

# Read units map for comparison.
df_units_map = read_file(joinpath(MAP_DIR, "units.csv"))

# Units in the output DataFrame that are not in the map::
unlisted_units = setdiff(df_units[:,:units], df_units_map[:,:to])
df_unlisted = filter(row -> row[:units] in unlisted_units, df_units)

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE