text:
"""
    Generator

Abstract type representing a discretised infinitesimal generator of a FLuidQueue.
Behaves much like a square matrix.
"""
abstract type Generator{T<:Mesh} <: AbstractMatrix{Float64} end

const UnionVectors = Union{StaticArrays.SVector,Vector{Float64}}
const UnionArrays = Union{Array{Float64,2},StaticArrays.SMatrix}
const BoundaryFluxTupleType = Union{
    NamedTuple{(:upper,:lower), Tuple{
        NamedTuple{(:in,:out),Tuple{Vector{Float64},Vector{Float64}}},
        NamedTuple{(:in,:out),Tuple{Vector{Float64},Vector{Float64}}}
        }
    },
    NamedTuple{(:upper, :lower), Tuple{
        NamedTuple{(:in, :out), Tuple{StaticArrays.SVector, StaticArrays.SVector}},
        NamedTuple{(:in, :out), Tuple{StaticArrays.SVector, StaticArrays.SVector}}
        }
    }
}

struct OneBoundaryFlux{T<:Union{StaticArrays.SVector,Vector{Float64}}}
    in::T
    out::T
end
struct BoundaryFlux{T<:Union{StaticArrays.SVector,Vector{Float64}}}
    upper::OneBoundaryFlux{T}
    lower::OneBoundaryFlux{T}
end

"""
    LazyGenerator{<:Mesh} <: Generator{<:Mesh}

A lazy representation of a block matrix with is a generator of a DiscretisedFluidQueue.

Lower memory requirements than FullGenerator but aritmetic operations and indexing may be slower.

# Arguments:
- `dq::DiscretisedFluidQueue`:
- `blocks::Tuple{Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2}}`:
    Block matrices describing the flow of mass within and between cells. `blocks[1]` is the lower
    diagonal block describing the flow of mass from cell k+1 to cell (k for phases
    with negative rate only). `blocks[2] (blocks[3])` is the
    diagonal block describing the flow of mass within a cell for a phase with positive (negative) rate.
    `blocks[4]` is the upper diagonal block describing the flow of mass from cell k to k+1 (for phases
    with positive rate only).
- `boundary_flux::BoundaryFlux`: A named tuple structure such that
        - `boundary_flux.lower.in`: describes flow of density into lower boundary
        - `boundary_flux.lower.out`: describes flow of density out of lower boundary
        - `boundary_flux.upper.in`: describes flow of density into upper boundary
        - `boundary_flux.upper.out`: describes flow of density out of upper boundary
- `D::Union{Array{Float64, 2}, LinearAlgebra.Diagonal{Bool, Array{Bool, 1}}}`: An array describing
    how the flow of density changes when the phase process jumps between phases with different memberships.
    This is the identity for FV and DG schemes.
"""
struct LazyGenerator{T}  <: Generator{T}
    dq::DiscretisedFluidQueue{T}
    blocks::NTuple{4,AbstractMatrix{Float64}}
    boundary_flux::BoundaryFlux
    D::AbstractMatrix{Float64}
    function LazyGenerator(
        dq::DiscretisedFluidQueue{T},
        blocks::NTuple{4,AbstractMatrix{Float64}},
        boundary_flux::BoundaryFlux,
        D::AbstractMatrix{Float64},
    ) where T
        s = size(blocks[1])
        for b in 1:4
            checksquare(blocks[b])
            !(s == size(blocks[b])) && throw(DomainError("blocks must be the same size"))
        end
        checksquare(D)
        !(s == size(D)) && throw(DomainError("blocks must be the same size as D"))

        return new{T}(dq,blocks,boundary_flux,D)
    end
end
function LazyGenerator(
    dq::DiscretisedFluidQueue,
    blocks::NTuple{3,AbstractMatrix{Float64}},
    boundary_flux::OneBoundaryFlux,
    D::AbstractMatrix{Float64},
)
    blocks = (blocks[1],blocks[2],blocks[2],blocks[3])
    boundary_flux = BoundaryFlux(boundary_flux, boundary_flux)
    return LazyGenerator(dq,blocks,boundary_flux,D)
end

function static_generator(lz)
    sz = size(lz.blocks[1], 1)
    ex_smatrix = StaticArrays.SMatrix{sz, sz, Float64}
    ex_svector = StaticArrays.SVector{sz, Float64}
    b1 = ex_smatrix(lz.blocks[1])
    b2 = ex_smatrix(lz.blocks[2])
    b3 = ex_smatrix(lz.blocks[3])
    b4 = ex_smatrix(lz.blocks[4])
    uprin = ex_svector(lz.boundary_flux.upper.in)
    uprout = ex_svector(lz.boundary_flux.upper.out)
    lwrin = ex_svector(lz.boundary_flux.lower.in)
    lwrout = ex_svector(lz.boundary_flux.lower.out)
    D = ex_smatrix(lz.D)
    dq = lz.dq

    # return
    out = LazyGenerator(dq,(b1,b2,b3,b4),
        BoundaryFlux(OneBoundaryFlux(uprin,uprout),OneBoundaryFlux(lwrin,lwrout)),
        D)
    return out
end

"""
    build_lazy_generator(dq::DiscretisedFluidQueue; v::Bool = false)

Build a lazy representation of the generator of a discretised fluid queue.
"""
function build_lazy_generator(dq::DiscretisedFluidQueue; v::Bool=false)
    throw(DomainError("Can construct LazyGenerator for DGMesh, FRAPMesh, only"))
end

"""
    size(B::LazyGenerator)
"""
function size(B::LazyGenerator)
    sz = total_n_bases(B.dq) + N₋(B.dq) + N₊(B.dq)
    return (sz,sz)
end
size(B::LazyGenerator, n::Int) =
    (n∈[1,2]) ? size(B)[n] : throw(DomainError("Lazy generator is a matrix, index must be 1 or 2"))

_check_phase_index(i::Int,model::Model) =
    (i∉phases(model)) && throw(DomainError("i is not a valid phase in model"))
_check_mesh_index(k::Int,mesh::Mesh) =
    !(1<=k<=n_intervals(mesh)) && throw(DomainError("k in not a valid cell"))
_check_basis_index(p::Int,mesh::Mesh) = !(1<=p<=n_bases_per_cell(mesh))

function _map_to_index_interior((k,i,p)::NTuple{3,Int},dq::DiscretisedFluidQueue)
    # i phase
    # k cell
    # p basis

    _check_phase_index(i,dq.model)
    _check_mesh_index(k,dq.mesh)
    _check_basis_index(p,dq.mesh)

    idx = (i-1)*n_bases_per_cell(dq) + (k-1)*n_bases_per_level(dq) + p
    return N₋(dq) + idx
end
function _map_to_index_boundary((k,i,p)::NTuple{3,Int},dq::DiscretisedFluidQueue)
    # i phase
    # k cell
    # p basis
    _check_phase_index(i,dq.model)
    !(p==1)&&throw(DomainError("p must be 1 at the boundary"))
    !((k==0)||(k==n_intervals(dq)+1))&&throw(DomainError("k not a boundary index of dq"))
    if (k==0)&&(_has_left_boundary(dq.model.S,i))
        idx = N₋(dq.model.S[1:i])
    elseif (k==n_intervals(dq)+1)&&_has_right_boundary(dq.model.S,i)
        idx = N₊(dq.model.S[1:i]) + total_n_bases(dq) + N₋(dq)
    else
        throw(DomainError(string((k,i,p))*" is not a valid index of dq"))
    end
    return idx
end
function _map_to_index((k,i,p)::NTuple{3,Int},dq::DiscretisedFluidQueue)
    if (k==0)||(k==n_intervals(dq)+1)
        idx = _map_to_index_boundary((k,i,p),dq)
    else
        idx = _map_to_index_interior((k,i,p),dq)
    end
    return idx
end
_getindex_from_tuple(B,(ks,is,ps)::Tuple,(ls,js,qs)::Tuple) = _getindex_from_tuple(B,B.dq,(ks,is,ps)::Tuple,(ls,js,qs)::Tuple)
function _build_index_from_tuple(dq,(ks,is,ps)::Tuple)
    (typeof(is)==Colon)&&(is=1:n_phases(dq))
    (typeof(ks)==Colon)&&(ks=0:n_intervals(dq)+1)
    (typeof(ps)==Colon)&&(ps=1:n_bases_per_cell(dq))

    idx = Int[]
    for k in ks
        for i in is
            for p in ps
                if (k==0)&&(p==1)&&_has_left_boundary(dq.model.S,i)
                    push!(idx,_map_to_index((0,i,1),dq))
                elseif (k==n_intervals(dq)+1)&&(p==1)&&_has_right_boundary(dq.model.S,i)
                    push!(idx,_map_to_index((n_intervals(dq)+1,i,1),dq))
                elseif (k!=0)&&(k!=n_intervals(dq)+1)
                    push!(idx,_map_to_index((k,i,p),dq))
                end
            end
        end
    end
    return idx
end
function _getindex_from_tuple(B,dq,(ks,is,ps)::Tuple,(ls,js,qs)::Tuple)
    # (typeof(js)==Colon)&&(js=1:n_phases(dq))
    # (typeof(ls)==Colon) ? (ls=0:n_intervals(dq)+1; ls_colon=true) : (ls_colon=false)
    # (typeof(qs)==Colon)&&(qs=1:n_bases_per_cell(dq))

    rows = _build_index_from_tuple(dq,(ks,is,ps))
    cols = _build_index_from_tuple(dq,(ls,js,qs))

    # cols = Int[]
    # for l in ls
    #     for j in js
    #         for q in qs
    #             if (l==0)&&(q==1)&&_has_left_boundary(dq.model.S,j)
    #                 push!(cols,_map_to_index((0,j,1),dq))
    #             elseif (l==n_intervals(dq)+1)&&(q==1)&&_has_right_boundary(dq.model.S,j)
    #                 push!(cols,_map_to_index((n_intervals(dq)+1,j,1),dq))
    #             elseif (l!=0)&&(l!=n_intervals(dq)+1)
    #                 push!(cols,_map_to_index((l,j,q),dq))
    #             end
    #         end
    #     end
    # end
    return B[rows,cols]
end

function Base.getindex(B::LazyGenerator,(ks,is,ps)::Tuple,(ls,js,qs)::Tuple)
    return _getindex_from_tuple(B,(ks,is,ps),(ls,js,qs))
end

_is_boundary_index(n::Int,B::LazyGenerator) = (n∈1:N₋(B.dq))||(n∈(size(B,1).-(0:N₊(B.dq)-1)))
function _map_from_index_interior(n::Int,B::LazyGenerator)
    # n matrix index to map to phase, cell, basis
    (!(1<=n<=size(B,1))||_is_boundary_index(n,B))&&throw(DomainError(n,"not a valid interior index"))

    n -= (N₋(B.dq)+1)
    P = n_bases_per_cell(B.dq)

    k = n÷n_bases_per_level(B.dq) + 1
    i = mod(n,n_bases_per_level(B.dq))÷n_bases_per_cell(B.dq) + 1
    p = mod(n,P) + 1
    return k, i, p
end
function _map_from_index_boundary(n::Int,B::LazyGenerator)
    # n matrix index to map to phase at boundary
    (!_is_boundary_index(n,B))&&throw(DomainError("not a valid boundary index"))

    if n>N₋(B.dq)
        i₊ = n-N₋(B.dq)-total_n_bases(B.dq)
        i = 1
        for j in phases(B.dq)
            (i₊==N₊(B.dq.model.S[1:j])) && break
            i += 1
        end
    else
        i₋ = n
        i = 1
        for j in phases(B.dq)
            (i₋==N₋(B.dq.model.S[1:j])) && break
            i += 1
        end
    end

    return i
end

function fast_mul(u::AbstractMatrix{Float64}, 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)"))

    if output_type <: SparseArrays.SparseMatrixCSC
        v = SparseArrays.spzeros(sz_u_1,sz_B_2)
    else
        v = zeros(sz_u_1,sz_B_2)
    end

    Kp = n_bases_per_cell(B.dq) # K = n_intervals(mesh), p = n_bases_per_cell(mesh)
    C = rates(B.dq)
    n₋ = N₋(B.dq)
    n₊ = N₊(B.dq)

    # boundaries
    # at lower
    _lmul_at_lwr_bndry!(v,u,B.dq.model.T,B.dq.model.S,n₋)
    # in to lower
    _lmul_into_lwr_bndry!(v,u,B.dq.model,B.dq.mesh,Kp,n₋,C,B.boundary_flux.lower.in)

    # out of lower
    _lmul_out_lwr_bndry!(v,u,B.dq.mesh,B.dq.model,Kp,n₋,B.boundary_flux.lower.in,B.boundary_flux.lower.out)

    # at upper
    _lmul_at_upr_bndry!(v,u,B.dq.model.T,B.dq.model.S,n₊)
    # in to upper
    _lmul_into_upr_bndry!(v,u,B.dq.model,B.dq.mesh,Kp,n₋,n₊,C,B.boundary_flux.upper.in)

    # out of upper
    _lmul_out_upr_bndry!(v,u,B.dq.mesh,B.dq.model,Kp,n₋,n₊,B.boundary_flux.upper.in,B.boundary_flux.upper.out)

    # innards
    _lmul_innards!(v,u,B.dq.model,B.dq.mesh,B.blocks,B.D,C,n₋)
    return v
end
function _lmul_innards_generic!(v,u,model,mesh,tmp_blocks,D,C,n₋,k)
    for j in 1:n_phases(model)
        i_idx = (k-1)*n_phases(model)*n_bases_per_cell(mesh) + n₋ .+ (1:n_bases_per_cell(mesh))
        j_idx = i_idx .+ (j-1)*n_bases_per_cell(mesh)
        (k>1) && ((C[j]>0.0) && (v[:,j_idx] += C[j]*u[:,j_idx.-n_bases_per_cell(mesh)*n_phases(model)]*tmp_blocks[4]))
        (k<n_intervals(mesh)) && ((C[j]<0.0) && (v[:,j_idx] += abs(C[j])*u[:,j_idx.+n_bases_per_cell(mesh)*n_phases(model)]*tmp_blocks[1]))
        for i in 1:n_phases(model)
            if (typeof(mesh)<:FRAPMesh)&&(membership(model.S,i)!=membership(model.S,j)) # + to - change
                v[:,j_idx] += u[:,i_idx]*D*model.T[i,j]
            elseif (i!=j)||(C[i]==0.0)
                v[:,j_idx] += u[:,i_idx]*model.T[i,j]
            elseif C[i]>0.0 # i==j
                v[:,j_idx] += u[:,i_idx]*C[i]*tmp_blocks[2] + u[:,i_idx]*model.T[i,j]
            elseif C[i]<0.0 # i==j
                v[:,j_idx] += u[:,i_idx]*abs(C[i])*tmp_blocks[3] + u[:,i_idx]*model.T[i,j]
            end
            i_idx = i_idx .+ n_bases_per_cell(mesh)
        end
    end
end
function _lmul_innards!(v,u,model,mesh::Mesh{Vector{Float64}},blocks,D,C,n₋)
    for k in 1:n_intervals(mesh)
        if 1<k<n_intervals(mesh)
            tmp_blocks = (blocks[1]/Δ(mesh,k+1), blocks[2]/Δ(mesh,k), blocks[3]/Δ(mesh,k), blocks[4]/Δ(mesh,k-1))
        elseif k==1
            tmp_blocks = (blocks[1]/Δ(mesh,k+1), blocks[2]/Δ(mesh,k), blocks[3]/Δ(mesh,k), zeros(0,0))
        else
            tmp_blocks = (zeros(0,0), blocks[2]/Δ(mesh,k), blocks[3]/Δ(mesh,k), blocks[4]/Δ(mesh,k-1))
        end
        _lmul_innards_generic!(v,u,model,mesh,tmp_blocks,D,C,n₋,k)
    end
end
function _lmul_innards!(v,u,model,mesh::Mesh{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}},blocks,D,C,n₋)
    tmp_blocks = (blocks[1]/Δ(mesh,1), blocks[2]/Δ(mesh,1), blocks[3]/Δ(mesh,1), blocks[4]/Δ(mesh,1))
    for k in 1:n_intervals(mesh)
        _lmul_innards_generic!(v,u,model,mesh,tmp_blocks,D,C,n₋,k)
    end
end

function _lmul_at_lwr_bndry!(v,u,T,S,n₋)
    v[:,1:n₋]+=u[:,1:n₋]*T[_has_left_boundary.(S),_has_left_boundary.(S)]
    return nothing
end
function _lmul_at_upr_bndry!(v,u,T,S,n₊)
    v[:,end-n₊+1:end]+=u[:,end-n₊+1:end]*T[_has_right_boundary.(S),_has_right_boundary.(S)]
    return nothing
end
function _lmul_into_lwr_bndry!(v,u,model::BoundedFluidQueue,mesh,Kp,n₋,C,bndry_flux_in_lwr)
    idxdown = n₋ .+ ((1:n_bases_per_cell(mesh)).+Kp*(findall(negative_phases(model)) .- 1)')[:]
    v[:,1:n₋] += u[:,idxdown]*LinearAlgebra.kron(
        # LinearAlgebra.diagm(0 => abs.(C[_has_left_boundary.(model.S)])),
        abs.(C[negative_phases(model)]).*model.P_lwr[:,_has_left_boundary.(model.S)],
        bndry_flux_in_lwr/Δ(mesh,1),
    )
    return nothing
end
function _lmul_into_upr_bndry!(v,u,model::BoundedFluidQueue,mesh,Kp,n₋,n₊,C,bndry_flux_in_upr)
    idxup = n₋ .+ ((1:n_bases_per_cell(mesh)) .+ Kp*(findall(positive_phases(model)) .- 1)')[:] .+
        (n_bases_per_phase(mesh)*n_phases(model).-n_bases_per_cell(mesh)*n_phases(model))
    v[:,end-n₊+1:end] += u[:,idxup]*LinearAlgebra.kron(
        # LinearAlgebra.diagm(0 => C[_has_right_boundary.(model.S)]),
        C[positive_phases(model)].*model.P_upr[:,_has_right_boundary.(model.S)],
        bndry_flux_in_upr/Δ(mesh,n_intervals(mesh)),
    )
    return nothing
end
function _lmul_out_lwr_bndry_generic!(v,u,mesh,model,Kp,n₋,bndry_flux_out_lwr)
    idxup = n₋ .+ (Kp*(findall(positive_phases(model)).-1)' .+ (1:n_bases_per_cell(mesh)))[:]
    v[:,idxup] += u[:,1:n₋]*LinearAlgebra.kron(model.T[_has_left_boundary.(model.S),positive_phases(model)],bndry_flux_out_lwr')
    return idxup
end
function _lmul_out_lwr_bndry!(v,u,mesh,model::BoundedFluidQueue,Kp,n₋,bndry_flux_in_lwr,bndry_flux_out_lwr)
    idxup = _lmul_out_lwr_bndry_generic!(v,u,mesh,model,Kp,n₋,bndry_flux_out_lwr)
    idxdown = n₋ .+ ((1:n_bases_per_cell(mesh)).+Kp*(findall(negative_phases(model)) .- 1)')[:]
    v[:,idxup] += u[:,idxdown]*LinearAlgebra.kron(
        # LinearAlgebra.diagm(0 => abs.(C[_has_left_boundary.(model.S)])),
        abs.(rates(model)[negative_phases(model)]).*model.P_lwr[:,positive_phases(model)],
        bndry_flux_in_lwr*bndry_flux_out_lwr'/Δ(mesh,1),
    )
    return nothing
end
function _lmul_out_upr_bndry_generic!(v,u,mesh,model,Kp,n₋,n₊,bndry_flux_out_upr)
    idxdown = n₋ .+ (Kp*(findall(negative_phases(model)).-1)' .+ (1:n_bases_per_cell(mesh)))[:] .+
        (n_bases_per_phase(mesh)*n_phases(model).-n_bases_per_cell(mesh)*n_phases(model))
    v[:,idxdown] += u[:,end-n₊+1:end]*LinearAlgebra.kron(model.T[_has_right_boundary.(model.S),negative_phases(model)],bndry_flux_out_upr')
    return idxdown
end
function _lmul_out_upr_bndry!(v,u,mesh,model::BoundedFluidQueue,Kp,n₋,n₊,bndry_flux_in_upr,bndry_flux_out_upr)
    idxdown = _lmul_out_upr_bndry_generic!(v,u,mesh,model,Kp,n₋,n₊,bndry_flux_out_upr)
    idxup = n₋ .+ ((1:n_bases_per_cell(mesh)) .+ Kp*(findall(positive_phases(model)) .- 1)')[:] .+
        (n_bases_per_phase(mesh)*n_phases(model).-n_bases_per_cell(mesh)*n_phases(model))
    v[:,idxdown] += u[:,idxup]*LinearAlgebra.kron(
        # LinearAlgebra.diagm(0 => C[_has_right_boundary.(model.S)]),
        abs.(rates(model)[positive_phases(model)]).*model.P_upr[:,negative_phases(model)],
        bndry_flux_in_upr*bndry_flux_out_upr'/Δ(mesh,n_intervals(mesh)),
    )
    return nothing
end

function fast_mul(B::LazyGenerator,u::AbstractMatrix{Float64})
    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_B_2 == sz_u_1) && throw(DomainError("Dimension mismatch, u*B, length(u) must be size(B,1)"))

    if output_type <: SparseArrays.SparseMatrixCSC
        v = SparseArrays.spzeros(sz_B_1,sz_u_2)
    else
        v = zeros(sz_B_1,sz_u_2)
    end
    model = B.dq.model
    mesh = B.dq.mesh
    Kp = n_bases_per_cell(B.dq) # K = n_intervals(mesh), p = n_bases_per_cell(mesh)
    C = rates(B.dq)
    n₋ = N₋(B.dq)
    n₊ = N₊(B.dq)

    # boundaries
    # at lower
    _rmul_at_lwr_bndry!(v,u,model.T,model.S,n₋)
    # in to lower
    _rmul_into_lwr_bndry!(v,u,model,mesh,Kp,n₋,C,B.boundary_flux.lower.in)

    # out of lower
    _rmul_out_lwr_bndry!(v,u,mesh,model,Kp,n₋,B.boundary_flux.lower.in,B.boundary_flux.lower.out)

    # at upper
    _rmul_at_upr_bndry!(v,u,model.T,model.S,n₊) # v[:,end-n₊+1:end] += u[:,end-n₊+1:end]*model.T[_has_right_boundary.(model.S),_has_right_boundary.(model.S)]
    # in to upper
    _rmul_into_upr_bndry!(v,u,model,mesh,Kp,n₋,n₊,C,B.boundary_flux.upper.in)

    # out of upper
    _rmul_out_upr_bndry!(v,u,mesh,model,Kp,n₋,n₊,B.boundary_flux.upper.in,B.boundary_flux.upper.out)

    # innards
    _rmul_innards!(v,u,B.dq.model,B.dq.mesh,B.blocks,B.D,C,n₋)
    return v
end
function _rmul_innards_generic!(v,u,model,mesh,tmp_blocks,D,C,n₋,k)
    for j in 1:n_phases(model)
        i_idx = (k-1)*n_phases(model)*n_bases_per_cell(mesh) + n₋ .+ (1:n_bases_per_cell(mesh))
        j_idx = i_idx .+ (j-1)*n_bases_per_cell(mesh)
        (k>1) && ((C[j]>0.0) && (v[j_idx.-n_bases_per_cell(mesh)*n_phases(model),:] += C[j]*tmp_blocks[4]*u[j_idx,:]))
        (k<n_intervals(mesh)) && ((C[j]<0.0) && (v[j_idx.+n_bases_per_cell(mesh)*n_phases(model),:] += abs(C[j])*tmp_blocks[1]*u[j_idx,:]))
        for i in 1:n_phases(model)
            if (typeof(mesh)<:FRAPMesh)&&(membership(model.S,i)!=membership(model.S,j)) # + to - change
                v[i_idx,:] += D*u[j_idx,:]*model.T[i,j]
            elseif (i!=j)||(C[i]==0.0)
                v[i_idx,:] += model.T[i,j]*u[j_idx,:]
            elseif C[i]>0.0 # i==j
                v[i_idx,:] += C[i]*tmp_blocks[2]*u[j_idx,:] + model.T[i,j]*u[j_idx,:]
            elseif C[i]<0.0 # i==j
                v[i_idx,:] += abs(C[i])*tmp_blocks[3]*u[j_idx,:] + model.T[i,j]*u[j_idx,:]
            end
            i_idx = i_idx .+ n_bases_per_cell(mesh)
        end
    end
end
function _rmul_innards!(v,u,model,mesh::Mesh{Vector{Float64}},blocks,D,C,n₋)
    for k in 1:n_intervals(mesh)
        if 1<k<n_intervals(mesh)
            tmp_blocks = (blocks[1]/Δ(mesh,k+1), blocks[2]/Δ(mesh,k), blocks[3]/Δ(mesh,k), blocks[4]/Δ(mesh,k-1))
        elseif k==1
            tmp_blocks = (blocks[1]/Δ(mesh,k+1), blocks[2]/Δ(mesh,k), blocks[3]/Δ(mesh,k), zeros(0,0))
        else
            tmp_blocks = (zeros(0,0), blocks[2]/Δ(mesh,k), blocks[3]/Δ(mesh,k), blocks[4]/Δ(mesh,k-1))
        end
        _lmul_innards_generic!(v,u,model,mesh,tmp_blocks,D,C,n₋,k)
    end
end
function _rmul_innards!(v,u,model,mesh::Mesh{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}},blocks,D,C,n₋)
    tmp_blocks = (blocks[1]/Δ(mesh,1), blocks[2]/Δ(mesh,1), blocks[3]/Δ(mesh,1), blocks[4]/Δ(mesh,1))
    for k in 1:n_intervals(mesh)
        _rmul_innards_generic!(v,u,model,mesh,tmp_blocks,D,C,n₋,k)
    end
end

function _rmul_at_lwr_bndry!(v,u,T,S,n₋)
    v[1:n₋,:] += T[_has_left_boundary.(S),_has_left_boundary.(S)]*u[1:n₋,:]
    return nothing
end
function _rmul_at_upr_bndry!(v,u,T,S,n₊)
    v[end-n₊+1:end,:] += T[_has_right_boundary.(S),_has_right_boundary.(S)]*u[end-n₊+1:end,:]
    return nothing
end
function _rmul_into_lwr_bndry!(v,u,model::BoundedFluidQueue,mesh,Kp,n₋,C,bndry_flux_in_lwr)
    idxdown = n₋ .+ ((1:n_bases_per_cell(mesh)).+Kp*(findall(negative_phases(model)) .- 1)')[:]
    v[idxdown,:] += LinearAlgebra.kron(
        abs.(C[negative_phases(model)]).*model.P_lwr[:,_has_left_boundary.(model.S)],
        bndry_flux_in_lwr/Δ(mesh,1),
    )*u[1:n₋,:]
    return nothing
end
function _rmul_into_upr_bndry!(v,u,model::BoundedFluidQueue,mesh,Kp,n₋,n₊,C,bndry_flux_in_upr)
    idxup = n₋ .+ ((1:n_bases_per_cell(mesh)).+Kp*(findall(positive_phases(model)) .- 1)')[:] .+
        (n_bases_per_phase(mesh)*n_phases(model).-n_bases_per_cell(mesh)*n_phases(model))
    v[idxup,:] += LinearAlgebra.kron(
        C[positive_phases(model)].*model.P_upr[:,_has_right_boundary.(model.S)],
        bndry_flux_in_upr/Δ(mesh,n_intervals(mesh)),
    )*u[end-n₊+1:end,:]
    return nothing
end
function _rmul_out_lwr_bndry_generic!(v,u,mesh,model,Kp,n₋,bndry_flux_out_lwr)
    idxup = n₋ .+ (Kp*(findall(positive_phases(model)).-1)' .+ (1:n_bases_per_cell(mesh)))[:]
    v[1:n₋,:] += LinearAlgebra.kron(model.T[_has_left_boundary.(model.S),positive_phases(model)],bndry_flux_out_lwr')*u[idxup,:]
    return idxup
end
function _rmul_out_lwr_bndry!(v,u,mesh,model::BoundedFluidQueue,Kp,n₋,bndry_flux_in_lwr,bndry_flux_out_lwr)
    idxup = _rmul_out_lwr_bndry_generic!(v,u,mesh,model,Kp,n₋,bndry_flux_out_lwr)
    idxdown = n₋ .+ ((1:n_bases_per_cell(mesh)).+Kp*(findall(negative_phases(model)) .- 1)')[:]
    v[idxdown,:] += LinearAlgebra.kron(
        # LinearAlgebra.diagm(0 => abs.(C[_has_left_boundary.(model.S)])),
        abs.(rates(model)[negative_phases(model)]).*model.P_lwr[:,positive_phases(model)],
        bndry_flux_in_lwr*bndry_flux_out_lwr'/Δ(mesh,1),
    )*u[idxup,:]
    return nothing
end
function _rmul_out_upr_bndry_generic!(v,u,mesh,model,Kp,n₋,n₊,bndry_flux_out_upr)
    idxdown = n₋ .+ (Kp*(findall(negative_phases(model)).-1)' .+ (1:n_bases_per_cell(mesh)))[:] .+
        (n_bases_per_phase(mesh)*n_phases(model).-n_bases_per_cell(mesh)*n_phases(model))
    v[end-n₊+1:end,:] += LinearAlgebra.kron(model.T[_has_right_boundary.(model.S),negative_phases(model)],bndry_flux_out_upr')*u[idxdown,:]
    return idxdown
end
function _rmul_out_upr_bndry!(v,u,mesh,model::BoundedFluidQueue,Kp,n₋,n₊,bndry_flux_in_upr,bndry_flux_out_upr)
    idxdown = _rmul_out_upr_bndry_generic!(v,u,mesh,model,Kp,n₋,n₊,bndry_flux_out_upr)
    idxup = n₋ .+ ((1:n_bases_per_cell(mesh)) .+ Kp*(findall(positive_phases(model)) .- 1)')[:] .+
        (n_bases_per_phase(mesh)*n_phases(model).-n_bases_per_cell(mesh)*n_phases(model))
    v[idxup,:] += LinearAlgebra.kron(
        # LinearAlgebra.diagm(0 => C[_has_right_boundary.(model.S)]),
        abs.(rates(model)[positive_phases(model)]).*model.P_upr[:,negative_phases(model)],
        bndry_flux_in_upr*bndry_flux_out_upr'/Δ(mesh,n_intervals(mesh)),
    )*u[idxdown,:]
    return nothing
end

fast_mul(A::LazyGenerator, B::LazyGenerator) = fast_mul(build_full_generator(A).B,B)
function fast_mul(A::LazyGenerator,x::Real)
    blocks = (x*A.blocks[i] for i in 1:4)
    boundary_flux = BoundaryFlux(
        OneBoundaryFlux(A.boundary_flux.upper.in*x,A.boundary_flux.upper.out*x),# upper
        OneBoundaryFlux(A.boundary_flux.lower.in*x,A.boundary_flux.lower.out*x) # lower
    )
    D = x*A.D
    return LazyGenerator(A.dq,blocks,boundary_flux,D)
end
fast_mul(x::Real,A::LazyGenerator) = fast_mul(A,x)

function show(io::IO, mime::MIME"text/plain", B::LazyGenerator)
    if VERSION >= v"1.6"
        show(io, mime, fast_mul(SparseArrays.SparseMatrixCSC{Float64,Int}(LinearAlgebra.I(size(B,1))),B))
    else
        show(io, mime, fast_mul(Matrix{Float64}(LinearAlgebra.I(size(B,1))),B))
    end
end

function getindex_interior(B::LazyGenerator,row::Int,col::Int)

    if abs(row-col) > n_bases_per_level(B.dq) + n_bases_per_cell(B.dq)
        return 0.0
    else
        k, i, p = _map_from_index_interior(row,B)
        l, j, q = _map_from_index_interior(col,B)
        C = rates(B.dq.model)
        if abs(k-l)>1 # cant move more than one level, skip free
            return 0.0
        elseif (k!=l)&&((i!=j)||(C[i]==0.0)) # can't move a level unless i=j
            return 0.0
        elseif (k==l+1)&&(C[i]<0.0) # lower diag block
            return abs(C[i])*B.blocks[1][p,q]/Δ(B.dq,k)
        elseif (k==l-1)&&(C[i]>0.0) # upper diag block
            return C[i]*B.blocks[4][p,q]/Δ(B.dq,k)
        elseif (k==l) # on diagonal
            if (i==j)
                if (C[i]<0.0)
                    (p==q) ? (return abs(C[i])*B.blocks[3][p,q]/Δ(B.dq,k)+B.dq.model.T[i,i]) : (return abs(C[i])*B.blocks[3][p,q]/Δ(B.dq,k))
                elseif (C[i]>0.0)
                    (p==q) ? (return C[i]*B.blocks[2][p,q]/Δ(B.dq,k)+B.dq.model.T[i,i]) : (return C[i]*B.blocks[2][p,q]/Δ(B.dq,k))
                elseif (p==q)
                    return B.dq.model.T[i,i]
                else
                    return 0.0
                end
            elseif (l==1)&&(k==1)&&(C[i]<0.0)&&(C[j]>0.0) # lwr boundary condition
                i_idx = i - sum(rates(B.dq.model.S[1:i-1]).>=0.0)
                return abs(C[i])*B.dq.model.P_lwr[i_idx,j]*B.boundary_flux.lower.in[p]*B.boundary_flux.lower.out[q]/Δ(B.dq,k) + B.dq.model.T[i,j]*B.D[p,q]
            elseif (l==n_intervals(B.dq))&&(k==n_intervals(B.dq))&&(C[i]>0.0)&&(C[j]<0.0) # upr boundary condition
                i_idx = i - sum(rates(B.dq.model.S[1:i-1]).<=0.0)
                return C[i]*B.dq.model.P_upr[i_idx,j]*B.boundary_flux.upper.in[p]*B.boundary_flux.upper.out[q]/Δ(B.dq,k) + B.dq.model.T[i,j]*B.D[p,q]
            elseif (typeof(B.dq.mesh)<:FRAPMesh)&&(membership(B.dq.model.S,i)!=membership(B.dq.model.S,j)) # + to - change
                return B.dq.model.T[i,j]*B.D[p,q]
            elseif (p==q)
                return B.dq.model.T[i,j]
            else
                return 0.0
            end
        else
            return 0.0
        end
    end
end

function getindex_out_boundary(B::LazyGenerator,row::Int,col::Int)
    (!_is_boundary_index(row,B))&&throw(DomainError(row,"row index does not correspond to a boundary"))
    l, j, q = _map_from_index_interior(col,B)

    model = B.dq.model
    C = rates(model)

    i = _map_from_index_boundary(row,B)
    if (l==1)&&(C[j]>0.0)&&_has_left_boundary(model.S,i)&&(row∈(1:N₊(B.dq)))
        v = model.T[i,j]*B.boundary_flux.lower.out[q]
    elseif (l==n_intervals(B.dq))&&(C[j]<0.0)&&_has_right_boundary(model.S,i)&&(row∈(size(B,2).-(N₋(B.dq)-1:-1:0)))
        v = model.T[i,j]*B.boundary_flux.upper.out[q]
    else
        v = 0.0
    end

    return v
end

function getindex_in_boundary(B::LazyGenerator,row::Int,col::Int)
    (!_is_boundary_index(col,B))&&throw(DomainError(col,"col index does not correspond to a boundary"))
    k, i, p = _map_from_index_interior(row,B)

    C = rates(B.dq)

    j = _map_from_index_boundary(col,B)
    if (k==1)&&(C[i]<0.0)&&(C[j]<=0.0)&&(col∈(1:N₋(B.dq)))
        i_idx = i - sum(rates(B.dq.model.S[1:i-1]).>=0.0)
        v = (B.dq.model.P_lwr[i_idx,j])*abs(C[i])*B.boundary_flux.lower.in[p]/Δ(B.dq,1)
    elseif (k==n_intervals(B.dq))&&(C[i]>0.0)&&(C[j]>=0.0)&&(col∈(size(B,2).-(N₊(B.dq)-1:-1:0)))
        i_idx = i - sum(rates(B.dq.model.S[1:i-1]).<=0.0)
        v = (B.dq.model.P_upr[i_idx,j])*abs(C[i])*B.boundary_flux.upper.in[p]/Δ(B.dq,k)
    else
        v = 0.0
    end

    return v
end

function getindex(B::LazyGenerator,row::Int,col::Int)
    sz = size(B)
    !((0<row<=sz[1])&&(0<col<=sz[2]))&&throw(BoundsError(B,(row,col)))

    if _is_boundary_index(row,B) && _is_boundary_index(col,B)
        i = _map_from_index_boundary(row,B)
        j = _map_from_index_boundary(col,B)
        if (row∈(1:N₋(B.dq)))&&(col∈(1:N₋(B.dq)))
            v = B.dq.model.T[i,j]
        elseif (row∈(size(B,1).-(N₊(B.dq)-1:-1:0)))&&(col∈(size(B,2).-(N₊(B.dq)-1:-1:0)))
            v = B.dq.model.T[i,j]
        else
            v = 0.0
        end
    elseif _is_boundary_index(col,B)
        v = getindex_in_boundary(B,row,col)
    elseif _is_boundary_index(row,B)
        v = getindex_out_boundary(B,row,col)
    else
        v = getindex_interior(B,row,col)
    end
    return v
end

# function getindex(B::LazyGenerator,i::Int)
#     !(0<i<=length(B))&&throw(BoundsError(B,i))
#     sz = size(B)
#     col = (i-1)÷sz[1]
#     row = i-col*sz[1]
#     col += 1
#     return B[row,col]
# end

max_stars_repo_name: angus-lewis/StochasticFluidQueues

max_issues_repo_name: angus-lewis/StochasticFluidQueues

max_forks_repo_name: angus-lewis/StochasticFluidQueues


text:
# # Tutorial
#
# Here we will see what are the basic building blocks of AlgebraOfGraphics, and how to
# combine them to create complex plots based on tables or other style formats.
#
# ## Basic building blocks
#
# The most important functions are `style`, and `spec`.
# `style` determines the mapping from data to plot. Its positional arguments correspond to
# the `x`, `y` or `z` axes of the plot, whereas the keyword arguments correspond to plot
# attributes that can vary continuously or discretely, such as `color` or `markersize`.
# The style is split according to the categorical attributes in it, and then styled using
# a default palette.
# Finally `spec` can be used to give style-independent specifications about the plot
# (plotting function or attributes).
#
# `style` and `spec` work in various context. In the following we will explore
# `DataContext`, which is introduced doing `data(df)` for any tabular style structure `df`.
# In this context, `style` and `group` accept symbols and integers, which correspond to
# columns of the data.
#
# ## Operations
#
# The outputs of `style`, `spec`, and `data` can be combined with `+` or `*`,
# to generate an `AlgebraicList` object, which can then be plotted using the
# function `draw`. The actual drawing is done by AbstractPlotting.
#
# The operation `+` is used to create separate layer. `a + b` has as many layers as `la + lb`,
# where `la` and `lb` are the number of layers in `a` and `b` respectively.
#
# The operation `a * b` create `la * lb` layers, where `la` and `lb` are the number of layers
# in `a` and `b` respectively. Each layer of `a * b` contains the combined information of
# the corresponding layer in `a` and the corresponding layer in `b`. In simple cases,
# however, both `a` and `b` will only have one layer, and `a * b` simply combines the
# information.
#
# ## Working with tables

using RDatasets: dataset
using AlgebraOfGraphics, AbstractPlotting, CairoMakie
mpg = dataset("ggplot2", "mpg");
cols = style(:Displ, :Hwy);
grp = style(color = :Cyl => categorical);
scat = spec(Scatter)
pipeline = cols * scat
data(mpg) * pipeline |> draw
AbstractPlotting.save("scatter.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](scatter.svg)
#
# Now let's simply add `grp` to the pipeline to style the color.

data(mpg) * grp * pipeline |> draw
AbstractPlotting.save("grouped_scatter.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](grouped_scatter.svg)
# Traces can be added together with `+`.

using AlgebraOfGraphics: linear
pipenew = cols * (scat + linear)
data(mpg) * pipenew |> draw
AbstractPlotting.save("linear.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](linear.svg)
# We can put grouping in the pipeline (we get a warning because of a degenerate group).

data(mpg) * grp * pipenew |> draw
AbstractPlotting.save("grouped_linear.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](grouped_linear.svg)
# This is a more complex example, where we split the scatter plot,
# but do the linear regression with all the style. Moreover, we pass weights to `linear`
# to compute the regression line with weighted least squares.

different_grouping = grp * scat + linear * style(wts=:Hwy)
data(mpg) * cols * different_grouping |> draw
AbstractPlotting.save("semi_grouped.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](semi_grouped.svg)
#
# Different analyses are also possible, always with the same syntax:

using AlgebraOfGraphics: smooth, density
data(mpg) * cols * grp * (scat + smooth(span = 0.8)) |> draw
AbstractPlotting.save("loess.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](loess.svg)

data(mpg) * cols * density |> draw
AbstractPlotting.save("density.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](density.svg)
#
# We can also add styling that only makes sense in one spec (e.g. `markersize`) by
# multiplying them:
#
newstyle = style(markersize = :Cyl) * spec(markersize = (0.1, 5))
data(mpg) * cols * (scat * newstyle + smooth(span = 0.8)) |> draw
AbstractPlotting.save("loess_markersize.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](loess_markersize.svg)
#
# ## Layout
#
# Thanks to the MakieLayout package it is possible to create plots where categorical variables
# inform the layout.

iris = dataset("datasets", "iris")
cols = style(:SepalLength, :SepalWidth)
grp = style(layout_x = :Species)
geom = spec(Scatter) + linear
data(iris) * cols * grp * geom |> draw
AbstractPlotting.save("layout.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](layout.svg)
#
# ## Non tabular style (slicing context)
#
# The framework is not specific to tables., but can be used in different contexts.
# For instance, `dims()` introduces a context where each entry of the array corresponds
# to a trace.

x = [-pi..0, 0..pi]
y = [sin cos] # We use broadcasting semantics on `tuple.(x, y)`.
dims() * style(x, y, color = dims(1), linestyle = dims(2)) |> draw
AbstractPlotting.save("functions.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](functions.svg)

using Distributions
IGpdf(μ, λ) = t -> pdf(InverseGaussian(μ, λ), t)
pdfs = IGpdf.(1:4, [6 10])
dims() * style(fill(0..5), pdfs, color = dims(1), linestyle = dims(2)) |> draw
AbstractPlotting.save("distributions.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](distributions.svg)
#
# More generally, one can pass arguments to `dims` to implement the
# "slices are series" approach.

s = dims(1) * style(rand(50, 3), rand(50, 3, 2))
grp = style(color = dims(2), layout_x = dims(3))
s * grp * spec(Scatter) |> draw
AbstractPlotting.save("arrays.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](arrays.svg)

# This approach can be used in combination with the tabular context to work
# with "wide" data, where grouping is done by column.

iris = dataset("datasets", "iris")
cols = style([:SepalLength, :SepalWidth], [:PetalLength :PetalWidth])
grp = style(layout_x = dims(1), layout_y = dims(2), color = :Species)
geom = spec(Scatter) + linear
data(iris) * cols * grp * geom |> draw
AbstractPlotting.save("layout_wide.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](layout_wide.svg)

ext: jl


text:
solvers_const = Dict(
    :ARCLDLt_abs => (HessDense, PDataLDLt, solve_modelARCDiagAbs, ()),
    :ARCLDLt_HO_Sham => (
        HessDense,
        PDataLDLt,
        (v, w, x, y, z) -> solve_modelARCDiag_HO(
            v,
            w,
            x,
            y,
            z,
            ho_correction = :Shamanskii,
            λfact = 100.0,
        ),
        (),
    ),
    :ARCLDLt_HO_vs_Nwt => (
        HessDense,
        PDataLDLt,
        (v, w, x, y, z) -> solve_modelARCDiag_HO_vs_Nwt(
            v,
            w,
            x,
            y,
            z,
            λfact = 100.0,
            nwt_res_fact = 0.8,
        ),
        (),
    ),
    :ARCLDLt => (HessDense, PDataLDLt, solve_modelARCDiag, ()),
    # :ARCMA57_abs => (HessSparse, PDataMA57, solve_modelARCDiagAbs, ()),
    # :ARCMA57_Sham_vs_Nwt_λ => (HessSparse, PDataMA57, (v, w, x, y, z) -> solve_modelARCDiag_HO_vs_Nwt(v, w, x, y, z, λfact = 1.0, nwt_res_fact = 0.25, ho_correction = :Shamanskii_MA57), ()),
    # :ARCMA57_Sham_vs_Nwt => (HessSparse, PDataMA57, solve_modelARCDiag, ()),
    # :ARCMA57_Sham_λ => (HessSparse, PDataMA57, (v, w, x, y, z) -> solve_modelARCDiag_HO(v, w, x, y, z, ho_correction = :Shamanskii_MA57, λfact = 1.0), ()),
    # :ARCMA57 => (HessSparse, PDataMA57, solve_modelARCDiag, ()),
    # :ARCMA97_abs => (HessDense, PDataMA97, solve_modelARCDiagAbs, ()),
    # :ARCMA97 => (HessDense, PDataMA97, solve_modelARCDiag, ()),
    :ARCqKdense => (
        HessDense,
        PDataKARC,
        solve_modelKARC,
        [(:shifts => 10.0 .^ (collect(-20.0:1.0:20.0)))],
    ),
    :ARCqKOp => (
        HessOp,
        PDataKARC,
        solve_modelKARC,
        [:shifts => 10.0 .^ (collect(-20.0:1.0:20.0))],
    ),
    :ARCqKsparse => (
        HessSparse,
        PDataKARC,
        solve_modelKARC,
        [:shifts => 10.0 .^ (collect(-20.0:1.0:20.0))],
    ),
    :ARCSpectral_abs => (HessDense, PDataSpectral, solve_modelARCDiagAbs, ()),
    :ARCSpectral => (HessDense, PDataSpectral, solve_modelARCDiag, ()),
    :ST_TRdense => (HessDense, PDataST, solve_modelST_TR, ()),
    :ST_TROp => (HessOp, PDataST, solve_modelST_TR, ()),
    :ST_TRsparse => (HessSparse, PDataST, solve_modelST_TR, ()),
    :TRKdense => (HessDense, PDataKTR, solve_modelKTR, ()),
    :TRKOp => (HessOp, PDataKTR, solve_modelKTR, ()),
    :TRKsparse => (HessSparse, PDataKTR, solve_modelKTR, ()),
    :TRLDLt_abs => (HessDense, PDataLDLt, solve_modelTRDiagAbs, ()),
    :TRLDLt_HO_vs_Nwt_λ => (
        HessDense,
        PDataLDLt,
        (v, w, x, y, z) -> solve_modelTRDiag_HO_vs_Nwt_λ(
            v,
            w,
            x,
            y,
            z;
            ho_correction = :Shamanskii,
            nwt_res_fact = 0.25,
            λfact = 1.0,
        ),
        (),
    ),
    :TRLDLt => (HessDense, PDataLDLt, solve_modelTRDiag, ()),
    # :TRMA57_abs => (HessSparse, PDataMA57, solve_modelTRDiagAbs, ()),
    # :TRMA57_Sham_vs_Nwt_λ => (HessSparse, PDataMA57, (v, w, x, y, z) -> solve_modelTRDiag_HO_vs_Nwt_λ(v, w, x, y, z; ho_correction = :Shamanskii_MA57, nwt_res_fact = 0.25, λfact = 1.0), ()),
    # :TRMA57 => (HessSparse, PDataMA57, solve_modelTRDiag, ()),
    # :TRMA97_abs => (HessDense, PDataMA97, solve_modelTRDiagAbs, ()),
    # :TRMA97 => (HessSparse, PDataMA97, solve_modelTRDiag, ()),
    :TRSpectral_abs => (HessDense, PDataSpectral, solve_modelTRDiagAbs, ()),
    :TRSpectral => (HessDense, PDataSpectral, solve_modelTRDiag, ()),
)


text:
type brack
    i::Vector{Int}
    w::Float64
end

convert(::Type{Tuple},x::Vector{Int})=ntuple(i->x[i],length(x))

*(a::brack,b::brack) = (x= deepcopy(a);push!(x.i,b.i[1]);x.w=a.w*b.w;x)

kron(a::Array{brack,1})= a

spzeros(n::Int) = SparseVector(n,Int[],Float64[])

function findbracket(v,x::Vector)
    i=searchsortedfirst(x,v)
    if i==1
        return [brack([1],.5);brack([1],.5)]
    else
        dx = x[i]-x[i-1]
        return [brack([i-1],(x[i]-v)/dx);brack([i],(v-x[i-1])/dx)]
    end
end

function updateT(M::Model)
    N = length(M.distribution.tG[1])
    Pf = Vector{Float64}[clamp!(M(s.name,hcat([g[:] for g in M.distribution.tG]...)),s.bounds...) for s in Endogenous(M.variables)]

    dims1  = ntuple(i->length(Endogenous(M.variables)[i].val),length(Endogenous(M.variables)))
    pdims1 = prod(dims1)
    Gdims = size(M.distribution.tG[1])

    # fill!(M.distribution.T,0.0)
    M.distribution.T.colptr = ones(Int,size(M.distribution.T,1)+1)
    M.distribution.T.nzval=Float64[]
    M.distribution.T.rowval=Float64[]

    Endo = Endogenous(M.variables)
    Stoc = Stochastic(M.variables)
    ne = length(Endo)
    for i = 1:N
        brackets = [findbracket(Pf[ie][i],Endo[ie].val.x) for ie = 1:length(Endo)]
        S1= collect(ind2sub(Gdims,i))[ne+1:end]
        abrackets = kron(brackets...)

        w = spzeros(pdims1)
        for b in abrackets
            w[sub2ind(dims1,b.i...)]+=b.w
        end
        for S2 = 1:length(Stoc)
            if isag(Stoc[S2])
                w = vcat([w[:]*p for p ∈ Stoc[S2].val.T[:,S1[S2]]]...)
            else
                w = vcat([w[:]*p for p ∈ Stoc[S2].val.T[S1[S2],:]]...)
            end
        end

        for j = 1:length(w.nzind)
            M.distribution.T[w.nzind[j],i] += w.nzval[j]
        end
    end
    return
end

function updateT1(M::Model)
    N = length(M.distribution.tG[1])
    Pf = Vector{Float64}[clamp!(M(s.name,hcat([g[:] for g in M.distribution.tG]...)),s.bounds...) for s in Endogenous(M.variables)]
    Gdims = size(M.distribution.tG[1])
    Sdim   = tuple([length(s) for s in Stochastic(M.variables)]...)

    M.distribution.T.colptr = ones(Int,size(M.distribution.T,1)+1)
    M.distribution.T.nzval=Float64[]
    M.distribution.T.rowval=Float64[]

    Endo = Endogenous(M.variables)
    Stoc = Stochastic(M.variables)
    ne = length(Endo)
    for i = 1:N
        brackets = [findbracket(Pf[ie][i],Endo[ie].val.x) for ie = 1:length(Endo)]
        S1 = Int[findfirst(s.val.x,M.distribution.tG[findfirst(names(State(M.variables)),s.name)][i]) for s ∈ Stoc]

        S1= collect(ind2sub(Gdims,i))[ne+1:end]

        abrackets = kron(brackets...)

        for j = 1:length(abrackets)
            for jj ∈ CartesianRange(Sdim)
                t = 1.0
                for S2 = 1:length(Stoc)
                    if isag(Stoc[S2])
                        t*= Stoc[S2].val.T[jj.I[S2],S1[S2]]
                    else
                        t*= Stoc[S2].val.T[S1[S2],jj.I[S2]]
                    end
                end
                M.distribution.T[sub2ind(Gdims,abrackets[j].i...,jj.I...),i]+=abrackets[j].w*t
            end
        end
    end
    return
end

function updatetd(M::Model)
    d0 = M.distribution.td[:]
    for j = 1:5000
        d1=M.distribution.T*d0
        mean(abs((d1[:]-d0[:])))<1e-10 && break
        d0[:]=d1[:]/sum(d1)
        if j==5000
            warn("Model distribution did not converge")
        end
    end

    M.distribution.td[:] = d0
    return
end

function updatedistribution(M::Model)
    updateT(M)
    updatetd(M)
end

function ∫(M::Model,v::Symbol,t=0)
    ∫(M,M[v,t])
end

function ∫(M::Model,v::Vector{Float64})
    ∫(M,reshape(M.G(v,hcat([vec(x) for x in M.distribution.tG]...)),size(M.distribution.tG[1])))
end

function ∫(M::Model,V::Array{Float64})
    d = V.*M.distribution.td
    if all(![isag(s) for s in State(M.variables)])
        return sum(d)
    else
        ns=length(State(M.variables))
        id = find([isag(s) for s in State(M.variables)])
        D= zeros(size(M.distribution.td,id...))
        for i = 1:prod(size(D))
            id1 = ind2sub(size(D),i)
            inds = Any[Colon() for i = 1:ns]
            cnt = 1
            for ii = 1:ns
                if [isag(s) for s in State(M.variables)][ii]
                    inds[ii] = id1[cnt]
                    cnt+=1
                end
            end

            D[id1...] = sum(d[inds...])/max(sum(M.distribution.td[inds...]),1e-16)
        end
        return D
    end
end

function updateAggregates(M::Model,ϕ=1.0)
    if any([isag(s) for s in State(M.variables)])
        agg=Aggregate(M.variables)
        for s ∈ Aggregate(M.variables)
            ag  = ∫(M,s.val.target[1],s.val.target[2])
            eag = filter(x->isag(x),State(M.variables))
            t = [vtype(e)==Endogenous ? -1 : 0 for e in eag]

            X   = ndgrid([e.val.x for e in eag]...)
            Xid = findfirst(State(M.variables,false),s)
            for ii = 1:length(X[1])
                id=BitArray(.*([M[eag[ie].name,t[ie]].==X[ie][ii] for ie = 1:length(eag)]...))
                M.X[id,Xid] *= (1-ϕ)
                M.X[id,Xid] += ϕ*ag[ii]*ones(sum(id))
            end
        end
    else
        for s ∈ Aggregate(M.variables)
            M.X[:,findfirst(State(M.variables,false),s)] *= (1-ϕ)
            M.X[:,findfirst(State(M.variables,false),s)] += ϕ*∫(M,s.val.target[1],s.val.target[2])
        end
    end
end

function updateA(M,ϕ=1.0)
    updatedistribution(M)
    M.temp.updateD(M)
    updateAggregates(M,ϕ)
end


text:
using Documenter, SLIMutils

makedocs(
    sitename = "SLIMutils Reference",
    modules = [SLIMutils]
)

deploydocs(
    repo = "github.com/slimgroup/SLIMutils.jl.git"
)


text:
abstract type ConstructionFlag end
struct TrustInput <: ConstructionFlag end

"""
`EventSeries` holds a vector of `values` toghether with its corresponding
`timestamps`. `EventSeries` is subtype `AbstractVector`, and of both `timestamps`
and `values` must be subtypes of `AbstractVector`.

EventSeries(timestamps::U, values::W, ::TrustInput)

    - Assumes that `timestamps` and `values` are of equal
      length, and that `timestamps` is sorted.

EventSeries(timestamps, values; drop_repeated=true, keep_end=true)

    - Checks that `timestamps` and `values` are of equal length
    - Checks that `timestamps` is sorted
    - Throws `AssertionError` if any of the above mentioned checks failss.
    - If `drop_repeated` is true (default) all repeated values in values,
      and corresponding timestamps, are removed before construction.
    - If `keep_end` is true (default), the last value in values (with its
      corresponding timestamp) is kept regardless of whether it equals the
      previous value. Keep_end will only be considered if `drop_repeated`
      is true.
"""
struct EventSeries{T, U, V, W} <: AbstractVector{Tuple{T, V}}
    timestamps::U
    values::W
end

function EventSeries(
    timestamps::U,
    values::W,
    ::TrustInput
    ) where {T, U<:AbstractVector{T}, V,  W<:AbstractVector{V}}
    EventSeries{T, U, V, W}(timestamps, values)
end

function EventSeries(timestamps, values; drop_repeated=true, keep_end=true)
    @assert issorted(timestamps)
    @assert length(timestamps) == length(values)
    if drop_repeated
        select = [true; [a!=b for (a,b) in neighbors(values)]]
        keep_end && (select[end] = true)
        timestamps, values = all(select) ? (timestamps, values) : (timestamps[select],  values[select])
    end
    EventSeries(timestamps, values, TrustInput())
end

Base.size(ts::EventSeries) = size(ts.timestamps)
Base.getindex(ts::EventSeries, i::Number) = Event(timestamp=ts.timestamps[i], value=ts.values[i])
Base.getindex(ts::EventSeries, I) = EventSeries(ts.timestamps[I], ts.values[I])
Base.IndexStyle(ts::EventSeries) = IndexStyle(ts.timestamps)

function Base.push!(ts::EventSeries, e::Event)
    isempty(ts.timestamps) || @assert last(ts.timestamps) <= e.timestamp
    push!(ts.timestamps, e.timestamp)
    push!(ts.values, e.value)
end

function Base.append!(ts1::EventSeries, ts2::EventSeries)
    @assert last(ts1.timestamps) <= first(ts2.timestamps)
    append!(ts1.timestamps, ts2.timestamps)
    append!(ts1.values, ts2.values)
end

"""
    timestamptype(y::EventSeries)

Returns time type of the EventSeries, which equals `eltype(y.timestamps)`.
"""
timestamptype(::EventSeries{T}) where {T} = T

"""
`TaggedEventSeries{T}` holds a mapping from `tag::T` to an `eventseries::EventSeries`.

Assumes that the `timestamptype(eventseries)` can be compared betweeen types.

# Constructors

TaggedEventSeries(d::Dict{Symbol, EventSeries})
TaggedEventSeries()

"""
struct TaggedEventSeries
    data::Dict{Symbol, EventSeries}
end
function TaggedEventSeries(;kwargs...)
    tes = TaggedEventSeries(Dict{Symbol, EventSeries}())
    for (k,v) in kwargs
        tes[k] = v
    end
    tes
end

Base.length(tts::TaggedEventSeries) = sum(length(v) for v in values(tts.data))
Base.size(tts::TaggedEventSeries) = (length(tts),)

Base.getindex(tts::TaggedEventSeries, args...) = getindex(tts.data, args...)
Base.setindex!(tts::TaggedEventSeries, args...) = setindex!(tts.data, args...)

function Base.push!(tts::TaggedEventSeries, e::TaggedEvent)
    if e.tag in keys(tts.data)
        push!(tts[e.tag], Event(e.timestamp, e.value))
    else
        tts[e.tag] = EventSeries([e.timestamp], [e.value])
    end
end

function EventSeries(tts::TaggedEventSeries)
    timestamps_ = timestamps(tts) |> collect
    values_ = fill_forward_value.(Ref(tts), timestamps_)
    EventSeries(timestamps_, values_)
end

"""
    sorted_tag_idx(tts::TaggedEventSeries)

Returns a time sorted iterator over all (tag, index_in_event_series) pairs in the
input TaggedEventSeries.
"""
function sorted_tag_idx(tts::TaggedEventSeries)
    Channel() do c
        state = Dict(t=>(tts[t].timestamps[1], 1, length(tts[t])) for t in keys(tts.data))
        while !isempty(state)
            (time, idx, len), tag = findmin(state)
            put!(c, (tag, idx))
            if idx == len
                delete!(state, tag)
            else
                idx = idx + 1
                state[tag] = (tts[tag].timestamps[idx], idx, len)
            end
        end
    end
end

"""
    tagged_events_itr(tts::TaggedEventSeries)

Returns a time sorted iterator over all TaggedEvents in the input TaggedEventSeries.

see also `tagged_events(tts::TaggedEventSeries)`.
"""
tagged_events_itr(tts::TaggedEventSeries) = (TaggedEvent(tag, tts[tag][idx]) for (tag, idx) in sorted_tag_idx(tts))
timestamps(tts::TaggedEventSeries) = unique(tts[tag].timestamps[idx] for (tag, idx) in sorted_tag_idx(tts))

"""
    tagged_events(tts::TaggedEventSeries)

Returns a time sorted vector of all TaggedEvents in the input TaggedEventSeries.

see also tagged_events_itr(tts::TaggedEventSeries)
"""
function tagged_events(tts::TaggedEventSeries)
    sort(vcat([[TaggedEvent(tag, e) for e in ts] for (tag, ts) in tts.data]...), by=x->x.timestamp)
end

function fill_forward_event(ts::EventSeries, time)
    t0, t1 =  first(ts.timestamps), last(ts.timestamps)
    if time <= t0
        return time < t0 ? nothing : ts[1]
    elseif t1 <= time
        return ts[end]
    else
        idx = searchsortedfirst(ts.timestamps, time)
        return time == ts.timestamps[idx] ? ts[idx] : ts[idx-1]
    end
end
fill_forward_event(tts::TaggedEventSeries, time) = _fill_forward_event(tts, time, sort(collect(keys(tts.data))))
_fill_forward_event(tts::TaggedEventSeries, time, skeys) = NamedTuple{Tuple(skeys)}(Tuple(fill_forward_event(tts[tag], time) for tag in skeys))

"""
    fill_forward_event(ts::EventSeries, time)

Returns the prior (with respect to `time`) `Event` in the input `EventSeries`).

    fill_forward_event(ts::TaggedEventSeries, time)

Returns a named tuple with the prior (with respect to `time`) `Event`s for each
tag in the input `TaggedEventSeries`.

If `time` is before the first timestamp in `EventSeries` / `TaggedEventSeries`,
then `nothing` is returned.
"""
fill_forward_event

fill_forward_value(ts::EventSeries, time) = value(fill_forward_event(ts, time))
fill_forward_value(tts::TaggedEventSeries, time) = _fill_forward_value(tts, time, sort(collect(keys(tts.data))))
_fill_forward_value(tts::TaggedEventSeries, time, skeys) = NamedTuple{Tuple(skeys)}(Tuple(fill_forward_value(tts[tag], time) for tag in skeys))

"""
    fill_forward_value(ts::EventSeries, time)

Returns the prior (with respect to `time`) `value` in the input `EventSeries`).

    fill_forward_value(ts::TaggedEventSeries, time)

Returns a named tuple with the prior (with respect to `time`) `values`s for each
tag in the input `TaggedEventSeries`.

If `time` is before the first timestamp in `EventSeries` / `TaggedEventSeries`,
then `nothing` is returned.
"""
fill_forward_value

"""
     select(y::EventSeries, t1, t2)

returns an EventSeries which is a subset of the input series, containing the `Event`s
in the time domain [tstart, tend]. The endpoint values are set by filling forward
Assumes that the input time domain `[t1, t2]` is contained in the input EventSeries.
"""
function select(y::EventSeries, t1, t2)
    @assert y.timestamps[1] <= t1 < t2 <= y.timestamps[end] "invaid time limits"
    v1 = fill_forward_value(y, t1)
    v2 = fill_forward_value(y, t2)
    select = [t1 < e.timestamp < t2 for e in y]
    EventSeries([t1; y.timestamps[select]; t2], [v1; y.values[select]; v2])
end

"""
    align(ys::EventSeries{T}...)

Returns a tuple of EventSeries containing subsets of the corresponding series input,
such that the time domain of each output series corresponds to the largest common time
domain of the input series. See `select`.
"""
function align(ys::EventSeries{T}...) where T
    tmin = maximum(y.timestamps[1] for y in ys)
    tmax = minimum(y.timestamps[end] for y in ys)
    Tuple(select(y, tmin, tmax) for y in ys)
end

"""
    cumtime(es::EventSeries{T}, val)

Returns cummulative time of EventSeries where the value equals `val`
"""
function cumtime(es::EventSeries{T}, fun) where T
    eventpairs = SparseTimeSeries.neighbors(es)
    durations = (e2.timestamp-e1.timestamp for (e1,e2) in eventpairs if fun(e1.value))
    isempty(durations) ? zero(T) : sum(durations)
end

cumtime(es::EventSeries{T}, fun, t1, t2) where {T} = cumtime(select(es,t1,t2), fun)

function fuse(;kwargs...)
    tes = TaggedEventSeries(;kwargs...)
    EventSeries(tes)
end

function splice(;kwargs...) where {T}
    tes = TaggedEventSeries(;kwargs...)
    taggedevents = tagged_events(tes)
    EventSeries([e.timestamp for e in taggedevents], [e.tag=>e.value for e in taggedevents])
end

function Base.filter(f, y::EventSeries)
    select = [f(e) for e in y]
    EventSeries(y.timestamps[select], y.values[select])
end

segments(y::EventSeries) = (Segment(e1.value, e1.timestamp, e2.timestamp) for (e1, e2) in neighbors(y))

function next(ts::EventSeries, time; include_boundary=false)
    t0, t1 =  first(ts.timestamps), last(ts.timestamps)
    if time > t1
        return nothing
    elseif time == t1
        return include_boundary ? ts[end] : nothing
    else
        idx = searchsortedfirst(ts.timestamps, time)
        idx = !include_boundary && in(time, ts.timestamps) ? idx+1 : idx
        ts[idx]
    end

end

function previous(ts::EventSeries, time; include_boundary=false)
    t0, t1 =  first(ts.timestamps), last(ts.timestamps)
    if time < t0
        return nothing
    elseif time == t0
        return include_boundary ? ts[1] : nothing
    else
        idx = searchsortedfirst(ts.timestamps, time)
        if in(time, ts.timestamps) && include_boundary
            return ts[idx]
        else
            return ts[idx-1]
        end
    end

end


text:
export Mesh, link_triangles_to_tetrahedra!, compute_size!, find_tetrahedron_containing_point, keep!, collect_lines!, unify!, extend_mesh, vtk_write
export vtk_write_tri, get_surface_points, get_normal_vectors, generate_field
export octosplit, color_domains


text:
fs1()

@testset "unixify" begin
    @test F.unixify("")         == "/"
    @test F.unixify("blah.txt") == "blah.txt"
    @test F.unixify("blah/")    == "blah/"
    @test F.unixify("foo/bar")  == "foo/bar/"
end

@testset "join_rpath" begin
    @test F.join_rpath("blah/blah") == joinpath("blah", "blah")
end

@testset "parse_rpath" begin
    F.PATHS[:folder] = "fld"
    # /[path]
    @test_throws F.RelativePathError F.parse_rpath("/")
    @test F.parse_rpath("/a") == "/a"
    @test F.parse_rpath("/a/b") == "/a/b"
    @test F.parse_rpath("/a/b", canonical=true) == joinpath("fld", "a", "b")
    @test F.parse_rpath("/a/../b", canonical=true) == joinpath("fld", "b")

    # ./[path]
    set_curpath(joinpath("pages", "pg1.md"))
    F.PATHS[:assets] = joinpath("fld", "assets")
    @test_throws F.RelativePathError F.parse_rpath("./")
    @test F.parse_rpath("./a") == "/assets/pages/pg1/a"
    @test F.parse_rpath("./a/b") == "/assets/pages/pg1/a/b"
    @test F.parse_rpath("./a/b", canonical=true) == joinpath("fld", "assets", "pages", "pg1", "a", "b")

    # [path]
    @test_throws F.RelativePathError F.parse_rpath("")
    @test F.parse_rpath("blah") == "/assets/blah"
    @test F.parse_rpath("blah", code=true) == "/assets/pages/pg1/code/blah"
    @test F.parse_rpath("blah", canonical=true) == joinpath("fld", "assets", "blah")
end

@testset "resolve_rpath" begin
    root = F.PATHS[:folder] = mktempdir()
    ass  = F.PATHS[:assets] = joinpath(root, "assets")
    mkdir(ass)
    write(joinpath(ass, "p1.jl"), "a = 5")
    write(joinpath(ass, "p2.png"), "gibberish")

    fp, d, fn = F.resolve_rpath("/assets/p1", "julia")
    @test fp == joinpath(ass, "p1.jl")
    @test d == ass
    @test fn == "p1"

    fp, d, fn = F.resolve_rpath("/assets/p2")
    @test fp == joinpath(ass, "p2.png")
    @test d == ass
    @test fn == "p2"

    @test_throws F.FileNotFoundError F.resolve_rpath("/assets/foo")
    @test_throws F.FileNotFoundError F.resolve_rpath("foo")
    @test_throws F.FileNotFoundError F.resolve_rpath("foo", "julia")
end

@testset "form_cpaths" begin
    root = F.PATHS[:folder] = mktempdir()
    ass  = F.PATHS[:assets] = joinpath(root, "assets")
    mkdir(ass)
    curpath = set_curpath(joinpath("pages", "pg1.md"))

    write(joinpath(ass, "p1.jl"), "a = 5")
    write(joinpath(ass, "p2.png"), "gibberish")

    sp, sd, sn, od, op, rp = F.form_codepaths("p1")
    @test sp == joinpath(ass, splitext(curpath)[1], "code", "p1.jl")
    @test sd == joinpath(ass, splitext(curpath)[1], "code")
    @test sn == "p1.jl"
    @test od == joinpath(ass, splitext(curpath)[1], "code", "output")
    @test op == joinpath(od, "p1.out")
    @test rp == joinpath(od, "p1.res")
end


text:
using Turing, ContinuousBenchmarks

include(joinpath(ContinuousBenchmarks.STAN_DATA_DIR, "bernoulli-stan.data.jl"))
include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "bernoulli.model.jl"))

@tbenchmark(HMC(10, 0.25, 5), bermodel, berstandata[1]["y"])
bench_res = @tbenchmark(HMC(10000, 0.25, 5), bermodel, berstandata[1]["y"])
logd = build_logd("Bernoulli Model", bench_res...)

#include(joinpath(ContinuousBenchmarks.BENCH_DIR, "bernoulli-stan.run.jl"))
#logd["stan"] = Dict("theta" => mean(theta_stan))
#logd["time_stan"] = ber_time

print_log(logd)
send_log(logd)


text:
cxx"""
static cv::UMat image;
static bool backprojMode = false;
static bool selectObject = false;
static int trackObject = 0;
static bool showHist = true;
static cv::Rect selection;
static int vmin = 10, vmax = 256, smin = 30;
//int argc = 0;
//char** argv;
//cv::String keys;
//cv::CommandLineParser parser(int argc, const char** argv, cv::String &keys);

static void onMouse(int event, int x, int y, int, void*)
{
    static cv::Point origin;

    if (selectObject)
    {
        selection.x = std::min(x, origin.x);
        selection.y = std::min(y, origin.y);
        selection.width = std::abs(x - origin.x);
        selection.height = std::abs(y - origin.y);

        selection &= cv::Rect(0, 0, image.cols, image.rows);
    }

    switch (event)
    {
    case cv::EVENT_LBUTTONDOWN:
        origin = cv::Point(x, y);
        selection = cv::Rect(x, y, 0, 0);
        selectObject = true;
        break;
    case cv::EVENT_LBUTTONUP:
        selectObject = false;
        if (selection.width > 0 && selection.height > 0)
            trackObject = -1;
        break;
    default:
        break;
    }
}

static void help()
{
    std::cout << "\nThis is a demo that shows mean-shift based tracking using Transparent API\n"
            "You select a color objects such as your face and it tracks it.\n"
            "This reads from video camera (0 by default, or the camera number the user enters\n"
            "Usage: \n"
            "   ./camshiftdemo [camera number]\n";

    std::cout << "\n\nHot keys: \n"
            "\tESC - quit the program\n"
            "\ts - stop the tracking\n"
            "\tb - switch to/from backprojection view\n"
            "\th - show/hide object histogram\n"
            "\tp - pause video\n"
            "\tc - use OpenCL or not\n"
            "To initialize tracking, select the object with mouse\n";
}

int ocl_camshift()
{

  help();

    cv::VideoCapture cap;
    cv::Rect trackWindow;
    int hsize = 16;
    float hranges[2] = { 0, 180 };

    //const char * const keys = { "{@camera_number| 0 | camera number}" };
    int camNum = 0; //parser.get<int>(0)

    cap.open(camNum);

    if (!cap.isOpened())
    {
        help();

        std::cout << "***Could not initialize capturing...***\n";
        std::cout << "Current parameter's value: \n";
        //parser.printMessage();

        return EXIT_FAILURE;
    }

    cv::namedWindow("Histogram", cv::WINDOW_NORMAL);
    cv::namedWindow("CamShift Demo", cv::WINDOW_NORMAL);
    cv::setMouseCallback("CamShift Demo", onMouse);
    cv::createTrackbar("Vmin", "CamShift Demo", &vmin, 256);
    cv::createTrackbar("Vmax", "CamShift Demo", &vmax, 256);
    cv::createTrackbar("Smin", "CamShift Demo", &smin, 256);

    cv::Mat frame, histimg(200, 320, CV_8UC3, cv::Scalar::all(0));
    cv::UMat hsv, hist, hue, mask, backproj;
    bool paused = false;

    for ( ; ; )
    {
        if (!paused)
        {
            cap >> frame;
            if (frame.empty())
                break;
        }

        frame.copyTo(image);

        if (!paused)
        {
            cv::cvtColor(image, hsv, cv::COLOR_BGR2HSV);

            if (trackObject)
            {
                int _vmin = vmin, _vmax = vmax;

                cv::inRange(hsv, cv::Scalar(0, smin, std::min(_vmin, _vmax)),
                        cv::Scalar(180, 256, std::max(_vmin, _vmax)), mask);

                int fromTo[2] = { 0,0 };
                hue.create(hsv.size(), hsv.depth());
                cv::mixChannels(std::vector<cv::UMat>(1, hsv), std::vector<cv::UMat>(1, hue), fromTo, 1);

                if (trackObject < 0)
                {
                    cv::UMat roi(hue, selection), maskroi(mask, selection);
                    cv::calcHist(std::vector<cv::Mat>(1, roi.getMat(cv::ACCESS_READ)), std::vector<int>(1, 0),
                                 maskroi, hist, std::vector<int>(1, hsize), std::vector<float>(hranges, hranges + 2));
                    cv::normalize(hist, hist, 0, 255, cv::NORM_MINMAX);

                    trackWindow = selection;
                    trackObject = 1;

                    histimg = cv::Scalar::all(0);
                    int binW = histimg.cols / hsize;
                    cv::Mat buf (1, hsize, CV_8UC3);
                    for (int i = 0; i < hsize; i++)
                        buf.at<cv::Vec3b>(i) = cv::Vec3b(cv::saturate_cast<uchar>(i*180./hsize), 255, 255);
                    cv::cvtColor(buf, buf, cv::COLOR_HSV2BGR);

                    {
                        cv::Mat _hist = hist.getMat(cv::ACCESS_READ);
                        for (int i = 0; i < hsize; i++)
                        {
                            int val = cv::saturate_cast<int>(_hist.at<float>(i)*histimg.rows/255);
                            cv::rectangle(histimg, cv::Point(i*binW, histimg.rows),
                                       cv::Point((i+1)*binW, histimg.rows - val),
                                       cv::Scalar(buf.at<cv::Vec3b>(i)), -1, 8);
                        }
                    }
                }

                cv::calcBackProject(std::vector<cv::UMat>(1, hue), std::vector<int>(1, 0), hist, backproj,
                                    std::vector<float>(hranges, hranges + 2), 1.0);
                cv::bitwise_and(backproj, mask, backproj);

                cv::RotatedRect trackBox = cv::CamShift(backproj, trackWindow,
                                    cv::TermCriteria(cv::TermCriteria::EPS | cv::TermCriteria::COUNT, 10, 1));
                if (trackWindow.area() <= 1)
                {
                    int cols = backproj.cols, rows = backproj.rows, r = (std::min(cols, rows) + 5)/6;
                    trackWindow = cv::Rect(trackWindow.x - r, trackWindow.y - r,
                                       trackWindow.x + r, trackWindow.y + r) &
                                  cv::Rect(0, 0, cols, rows);
                }

                if (backprojMode)
                    cv::cvtColor(backproj, image, cv::COLOR_GRAY2BGR);

                {
                    cv::Mat _image = image.getMat(cv::ACCESS_RW);
                    cv::ellipse(_image, trackBox, cv::Scalar(0, 0, 255), 3, cv::LINE_AA);
                }
            }
        }
        else if (trackObject < 0)
            paused = false;

        if (selectObject && selection.width > 0 && selection.height > 0)
        {
            cv::UMat roi(image, selection);
            cv::bitwise_not(roi, roi);
        }

        cv::imshow("CamShift Demo", image);
        if (showHist)
            cv::imshow("Histogram", histimg);

        char c = (char)cv::waitKey(10);
        if (c == 27)
            break;

        switch(c)
        {
        case 'b':
            backprojMode = !backprojMode;
            break;
        case 't':
            trackObject = 0;
            histimg = cv::Scalar::all(0);
            break;
        case 'h':
            showHist = !showHist;
            if (!showHist)
                cv::destroyWindow("Histogram");
            else
                cv::namedWindow("Histogram", cv::WINDOW_AUTOSIZE);
            break;
        case 'p':
            paused = !paused;
            break;
        case 'c':
            cv::ocl::setUseOpenCL(!cv::ocl::useOpenCL());
        default:
            break;
        }
    }
    return EXIT_SUCCESS;
}
"""


text:
using QuantumOpticsBase
using Test
using Random

@testset "spinors" begin

Random.seed!(0)

D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2)))

b1a = GenericBasis(2)
b1b = GenericBasis(3)
b2a = GenericBasis(1)
b2b = GenericBasis(4)
b3a = GenericBasis(1)
b3b = GenericBasis(5)

b_l = b1a⊕b2a⊕b3a
b_r = b1b⊕b2b⊕b3b

op1a = randoperator(b1a, b1b)
op1b = randoperator(b1a, b1b)
op2a = randoperator(b2a, b2b)
op2b = randoperator(b2a, b2b)
op3a = randoperator(b3a, b3b)
op123 = op1a ⊕ op2a ⊕ op3a
@test op123.basis_l == b_l
@test op123.basis_r == b_r

# Associativity
@test 1e-13 > D((op1a ⊕ op2a) ⊕ op3a, op1a ⊕ (op2a ⊕ op3a))
@test 1e-13 > D(op1a ⊕ op2a ⊕ op3a, op1a ⊕ (op2a ⊕ op3a))

# Mixed-product property
@test 1e-13 > D((op1a ⊕ op2a) * dagger(op1b ⊕ op2b), (op1a*dagger(op1b)) ⊕ (op2a*dagger(op2b)))

# Transpose
@test 1e-13 > D(dagger(op1a ⊕ op2a), dagger(op1a) ⊕ dagger(op2a))
@test 1e-13 > D(dagger(op1a ⊕ op2a), dagger(op1a) ⊕ dagger(op2a))

# Sparse version
@test isa(sparse(op1a)⊕sparse(op2a)⊕sparse(op3a), SparseOpType)

# Test lazy implementation
L = LazyDirectSum(op1a,op2a)
@test 1e-13 > D(op1a⊕op2a, L)
@test 1e-13 > D(2*(op1a⊕op2a), L+L)
@test isa(L ⊕ op3a, LazyDirectSum)
@test 1e-13 > D(op1a ⊕ op2a ⊕ op3a, L ⊕ op3a)
@test 1e-13 > D((op1a⊕op2a)*dagger(op1a⊕op2a), L*L')

# Test embedding
@test embed(b_l,b_r,1,op1a) == op1a ⊕ SparseOperator(b2a,b2b) ⊕ SparseOperator(b3a,b3b)
@test embed(b_l,b_r,[1,3],op1a ⊕ op3a) == embed(b_l,b_r,[1,3],(op1a,op3a)) == op1a ⊕ SparseOperator(b2a,b2b) ⊕ op3a

# Test with off-diagonal blocks
op = op1a ⊕ op3a
op_upper = randoperator(b1a,b3b)
setblock!(op, op_upper, 1, 2)
op_lower = randoperator(b1b,b3a)
setblock!(op, op_lower', 2, 1)

op_tot = op1a ⊕ SparseOperator(b2a,b2b) ⊕ op3a
setblock!(op_tot, op_upper, 1, 3)
setblock!(op_tot, op_lower', 3, 1)
@test op_tot == embed(b_l,b_r,[1,3],op)

# Lazy embed
@test embed(b_l,b_r,[1,2],L) == L ⊕ SparseOperator(b3a,b3b)

### Test example

# Define x-space
Nunitcells=2
Lx=Nunitcells/2
x_min=0
x_max = Lx
x_steps = 4

# Define y-space
Ly=Nunitcells/2
y_min = 0
y_max = Ly
y_steps = 8

# Define Bases
b_x = PositionBasis(x_min, x_max, x_steps)
b_y = PositionBasis(y_min, y_max, y_steps)
bcomp_x = b_x ⊗ b_y
b_sum = bcomp_x ⊕ bcomp_x

b_mom = MomentumBasis(b_x)
b_mom_y = MomentumBasis(b_y)
bcomp_p = MomentumBasis(b_x) ⊗ MomentumBasis(b_y)
b_sum_p = bcomp_p ⊕ bcomp_p

# Define Operators
H_up = randoperator(bcomp_x)
H_down = randoperator(b_x⊗b_y)
Hkin_x = randoperator(bcomp_x)
H = (H_up ⊕ H_down)

# Off-diagonal blocks - assign by hand
Ω_R = randoperator(bcomp_x)
nn = length(H_up.basis_l)
H.data[1:nn,nn+1:end] = Ω_R.data
H.data[nn+1:end,1:nn] = Ω_R.data'

H_h_upperline = hcat(Matrix(H_up.data),Matrix(Ω_R.data))
H_h_lowerline = hcat(Matrix(Ω_R.data'),Matrix(H_down.data))
H_matr=vcat(H_h_upperline,H_h_lowerline)
@test H.data ≈ H_matr
@test getblock(H, 1,2) == Ω_R

# States
ψ0_x=randstate(b_x)
ψ0_y=randstate(b_y)
ψ0_up=ψ0_x⊗ψ0_y
# Build with directsum (\oplus)
ψ0=(ψ0_x⊗ψ0_y)⊕(ψ0_x⊗ψ0_y)

# Test initial state
@test ψ0.data ≈ [ψ0_up.data;ψ0_up.data]
@test getblock(ψ0, 1) == ψ0_up == getblock(ψ0, 2)

# Test FFTs
Txp = transform(bcomp_x, bcomp_p; ket_only=true) # "normal" FFT operators as implemented
Tpx = transform(bcomp_p, bcomp_x; ket_only=true)
ψ0_p1 = (Tpx*ψ0_up)⊕(Tpx*ψ0_up) # Build with usual FFTs
ψ0_p2 = (Tpx⊕Tpx)*ψ0
@test ψ0_p1.data ≈ ψ0_p2.data

# Test off-diagonal blocks
H_test = copy(H)
setblock!(H_test, Ω_R, 1, 2)
setblock!(H_test, Ω_R', 2, 1)
@test H == H_test

# Lazy formulation
Txp_tot = Txp ⊕ Txp
Tpx_tot = Tpx ⊕ Tpx
Hkin_p = Tpx*Hkin_x*Txp
Hkin_tot = Hkin_p ⊕ Hkin_p
H_lazy = LazySum(H, LazyProduct(Txp_tot, Hkin_tot, Tpx_tot))

Hkin_x_tot = Txp_tot*Hkin_tot*Tpx_tot
@test dense(H_lazy).data ≈ dense(H + Hkin_x_tot).data
@test (H_lazy*ψ0).data ≈ (H*ψ0 + Hkin_x_tot*ψ0).data

# # Test time evolution
# tout, ψt = timeevolution.schroedinger([0.0:0.1:1.0;], ψ0, H)
# tout, ψt_lazy = timeevolution.schroedinger([0.0:0.1:1.0;], ψ0, H_lazy)
#
# for i=1:length(tout)
#     @test ψt[i].data ≈ ψt_lazy[i].data
# end

end # testset


text:
using HPAT

@acc hpat function udf_test(file_name)
    df = DataSource(DataTable{:id=Int64,:x=Float64,:y=Float64}, HDF5, file_name)
    df2 = aggregate(df, :id, :sx = sum(2*:x), :sy = sum(2*:y))
    return sum(df2[:sx])
end

udf_test(ENV["HOME"]*"/tmp/udf.hdf5")

max_stars_repo_name: ehsantn/data-analytics-benchmarks

max_issues_repo_name: ehsantn/data-analytics-benchmarks

max_forks_repo_name: ehsantn/data-analytics-benchmarks


text:
struct ARIMAParams{T} <: AbstractParams
    c::Bool
    p::Vector{T}
    d::T
    q::Vector{T}
    function ARIMAParams(c::Bool, p::AbstractVector{T}, d::T, q::AbstractVector{T}) where T <: Integer
        return new{T}(c, p, d, q)
    end
end

ARIMAParams(c::Bool,p::Integer,d::Integer,q::Integer) = ARIMAParams(c, collect(1:p), d, collect(1:q))
ARIMAParams(p::Integer,d::Integer,q::Integer) = ARIMAParams(true, p, d, q)

function toarma(params::ARIMAParams, arma::ARMAModel, μ)
    P = isempty(params.p) ? 0 : maximum(params.p)
    Q = isempty(params.q) ? 0 : maximum(params.q)
    d = params.d
    ϕdp = Polynomial([1;.-arma.ϕ]) * Polynomial(1 - variable())^d
    ϕ0 = .-coeffs(ϕdp)[2:end]
    μ0 = params.c ? μ * (1 - sum(ϕ0)) : zero(typeof(μ))
    return ARMAModel{P + d,Q}(μ0, ϕ0, arma.θ, arma.σ2)
end

toarma(params::ARIMAParams, arma::ARMAModel) = toarma(params, arma, arma.μ)

function fit(params::ARIMAParams, z::AbstractVector; n::Integer=2)
    length(z) > params.d || throw(DomainError("ARIMA requires at least d data points"))
    arima = fit(ARMAParams(params.c, params.p, params.q),  difference(z, d=params.d), n=n)
    return toarma(params, arima)
end


text:
import Theta: fgsm_s

@testset "FGSM" begin
    τ = [0.17012929+3.36021399im 0.42793573+1.65263414im 0.55503983+2.1393021im 0.49182605+2.17345553im 0.93865305+1.81009273im; 0.42793573+1.65263414im 0.87402978+1.33203047im 0.53489381+0.68238397im 0.57601957+1.20509707im 0.50028888+1.11605622im; 0.55503983+2.1393021im 0.53489381+0.68238397im 0.21882206+2.28319153im 0.74421992+1.51352273im 0.41855281+0.87400937im; 0.49182605+2.17345553im 0.57601957+1.20509707im 0.74421992+1.51352273im 0.51070570+1.62305859im 0.52380319+1.15073893im; 0.93865305+1.81009273im 0.50028888+1.11605622im 0.41855281+0.87400937im 0.52380319+1.15073893im 0.79948032+1.38925759im];
    R = RiemannMatrix(τ);
    @test begin
        z = zeros(5);
        theta_constants = [[theta(z, R, char=c) for c in chars_s] for chars_s in fgsm_chars()[1]];
        p = [sqrt(prod(theta_constant)) for theta_constant in theta_constants];
        abs(sum(p)) >= fgsm_s(R, fgsm_chars()[1])
        end
    @test fgsm(τ) == fgsm(R)
    @test fgsm_chars()[1][1] == [[[0,0,0,0,0], [0,0,0,0,0]], [[0,0,0,0,0], [1,1,1,1,1]], [[0,0,1,1,0], [0,1,0,0,0]], [[0,0,1,1,0], [1,0,1,1,1]], [[1,1,0,0,0], [0,0,0,1,0]], [[1,1,0,0,0], [1,1,1,0,1]], [[1,1,1,1,0], [0,1,0,1,0]], [[1,1,1,1,0], [1,0,1,0,1]]]
end


text:
function testcall(userinput; verbose=true)

    testname = match(r"t(?:est)?\((.*)\)", userinput).captures[1]

    breaktext(x) = strip.(split(x, "\n")) |> (x -> x[x .!= ""])

    (testname == ":") && (testname = ",")

    testrange = 1:10

    if occursin(",", testname)
        left, right = match(r"(.*),(.*)", testname).captures
        (left == "")  && (left = string(minimum(testrange)))
        (right == "") && (right = string(maximum(testrange)))

        for i in integer(left):integer(right)
            println("\n");
            println(Crayons.Box.YELLOW_FG("abstractlogic> t($i)"))
            abstractlogic("t($i)")
        end
        testcalls = String[]
    elseif testname ∈ ["1"]
        testcalls = ["clear;; a,b,c ∈ 1:4; a|c = 1; a > b; b > c"]
    elseif testname ∈ ["hp", "2", "harrypotter"]
        testcalls =  """
        clear;; a, b, c, d, e, f, g  ∈  NW, MA, MB, PO
        {{i}} == 'NW' {{2}} ; {{i}} == 'MA' {{1}}
        {{i}} == 'MB' {{1}} ; {{i}} == 'PO' {{3}}
        {{i!}} == 'NW' ==> {{i-1}} == 'PO'
        a != g; a,g != 'MA'; c,f != 'PO'; b == f
        """ |> breaktext
    elseif  testname ∈ ["lsatq11", "3"]
        testcalls = """
        a, b, c, d, e  ∈  0:4 || unique [clear]
        b == a*2; c < e
        prove: e == 1 ==> a==2
        range
        prove: {{i}} == 1 &&& {{>i}} == 2 {{1}}
        prove: c != 0 ==> b == 2
        any: c == 0 &&& d == 1
        show
        prove: {{i}} == 1|2 ==> {{i+2}} == 1|2
        prove: {{i}} == 1|3 ==> {{i+4}} != 1,3
        """ |> breaktext
    elseif  testname ∈ ["lsatq12", "4"]
        println("This problem has been simplified for testing purposes")
        testcalls = """
        t.1, t.2, t.3, f.1, f.2, f.3  ∈  _, Greed, Harvest, Limelight [clear]
        {{j}}.3 != '_'
        {{j}}.1 != '_' ==> {{j}}.2 != '_'
        {{i}} == 'Greed' {{1,}}
        {{i}} == 'Harvest' {{1,}}
        {{i}} == 'Limelight' {{1,}}
        {{j}}.1 != '_' ==> {{j}}.1 != {{j}}.2 &&& {{j}}.1 != {{j}}.3 &&& {{j}}.2 != {{j}}.3
        t.3 = 'Harvest'
        f.3 = 'Greed'|'Limelight'; f.2,f.1 != 'Greed','Limelight'
        export lsat3
        check t.2 = 'Limelight' ; t.3 = 'Harvest' ; f.3 = 'Limelight'
        any {{i}} == 'Limelight' {{3,}}
        any f.1,f.2,f.3 |= 'Greed'
        t.0, f.0 ∈ Greed, Harvest, Limelight
        """  |> breaktext
    elseif  testname ∈ ["lsatq13", "5"]
        testcalls = """
        w1, w2, w3, w4, w5, w6, w7  ∈ Guadeloupe, Jamaica, Martinique, Trinidad [clear]
        {{i}} = 'Guadeloupe' {{1,}}
        {{i}} = 'Jamaica' {{1,}}
        {{i}} = 'Martinique' {{1,}}
        {{i}} = 'Trinidad' {{1,}}
        w4 != 'Jamaica'
        w7 = 'Trinidad'
        {{i}} = 'Martinique' {{2}}
        {{<i}} = 'Martinique' &&& {{>i}} = 'Martinique' &&& {{i}} = 'Guadeloupe' {{1,}}
        {{i!}} = 'Jamaica' ==> {{i-1}} = 'Guadeloupe'
        {{i}} != {{i+1}}
        export lsat4
        check w1 = 'Guadeloupe'; w2 = 'Jamaica'; w3 = 'Martinique'; w4 = 'Trinidad'; w5 = 'Guadeloupe'; w6 = 'Martinique'; w7 = 'Trinidad'
        check w1 = 'Guadeloupe'; w2 = 'Martinique'; w3 = 'Trinidad'; w4 = 'Martinique'; w5 = 'Guadeloupe'; w6 = 'Jamaica'; w7 = 'Trinidad'
        check w1 = 'Jamaica'; w2 = 'Martinique'; w3 = 'Guadeloupe'; w4 = 'Martinique'; w5 = 'Guadeloupe'; w6 = 'Jamaica'; w7 = 'Trinidad'
        check w1 = 'Martinique'; w2 = 'Trinidad'; w3 = 'Guadeloupe'; w4 = 'Jamaica'; w5 = 'Martinique'; w6 = 'Guadeloupe'; w7 = 'Trinidad'
        check w1 = 'Martinique'; w2 = 'Trinidad'; w3 = 'Guadeloupe'; w4 = 'Trinidad'; w5 = 'Guadeloupe'; w6 = 'Jamaica'; w7 = 'Martinique'
        check w6 = 'Trinidad';; check w5 = 'Martinique'
        check w6 = 'Jamaica';; check w3 = 'Jamaica'
        check w3 = 'Guadeloupe'
        w5 = 'Trinidad'
        check w1 = 'Trinidad';; check w2 = 'Martinique'
        check w3 = 'Guadeloupe';; check w4 = 'Martinique'
        check w6 = 'Jamaica'
        import lsat4
        w1 = 'Guadeloupe'; w5 = 'Jamaica'
        check w2 = 'Jamaica';; check w2 = 'Trinidad'
        check w3 = 'Martinique';; check w6 = 'Guadeloupe'
        check w6 = 'Martinique'
        import lsat4
        w1 = 'Guadeloupe'; w2 = 'Trinidad'
        prove w3 = 'Martinique';; prove w4 = 'Martinique'
        prove w5 = 'Martinique';; prove w3 = 'Guadeloupe'
        prove w5 = 'Guadeloupe'
        import lsat4
        w3 = 'Martinique'
        any w4 = 'Guadeloupe'; w5 = 'Trinidad'
        any w4 = 'Jamaica'; w5 = 'Guadeloupe'
        any w4 = 'Martinique'; w5 = 'Trinidad'
        any w4 = 'Trinidad'; w5 = 'Jamaica'
        any w4 = 'Trinidad'; w5 = 'Martinique'
        import lsat4
        check w1|w2 = 'Guadeloupe'
        w1,w2 != 'Guadeloupe'
        check w2|w3 = 'Martinique'
        check {{i}} = 'Guadeloupe' {{,2}}
        check {{i}} = 'Jamaica' {{,2}}
        check {{i}} = 'Trinidad' {{,2}}
        """  |> breaktext
    elseif  testname ∈ ["lsatonlineq1", "6"]
        testcalls = """
        G, L, M, N, P, R, S, W ∈ red, not [clear]
        G, S = 'red' ==> W = 'red'
        N = 'red' ==> R, S = 'not'
        P = 'red' ==> L = 'not'
        L|M|R = 'red' {2}
        M, R = 'red'
        all: G, L = 'not'
        all: G, N = 'not'
        all: L, N = 'not'
        all: L, P = 'not'
        all: P, S = 'not'
        T, U, V, W, X, Y, Z ∈ unique [clear]
        X != 1,2
        W > X
        T,Y != 7
        Y|Z = W + 1
        V = U + 1 ||| V = U - 1
        """ |> breaktext
    elseif  testname ∈ ["lsatonlineq2", "7"]
        testcalls = """
        T, U, V, W, X, Y, Z ∈ unique [clear]
        X != 1,2
        W > X
        T,Y != 7
        Y|Z = W + 1
        V = U + 1 ||| V = U - 1
        export preserver
        V = 1
        all: T = 6
        all: X = 3
        all: Z = 7
        all: T = Y + 1
        all: W = X + 1
        import preserver
        U = 3
        range Y
        """ |> breaktext
    elseif  testname ∈ ["lsatonlineq3", "8"]
        println("This problem has been simplified for testing purposes")
        testcalls = """
        q1.m, q1.t, q1.w, q1.y ∈ 0:1 [clear]
        q2.m, q2.t, q2.w, q2.y ∈ 0:1
        q3.m, q3.t, q3.w, q3.y ∈ 0:1
        {{>=:q1.m,<=:q1.y}} == 1 {{1,}}
        {{>=:q2.m,<=:q2.y}} == 1 {{1,}}
        {{>=:q3.m,<=:q3.y}} == 1 {{1,}}
        {{j}}.m = 1 {{1,}}
        {{j}}.t = 1 {{1,}}
        {{j}}.w = 1 {{1,}}
        {{j}}.y = 1 {{1,}}
        {{i}} = 1 {{,6}}
        {{j}}.m = 1 ==> {{j+1}}.m != 1; {{j}}.t = 1 ==> {{j+1}}.t != 1
        {{j}}.w = 1 ==> {{j+1}}.w != 1; {{j}}.y = 1 ==> {{j+1}}.y != 1
        {{j}}.m = 1 {{2}}
        q2.w = 1
        export preserver
        {{j}}.w = 1 &&& {{j}}.y = 1 {{1,}}
        any: q2.m = 1
        J, K, L, M, N, P, Q ∈ 0:1 [clear]
        {{i}} = 1 {{4}}
        J ^ K
        N ^ P
        L! ==> N!
        K! ==> !Q
        export preserver
        P = 1
        show
        import preserver
        check: J == 1 &&& L == 1
        """ |> breaktext
    elseif  testname ∈ ["lsatonlineq4", "9"]
        testcalls = """
        J, K, L, M, N, P, Q ∈ 0:1 [clear]
        {{i}} = 1 {{4}}
        J ^ K
        N ^ P
        L! ==> N!
        K! ==> !Q
        export preserver
        P = 1
        show
        import preserver
        check: J == 1 &&& L == 1
        check: K == 1 &&& M == 1
        check: L == 1 &&& N == 1
        check: L == 1 &&& Q == 1
        check: M == 1 &&& Q == 1
        """ |> breaktext
    elseif  testname ∈ ["lsatonlineq5", "10"]
        testcalls = """
        s, h, j, k, l, m ∈ unique [clear]
        s = h + 1
        k < s
        m < l
        m = j + 1 ||| m = j - 1
        l < k
        range h
        """ |> breaktext
    else
        (testname != "") && replthrow("Test " * testname * " not found!")
        printmarkdown("Tests available: " * join(testrange, ", "))
        testcalls = String[]
    end

    for v in testcalls
        print(Crayons.Box.GREEN_FG("\nabstractlogic> "))
        println(v)
        abstractlogic(v)
        replerror && throw("Testcall Error")
    end
end

# function testcall(userinput; verbose=true)
#
#     testname = match(r"t(?:est)?\((.*)\)", userinput).captures[1]
#
#     folders = readdir()
#     ("examples" ∈ folders) && (testfiles = readdir("examples/repl/"))
#     !("examples" ∈ folders) && (testfiles = readdir("../examples/repl/"))
#
#     testfiles = testfiles .|> x -> replace(x, r".jl$"=>"") |> lowercase
#
#     occursin(r"^[0-9]{1,2}$", testname) && (testname = testfiles[parse(Int64, testname)])
#
#     (testname == "") && return printmarkdown("Tests available: `" * join(testfiles, "`, `") *
#         "`\n \n `abstractlogic> t(testname)`")
#
#     !(testname ∈ testfiles) &&
#       return replthrow("\"$testname\" not found. Tests available: " * join(testfiles, ", "))
#
#     ("examples" ∈ folders) && (testread = readlines("examples/repl/$testname.jl"))
#     !("examples" ∈ folders) && (testread = readlines("../examples/repl/$testname.jl"))
#
#     testread = (testread.|> x -> strip(x)) |> x -> x[x .!= ""] |> x -> x[(x .|> y -> y[1]) .!= '#']
#
#     codestart = findall(y -> occursin(r"^Start the repl", y), testread)
#
#     (length(codestart) == 0) && return replthrow("\"Start the repl\" not found!")
#     (length(codestart) > 1 ) && return replthrow("Multiple \"Start the repl\" found!")
#
#     testlines = testread[(codestart[1]+1):end]
#
#     for v in testlines
#         println(Crayons.Box.GREEN_FG("\nabstractlogic> " * replace(v, "abstractlogic> "=> "")))
#         abstractlogic(v)
#         replerror && throw("")
#     end
# end


text:
# Recursively list files under a directory.
#
# Args:
#   root: Path to descend.
#
# Returns:
#   A vector of paths relative to root.
#
function walkdir(root::String)
    root = abspath(root)
    contents = String[]
    stack = String[]
    push!(stack, root)
    while !isempty(stack)
        path = pop!(stack)
        for f in readdir(path)
            fullpath = joinpath(path, f)
            if isdir(fullpath)
                push!(stack, fullpath)
            else
                push!(contents, fullpath)
            end
        end
    end
    contents
end


text:
"""
`module  JAC.PairAnnihilation1Photon`  ... a submodel of JAC that contains all methods for computing positron-bound-electron pair annihilation
                                           (PEPA) with two-photon emission cross sections and rates; e^+ + |i(N)> --> |f(N-1)> + photon + photon;
                                           it is using JAC, JAC.ManyElectron.
"""
module PairAnnihilation2Photon

    using JAC, JAC.ManyElectron

    """
    `struct  PairAnnihilation2Photon.Settings`  ... defines a type for the details and parameters in computing positron-bound-electron pair
             annihilation (PEPA) with single-photon emission cross sections and rates; e^+ + |i(N)> --> |f(N-1)> + photon1 + photon2.

        + multipoles              ::Array{EmMultipole}           ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge}              ... Specifies the gauges to be included into the computations.
        + positronEnergies        ::Array{Float64,1}             ... List of positron energies.
        + printBeforeComputation  ::Bool                         ... True, if all energies and lines are printed before their evaluation.
        + selectLines             ::Bool                         ... True, if lines are selected individually for the computations.
        + selectedLines           ::Array{Tuple{Int64,Int64},1}  ... List of lines, given by tupels (inital-level, final-level).
    """
    struct Settings
        multipoles                ::Array{EmMultipole}
        gauges                    ::Array{UseGauge}
        positronEnergies          ::Array{Float64,1}
        printBeforeComputation    ::Bool
        selectLines               ::Bool
        selectedLines             ::Array{Tuple{Int64,Int64},1}
    end

    """
    `JAC.PairAnnihilation2Photon.Settings()`  ... constructor for the default values of pair-annihilation photon line computations
    """
    function Settings()
        Settings(EmMultipole[], UseGauge[], Float64[], false, false, Tuple{Int64,Int64}[])
    end

    # `Base.show(io::IO, settings::PairAnnihilation2Photon.Settings)`
    #		... prepares a proper printout of the variable settings::PairAnnihilation2Photon.Settings.
    function Base.show(io::IO, settings::PairAnnihilation2Photon.Settings)
        println(io, "multipoles:               $(settings.multipoles)  ")
        println(io, "gauges:                   $(settings.gauges)  ")
        println(io, "positronEnergies:         $(settings.positronEnergies)  ")
        println(io, "printBeforeComputation:   $(settings.printBeforeComputation)  ")
        println(io, "selectLines:              $(settings.selectLines)  ")
        println(io, "selectedLines:            $(settings.selectedLines)  ")
    end

    """
    `struct  PairAnnihilation2Photon.Channel`  ... defines a type for a positron-bound-electron pair annihilation (PEPA) with two-photon
                                                   emission channel that specifies all quantum numbers, phases and amplitudes.

        + multipole1     ::EmMultipole          ... Multipole of the emitted photon1.
        + multipole2     ::EmMultipole          ... Multipole of the emitted photon2.
        + gauge          ::EmGauge              ... Gauge for dealing with the (coupled) radiation field.
        + kappa          ::Int64                ... partial-wave of the incoming free positron
        + symmetry       ::LevelSymmetry        ... total angular momentum and parity of the scattering state
        + phase          ::Float64              ... phase of the partial wave
        + amplitude      ::Complex{Float64}     ... PairAnnihilation2PhotonChannel amplitude associated with the given channel.
   """
    struct  Channel
        multipole1       ::EmMultipole
        multipole2       ::EmMultipole
        gauge            ::EmGauge
        kappa            ::Int64
        symmetry         ::LevelSymmetry
        phase            ::Float64
        amplitude        ::Complex{Float64}
    end

    """
    `struct  PairAnnihilation2Photon.Line`  ... defines a type for a positron-bound-electron pair-annihilation (two-photon) line that may
                                                include the definition of channels.

        + initialLevel   ::Level                 ... initial-(state) level
        + finalLevel     ::Level                 ... final-(state) level
        + positronEnergy ::Float64                ... Energy of the (incoming free) positron.
        + photon1Energy  ::Float64                ... Energy of the emitted photon1.
        + photon2Energy  ::Float64                ... Energy of the emitted photon2.
        + crossSection   ::EmProperty             ... Cross section for this pair-annihilation (photon) line.
        + hasChannels    ::Bool                   ... Determines whether the individual (sub-) channels are defined in terms of their
                                                      free-positron energy, kappa, multipole, etc., or not.
        + channels       ::Array{PairAnnihilation2Photon.Channel,1}  ... List of PairAnnihilation2Photon.Channels of this line.
    """
    struct  Line
        initialLevel     ::Level
        finalLevel       ::Level
        positronEnergy   ::Float64
        photon1Energy    ::Float64
        photon2Energy    ::Float64
        crossSection     ::EmProperty
        hasChannels      ::Bool
        channels         ::Array{PairAnnihilation2Photon.Channel,1}
    end

    """
    `JAC.PairAnnihilation2Photon.Line()`  ... 'empty' constructor for a pair-annihilation (photon) line between a specified initial and
                                              final level.
    """
    function Line()
        Line(Level(), Level(), 0., 0., 0., EmProperty(0., 0.), false, PairAnnihilation2Photon[] )
    end

    # `Base.show(io::IO, line::PairAnnihilation2Photon.Line)`  ... prepares a proper printout of the variable line::PairAnnihilation2Photon.Line.
    function Base.show(io::IO, line::PairAnnihilation1Photon.Line)
        println(io, "initialLevel:      $(line.initialLevel)  ")
        println(io, "finalLevel:        $(line.finalLevel)  ")
        println(io, "positronEnergy:    $(line.positronEnergy)  ")
        println(io, "photon1Energy:     $(line.photon1Energy)  ")
        println(io, "photon2Energy:     $(line.photon2Energy)  ")
        println(io, "crossSection:      $(line.crossSection)  ")
        println(io, "hasChannels:       $(line.hasChannels)  ")
        println(io, "channels:          $(line.channels)  ")
    end

end # module


text:
# lid.jl
#
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

export lid_train, lid_save, lid_load, lid

# function to train the model using training data file
function lid_train(filename::String)
    train_data  = map(l -> split(chomp(l), '\t')[2], filelines(filename))
    train_truth = map(l -> split(chomp(l), '\t')[1], filelines(filename))
    bkgmodel, fextractor, model = tc_train(train_data, train_truth, lid_iterating_tokenizer, mincount = 2, cutoff = 1e10,
                                           trainer = (fvs, truth, init_model) -> train_mira(fvs, truth, init_model, iterations = 20, k = 2, C = 0.01, average = true),
                                           iteration_method = :eager)
    return bkgmodel, fextractor, model
end

# function to save the model to file
function lid_save(filename::String, bkgmodel::Text.BKG, model::LinearModel{String})
    jldopen(filename, "w", compress=true) do file
        g = g_create(file, "lid")
        g["bkgmodel_dict"]  = Dict(bkgmodel.dict)
        g["bkgmodel_index"] = bkgmodel.index
        g["bkgmodel_stats"] = bkgmodel.stats
        g["model"]          = model
    end
end

# function to load the model from file
function lid_load(filename::String)
    bkgmodel, model = None, None
    jldopen(filename, "r") do file
        g = file["lid"]
        bkgmodel_dict   = read(g["bkgmodel_dict"])
        bkgmodel_index  = read(g["bkgmodel_index"])
        bkgmodel_stats  = read(g["bkgmodel_stats"])
        bkgmodel    = Text.BKG(DefaultDict(1, bkgmodel_dict), bkgmodel_index, bkgmodel_stats)
        model       = read(g["model"])
    end
    return bkgmodel, ptext -> Text.tc_features(ptext, bkgmodel), model
end

# function to detect language of text using model
function detect_language(text::String, fextractor::Function, model::LinearModel{String})
    fv      = fextractor(lid_iterating_tokenizer(text))
    scores  = score(model, fv)
    bidx, b = best(scores)
    return model.index_class[bidx]
end

# function to detect language from HTTP Request
function lid(req::Request, result::Dict, fextractor::Function, model::LinearModel{String})
    if req.method == "PUT"
        text    = bytestring(req.data)
        lang    = detect_language(text, fextractor, model)
        result["text"]      = text
        result["language"]  = lang
    elseif req.method == "GET"
        result["all_languages"] = model.index_class
    else
        result["language"]  = "error"
    end
end

max_forks_repo_name: trevorlewis/TEXT-Language-REST


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Aria2_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Aria2")
JLLWrappers.@generate_main_file("Aria2", UUID("9ab3bdc3-1250-5043-8fac-ac7e82d2cbc9"))
end  # module Aria2_jll


text:
function singleModel!(data::alpheusData)
    options = data.options
    input = data.input
    preprocessor = data.preprocessor

    ### initialize model ###
    if options.solver == bonmin
        model = Model(() -> AmplNLWriter.Optimizer(Bonmin_jll.amplexe));
    elseif options.solver == couenne
        model = Model(() -> AmplNLWriter.Optimizer(Couenne_jll.amplexe));
    elseif options.solver == ipopt
        model = Model(Ipopt.Optimizer)
    end

    ### set parameters ###
    for (key, value) in options.parameters
        set_optimizer_attribute(model, key, value)
    end

    ### create variables ###
    @variable(model, 0.0 <= H[1:J] <= deepest, start = 1.0);
    @variable(model, 0.0 <= V[1:J] <= g*deepest, start = 1.0);
    @variable(model, z[1:J], Bin);

    ### create expressions ###
    @NLexpression(model, W[j in 1:J], z[j]*input.N_turbine*(V[j]/input.Vmax^2)^(3/2)*preprocessor.K);
    @NLexpression(model, A[j in 1:J], input.B*H[j]);
    @NLexpression(model, R[j in 1:J], A[j]/(input.B+2*H[j]));
    @NLexpression(model, S[j in 1:J], V[j]*(input.n0^2+input.nt^2*z[j]+2*input.n0*input.nt*z[j])/(R[j]^(4/3)));
    @NLexpression(model,SM[j in 1:J-1],(S[j] + S[j+1])/2);

    ### create constraint ###
    if !options.upstream
        @constraint(model, H[1] <= input.Hu);
    end
    if !options.downstream
        @constraint(model, H[J] >= input.Hd);
    end
    @constraint(model, sum(z) <= input.f);

    @constraint(model, [j in 1:J], V[j] >= input.Vmin*z[j]);
    @constraint(model, [j in 1:J], V[j] <= (input.Vmax - g*deepest)*z[j] + g*deepest);

    @constraint(model, z[1] == 0);
    @constraint(model, z[J] == 0);
    @constraint(model, [j in 1:J], H[j] >= input.Gt*z[j]);
    @constraint(model, [j in 1:J-1], -(1-z[j])*input.M*H[j]-z[j]*deepest <= H[j+1] - H[j]);
    @constraint(model, [j in 1:J-1], H[j+1] - H[j] <= (1-z[j])*input.M*H[j]+z[j]*deepest);
    @constraint(model, [j in 1:J], V[j] <= g*H[j]);
    @NLconstraint(model,[j in 1:J], V[j]*A[j]^2 == input.Q^2);
    @NLconstraint(model, [j in 1:J-1], V[j]/(2*g) + H[j] + preprocessor.E[j] == V[j+1]/(2*g) + H[j+1] + preprocessor.E[j+1] + SM[j]*preprocessor.Δx + W[j]);

    ### create objective ###

    @NLobjective(model, Max,sum(W[j] for j in 1:J) -(input.Hu-H[1])*input.ηu*options.upstream + (H[J]-input.Hd)*input.ηd*options.downstream)

    optimizeModel!(model)
    fillOutput!(model,data)
    nothing
end

function optimizeModel!(model::Model)
    @time status = optimize!(model)
    nothing
end

function fillOutput!(model::Model,data::alpheusData)
    output = data.output

    output.W = value.(model[:W]);
    output.H = value.(model[:H]);
    output.z = value.(model[:z]);
    V = value.(model[:V]);
    output.U = sqrt.(V)

    nothing
end

max_stars_repo_name: Thiago-NovaesB/alpheus

max_issues_repo_name: Thiago-NovaesB/alpheus

max_forks_repo_name: Thiago-NovaesB/alpheus


text:
using StatisticalRethinking, CmdStan
#gr(size=(600,600));

ProjDir = rel_path("..", "scripts", "03")
cd(ProjDir)

binomialstanmodel = "
// Inferring a Rate
data {
  int N;
  int<lower=0> k[N];
  int<lower=1> n[N];
}
parameters {
  real<lower=0,upper=1> theta;
  real<lower=0,upper=1> thetaprior;
}
model {
  // Prior Distribution for Rate Theta
  theta ~ beta(1, 1);
  thetaprior ~ beta(1, 1);

  // Observed Counts
  k ~ binomial(n, theta);
}
";

stanmodel = Stanmodel(name="binomial", monitors = ["theta"], model=binomialstanmodel,
  output_format=:mcmcchains);

N2 = 4
n2 = Int.(9 * ones(Int, N2))
k2 = [6, 5, 7, 6]

binomialdata = Dict("N" => length(n2), "n" => n2, "k" => k2);

rc, chn, cnames = stan(stanmodel, binomialdata, ProjDir, diagnostics=false,
  CmdStanDir=CMDSTAN_HOME);

MCMCChains.describe(chn)

MCMCChains.hpd(chn)

if rc == 0
  mixeddensity(chn, xlab="height [cm]", ylab="density")
  bnds = hpd(chn[:,1,1])
  vline!([bnds[:lower]], line=:dash)
  vline!([bnds[:upper]], line=:dash)
end

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
using CSTParser: UnarySyntaxOpCall, TimesOp

struct CCallAmpersand; end
register(CCallAmpersand, Deprecation(
    "The amersand (`&`) in ccall has been deprecated. Use Ref instead",
    "julia",
    v"0.7.0-DEV", v"0.7.0-DEV", typemax(VersionNumber)
))

match(CCallAmpersand, CSTParser.Call) do x
    # Match ccall only
    dep, expr, resolutions, context = x
    is_identifier(children(expr)[1], "ccall") || return
    has_cconv = any(x->is_identifier(children(expr)[5], x), ["cdecl", "stdcall", "fastcall", "thiscall", "llvmcall"])
    argtypes = children(expr)[7+(has_cconv ? 2 : 0)]
    expr′ = ChildReplacementNode(nothing, collect(children(expr)), expr)
    argtypes′ = children(expr′)[7+(has_cconv ? 2 : 0)] = ChildReplacementNode(expr′, collect(children(argtypes)), argtypes)
    found_any = false
    # Every other node is punctuation
    for (i, idx) in enumerate((9+(has_cconv ? 2 : 0)):2:length(children(expr)))
        node = children(expr)[idx]
        # Check if this is ampersand syntax
        if isexpr(node, UnarySyntaxOpCall) && isexpr(children(node)[1],
                OPERATOR, Tokens.AND)
            # Check if the argtype is `Ptr`.
            argT = children(argtypes)[2*i] # Every other node is punctuation
            isexpr(argT, Curly) || continue
            id = children(argT)[1]
            is_identifier(id, "Ptr") || continue
            # Rewrite to `Ref`
            children(argtypes′)[2*i] = replace_node(argT, id, ReplacementNode("Ref", id))
            # Remove ampersand
            children(expr′)[idx] = children(node)[2]
            found_any = true
        end
    end
    if found_any
        buf = IOBuffer()
        print_replacement(buf, expr′, false, false)
        push!(resolutions, TextReplacement(dep, expr.span, String(take!(buf))))
    end
end


text:
## temporary compatibility hack
if VERSION < v"0.4.0-dev"
    Base.Dict(z::Base.Zip2) = Dict(z.a, z.b)
    typealias AbstractTriangular Triangular
else
    import Base.LinAlg.AbstractTriangular
end

## inspired by DataFrames
if VERSION < v"0.5.0-dev+2023"
    displaysize(io::IO) = Base.tty_size()
else
    using Combinatorics
end

if VERSION < v"0.5.0-dev"
    String = ASCIIString
end


text:
# time-dot.jl
# some timing tests

# times reported below are all with julia 1.5
# on 2017 iMac (4.2GHz Intel Core i7) with Mojave 10.14.6

using ConstantArrays
#using FillArrays: Ones
using FillArrays: Ones, Trues # for 0.9
using BenchmarkTools

dim = (2^8,2^8,2^4)
x = randn(dim)
y = similar(x)

mask0 = trues(dim)
maskc = ConstantArray(dim)
maskf = Ones{Bool}(dim)
#const Trues = Ones{Bool, N, Axes} where {N, Axes} # lazy trues()
maskt = Trues(dim)
@assert mask0 == maskt

@assert mask0 == maskc == maskf == maskt

# function wrappers for timing tests

function dotter!(y, mask, x) # .*
	return y .= mask .* x
end

if true # .*
	@info "dotter"
	@assert dotter!(y, mask0, x) ==
		dotter!(y, maskc, x) ==
		dotter!(y, maskf, x) ==
		dotter!(y, maskt, x) == x
	@btime dotter!($y, $mask0, $x) # 1.0 ms
	@btime dotter!($y, $maskc, $x) # 0.6 ms :)
	@btime dotter!($y, $maskf, $x) # 0.6 ms :)
	@btime dotter!($y, $maskt, $x) # 0.6 ms :)
end

nothing


text:
is_dag(g) = is_directed(g) && !is_cyclic(g)
dag_or_error(g) = is_dag(g) || throw(DomainError(g, "Only Directed Acylic Graphs are supported"))

sources(graph) = findall(iszero, indegree(graph))
sinks(graph) = findall(iszero, outdegree(graph))

function longest_paths(graph, roots)
    dag_or_error(graph)
    dists = zeros(nv(graph))
    pending = [0 => r for r in roots]
    while(!isempty(pending))
        depth, node = pop!(pending)
        dists[node] = max(dists[node], depth)
        append!(pending, (depth+1) .=> outneighbors(graph, node))
    end
    return dists
end