text:
##########################################################
"""
$(TYPEDEF)

Abstract type for finite volume system structure.
"""
abstract type AbstractSystem{Tv<:Number, Ti <:Integer, Tm <:Integer} end

##################################################################
"""
$(TYPEDEF)

Structure holding data for finite volume system solution.

Information on species distribution is kept in sparse or dense matrices
matrices and, correspondingly,  the solution array is of type SparseSolutionArray
or matrix, respectively.

In the case of sparse unknown storage, the system matrix handles exactly those
degrees of freedom which correspond to unknowns. However, handling
of the sparse matrix structures for the bookkeeping of the unknowns
creates overhead.

$(TYPEDFIELDS)
"""
mutable struct System{Tv,Ti, Tm, TSpecMat<:AbstractMatrix, TSolArray<:AbstractMatrix} <: AbstractSystem{Tv,Ti, Tm}

    """
    Grid
    """
    grid::ExtendableGrid{Tv,Ti}

    """
    Physics data
    """
    physics::Physics

    """
    Array of boundary values
    """
    boundary_values::Array{Tv,2}

    """
    Array of boundary factors
    """
    boundary_factors::Array{Tv,2}

    """
    Matrix containing species numbers for inner regions
    """
    region_species::TSpecMat

    """
    Matrix containing species numbers for boundary regions
    """
    bregion_species::TSpecMat

    """
    Matrix containing degree of freedom numbers for each node
    """
    node_dof::TSpecMat

    """
    Jacobi matrix for nonlinear problem
    """
    matrix::ExtendableSparseMatrix{Tv,Tm}

    """
    Matrix factorization
    """
    factorization::Union{Nothing,ExtendableSparse.AbstractFactorization{Tv,Tm}}

    """
    Flag which says if the number of unknowns per node is constant
    """
    species_homogeneous::Bool

    """
    Number of quantities defined on system
    """
    num_quantities::Ti

    """
    Solution vector holding Newton update
    """
    update::TSolArray

    """
    Solution vector holding Newton residual
    """
    residual::TSolArray

    """
    Precomputed geometry factors for cell nodes
    """
    cellnodefactors::Array{Tv,2}

    """
    Precomputed geometry factors for cell edges
    """
    celledgefactors::Array{Tv,2}

    """
    Precomputed geometry factors for boundary nodes
    """
    bfacenodefactors::Array{Tv,2}

    """
    Precomputed geometry factors for boundary edges
    """
    bfaceedgefactors::Array{Tv,2}

    """
    Sparse matrix for generic operator handling
    """
    generic_matrix::SparseMatrixCSC

    """
    Sparse matrix colors for generic operator handling
    """
    generic_matrix_colors::Vector

    """
    Hash value of latest unknowns vector the assembly was called with
    """
    uhash::UInt64

    """
    Data for allocation check
    """
    allocs::Int

    System{Tv,Ti,Tm, TSpecMat, TSolArray}() where {Tv,Ti,Tm, TSpecMat, TSolArray} = new()
end

"""
    const DenseSystem

Type alias for system with dense matrix based species management
"""
const DenseSystem = System{Tv,Ti,Tm,Matrix{Ti},Matrix{Tv}} where {Tv, Ti, Tm}

"""
    const SparseSystem

Type alias for system with sparse matrix based species management
"""
const SparseSystem = System{Tv,Ti,Tm,SparseMatrixCSC{Ti,Ti},SparseSolutionArray{Tv,Ti} } where {Tv, Ti, Tm}

default_check_allocs()= haskey(ENV,"VORONOIFVM_CHECK_ALLOCS") ? parse(Bool,ENV["VORONOIFVM_CHECK_ALLOCS"]) :  false

##################################################################
"""
````
function System(grid,physics;
                unknown_storage=:dense;
                matrixindextype=Int32)
````

Create Finite Volume System.

- `grid`: 1D/2D/3D discretization grid
- `physics`: Physics struct containing node and edge callbacks
- `unknown_storage`: string or symbol:
     - `:dense` :  solution vector is an  `nspecies` x `nnodes`  dense matrix
     - `:sparse` :  solution vector is an `nspecies` x `nnodes`  sparse matrix
- `matrixindextype` : Index type for sparse matrices created in the system
"""
function System(grid,physics::Physics; unknown_storage=:dense, matrixindextype=Int64, check_allocs=default_check_allocs())
    system=System(grid,unknown_storage=unknown_storage, matrixindextype=matrixindextype, check_allocs=check_allocs)
    physics!(system,physics)
end

function physics!(system::System,physics)
    system.physics=physics
    system
end

function System(grid;unknown_storage=:dense, matrixindextype=Int64, check_allocs=default_check_allocs())
    Tv=coord_type(grid)
    Ti=index_type(grid)
    Tm=matrixindextype

    if Symbol(unknown_storage)==:dense
        system=System{Tv,Ti,Tm,Matrix{Ti}, Matrix{Tv}}()
    elseif Symbol(unknown_storage)==:sparse
        system=System{Tv,Ti,Tm,SparseMatrixCSC{Ti,Ti},SparseSolutionArray{Tv,Ti}}()
    else
        throw("specify either unknown_storage=:dense  or unknown_storage=:sparse")
    end

    maxspec=0
    system.grid=grid
    system.region_species=spzeros(Ti,Int16,maxspec,num_cellregions(grid))
    system.bregion_species=spzeros(Ti,Int16,maxspec,num_bfaceregions(grid))
    system.node_dof=spzeros(Ti,Tm,maxspec,num_nodes(grid))
    system.boundary_values=zeros(Tv,maxspec,num_bfaceregions(grid))
    system.boundary_factors=zeros(Tv,maxspec,num_bfaceregions(grid))
    system.species_homogeneous=false
    system.num_quantities=0
    system.uhash=0x0
    system.allocs=-1000
    system.factorization=nothing
    check_allocs!(system,check_allocs)
    return system
end

"""
$(SIGNATURES)

Constructor for DenseSystem.
"""
DenseSystem(grid,physics::Physics; matrixindextype=Int64)=System(grid,physics,matrixindextype=matrixindextype,unknown_storage=:dense)

"""
$(SIGNATURES)

Constructor for SparseSystem.
"""
SparseSystem(grid,physics::Physics; matrixindextype=Int64)=System(grid,physics,matrixindextype=matrixindextype,unknown_storage=:sparse)

###################################################################################################
# Test if allocated is zero and if allocation check is enabled
# However we need to be aware that @allocated reports allocations
# during the compilation phase. So we need to wait for at least
# one run of the system in order enact the checking.
function _check_allocs(system, allocated)
       if system.allocs>=0 # we had a couple of runs before to bridge the compilation phase
           system.allocs=allocated
           return system.allocs==0
       elseif system.allocs > -100 # probably still in compiling phase
           system.allocs=system.allocs+1
           return true
       else
           # otherwise, checking has been switched off.
           return true
       end
end

"""
```
check_allocs!(system,true_or_false)
```

Enable/disable checking for time-consuming allocations in the assembly loop.
By default,  this check  is switched off.  By setting  the environment
variable `ENV["VORONOIFVM_CHECK_ALLOCS"]="true"`, this  default can be
changed.

Unless  the   matrix  pattern  changes,  there   shouldn't  occur  any
allocations in this loop. The check  method is aware of matrix pattern
changes. As a consequence, allocations in the assembly loop are mostly
due to type instabilities in physics callbacks, see the the discussion
[here](../runexamples/#Performance-with-closures).  Type instabilities
can be debugged via the `@time`  macro applied to expressions in a
physics callback.

The following  cases provide some ideas  where to look for  reasons of
the problem and possible remedies:

Case 1: a parameter changes its value, and Julia is not sure about the type.
```julia
eps=1.0

flux(f,_u,edge)
    u=unkowns(edge,_u)
    f[1]=eps*(u[1,1]-[1,2])
end
... solve etc ...
eps=2.0
```
Remedy: use a type annotation `eps::Float64=...` to signalize your intent to Julia.
This behaviour is explained in the [Julia documentation](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured).

Case 2: variables in the closure have the same name as a variable
introduced in a callback.
```julia
flux(f,_u,edge)
    u=unkowns(edge,_u)
    f[1]=(u[1,1]-[1,2])
end

... create etc ...

u=solve(...)
```
Remedy: rename e.g. `u=solve()` to `sol=solve()`

"""
function check_allocs!(system::AbstractSystem,chk::Bool)
    if chk
        system.allocs=-2
    else
        system.allocs=-1000
    end
    system
end

##################################################################

# Constant to be used as boundary condition factor
# to mark Dirichlet boundary conditons.
const Dirichlet=1.0e30

#################################################################

"""
    addzrows(matrix,maxrow)

Return matrix with number of rows increased to maxrow, and set
the new elements to zero.
"""
function addzrows(matrix::Matrix,maxrow)
    nrow,ncol=size(matrix)
    if maxrow<=nrow
        return matrix
    end
    newmatrix=zeros(eltype(matrix),maxrow,ncol)
    for icol=1:ncol
        for irow=1:nrow
            newmatrix[irow,icol]=matrix[irow,icol]
        end
    end
    newmatrix
end

function addzrows(matrix::SparseMatrixCSC,maxrow)
    nrow,ncol=size(matrix)
    if maxrow<=nrow
        return matrix
    end
    SparseMatrixCSC(maxrow, matrix.n,matrix.colptr,matrix.rowval,matrix.nzval)
end

"""
     increase_num_species!(system,maxspec)

Increase number of species in system to maxspec by adding new rows to all  relevant
matrices.
"""
function increase_num_species!(system,maxspec)
    if maxspec<=num_species(system)
        return
    end

    if isdefined(system,:matrix)
        error("Unable to increase number of species to $(maxspec).\nPlease add species before first solver run.")
    end

    system.region_species=addzrows(system.region_species,maxspec)
    system.bregion_species=addzrows(system.bregion_species,maxspec)
    system.node_dof=addzrows(system.node_dof,maxspec)
    system.boundary_values=addzrows(system.boundary_values,maxspec)
    system.boundary_factors=addzrows(system.boundary_factors,maxspec)
end

##################################################################
"""
````
is_boundary_species(AbstractSystem, ispec) -> Bool
````

Check if species number corresponds to a boundary species.
"""
function is_boundary_species(system::AbstractSystem, ispec)
    isbspec=false
    for ibreg=1:num_bfaceregions(system.grid)
        if system.bregion_species[ispec,ibreg]>0
            isbspec=true
        end
    end
    return isbspec
end

##################################################################
"""
````
is_bulk_species(AbstractSystem, ispec) -> Bool
````

Check if species number corresponds to a bulk species.
"""
function is_bulk_species(system::AbstractSystem, ispec)
    isrspec=false
    for ixreg=1:num_cellregions(system.grid)
        if system.region_species[ispec,ixreg]>0
            isrspec=true
        end
    end
    return isrspec
end

##################################################################
"""
````
enable_species!(system,ispec,regions)
````

Add species `ispec` to a list of bulk regions. Species numbers for
bulk and boundary species have to be distinct.
"""
function enable_species!(system::AbstractSystem,ispec::Integer, regions::AbstractVector)
    increase_num_species!(system,ispec)

    if is_boundary_species(system,ispec)
        throw(DomainError(ispec,"Species is already boundary species"))
    end
    _cellregions=cellregions(system.grid)
    _cellnodes=cellnodes(system.grid)
    for i in eachindex(regions)
        ireg=regions[i]
        system.region_species[ispec,ireg]=ispec
        for icell=1:num_cells(system.grid)
            if _cellregions[icell]==ireg
                for iloc=1:num_targets(_cellnodes,icell)
                    iglob=_cellnodes[iloc,icell]
                    system.node_dof[ispec,iglob]=ispec
                end
            end
        end
    end
end

##################################################################
"""
````
enable_boundary_species!(system,ispec,regions)
````

Add species `ispec` to a list of boundary regions. Species numbers for
bulk and boundary species have to be distinct.

"""
function enable_boundary_species!(system::AbstractSystem, ispec::Integer, bregions::AbstractVector)
    increase_num_species!(system,ispec)

    if is_bulk_species(system,ispec)
        throw(DomainError(ispec,"Species is already bulk species"))
    end
    _bfaceregions=bfaceregions(system.grid)
    _bfacenodes=bfacenodes(system.grid)

    for i in eachindex(bregions)
        ireg=bregions[i]
        system.bregion_species[ispec,ireg]=ispec
        for ibface=1:num_bfaces(system.grid)
            if _bfaceregions[ibface]==ireg
                for iloc=1:size(_bfacenodes,1)
                    iglob=_bfacenodes[iloc,ibface]
                    system.node_dof[ispec,iglob]=ispec
                end
            end
        end
    end
end

# Create matrix in system and figure out if species
# distribution is homgeneous
function _complete!(system::AbstractSystem{Tv,Ti, Tm};create_newtonvectors=false) where {Tv,Ti, Tm}

    if isdefined(system,:matrix)
        return
    end
    system.matrix=ExtendableSparseMatrix{Tv,Tm}(num_dof(system), num_dof(system))
    system.species_homogeneous=true
    species_added=false
    for inode=1:size(system.node_dof,2)
        for ispec=1:size(system.node_dof,1)
            if system.node_dof[ispec,inode]==ispec
                species_added=true
            else
                system.species_homogeneous=false
            end
        end
    end

    if (!species_added)
        error("No species enabled.\n Call enable_species(system,species_number, list_of_regions) at least once.")
    end
    if create_newtonvectors
        system.residual=unknowns(system)
        system.update=unknowns(system)
    end
    update_grid!(system)
    if has_generic_operator(system)
        if has_generic_operator_sparsity(system)
            system.generic_matrix=system.physics.generic_operator_sparsity(system)
        else
            generic_operator(f,u)=system.physics.generic_operator(f,u,system)
            input=rand(num_dof(system))
            output=similar(input)
            tdetect=@elapsed begin
                sparsity_pattern = jacobian_sparsity(generic_operator,output,input)
                system.generic_matrix = Float64.(sparse(sparsity_pattern))
            end
            println("sparsity detection for generic operator: $(tdetect) s")
            if nnz(system.generic_matrix)==0
                error("Sparsity detection failed: no pattern found")
            end
        end
        tdetect=@elapsed begin
            system.generic_matrix_colors = matrix_colors(system.generic_matrix)
        end
        println("matrix coloring for generic operator: $(tdetect) s")
    end
end

"""
$(SIGNATURES)
Set generic operator sparsity, in the case where a generic operator has been
defined in physics.
"""
function generic_operator_sparsity!(system::AbstractSystem, sparsematrix::SparseMatrixCSC)
    system.generic_matrix=sparsematrix
end

"""
````
update_grid!(system; grid=system.grid)
````

Update grid (e.g. after rescaling of coordinates).
"""
function update_grid!(system::AbstractSystem{Tv,Ti,Tm}; grid=system.grid) where {Tv,Ti,Tm}

    geom        = grid[CellGeometries][1]
    csys        = grid[CoordinateSystem]
    coord       = grid[Coordinates]
    cellnodes   = grid[CellNodes]
    cellregions = grid[CellRegions]
    bgeom       = grid[BFaceGeometries][1]
    bfacenodes  = grid[BFaceNodes]
    nbfaces     = num_bfaces(grid)
    ncells      = num_cells(grid)

    system.cellnodefactors  = zeros(Tv, num_nodes(geom), ncells)
    system.celledgefactors  = zeros(Tv, num_edges(geom), ncells)
    system.bfacenodefactors = zeros(Tv, num_nodes(bgeom), nbfaces)
    system.bfaceedgefactors = zeros(Tv, num_edges(bgeom), nbfaces)

    for icell=1:ncells
        @views cellfactors!(geom, csys, coord, cellnodes, icell,
                            system.cellnodefactors[:,icell],system.celledgefactors[:,icell])
    end

    for ibface=1:nbfaces
        @views bfacefactors!(bgeom,csys,coord,bfacenodes,ibface,
                             system.bfacenodefactors[:,ibface], system.bfaceedgefactors[:, ibface])
    end
end

##################################################################
"""
$(SIGNATURES)

Check if degree of freedom is defined.
"""
isdof(system::AbstractSystem,ispec,inode)= system.node_dof[ispec,inode]==ispec

##################################################################
"""
$(SIGNATURES)

Number of species in system
"""
num_species(system::AbstractSystem) = size(system.node_dof,1)

##################################################################
"""
$(SIGNATURES)

Retrieve user data record.
"""
data(system::AbstractSystem) = system.physics.data

##################################################################
"""
$(SIGNATURES)

Set Dirichlet boundary condition for species ispec at boundary ibc:

``u_{ispec}=v`` on ``\\Gamma_{ibc}``
"""
function boundary_dirichlet!(system::AbstractSystem, ispec::Integer, ibc, v)
    increase_num_species!(system,ispec)
    system.boundary_factors[ispec,ibc]=Dirichlet
    system.boundary_values[ispec,ibc]=v
end

##################################################################
"""
$(SIGNATURES)

Set Neumann boundary condition for species ispec at boundary ibc:

``\\mathrm{flux}_{ispec}\\cdot \\vec n=v`` on ``\\Gamma_{ibc}``
"""
function boundary_neumann!(system::AbstractSystem, ispec, ibc, v)
    increase_num_species!(system,ispec)
    system.boundary_factors[ispec,ibc]=0.0
    system.boundary_values[ispec,ibc]=v
end

##################################################################
"""
$(SIGNATURES)

Set Robin boundary condition for species ispec at boundary ibc:

``\\mathrm{flux}_{ispec}\\cdot \\vec n + \\alpha u_{ispec}=v`` on ``\\Gamma_{ibc}``

"""
function boundary_robin!(system::AbstractSystem, ispec, ibc, α, v)
    increase_num_species!(system,ispec)
    system.boundary_factors[ispec,ibc]=α
    system.boundary_values[ispec,ibc]=v
end

##################################################################
"""
$(SIGNATURES)

Number of species (size of first dimension) of solution array.
"""
num_species(a::AbstractArray)=size(a,1)

#
# Initialize Dirichlet BC
#
function _initialize_dirichlet!(U::AbstractMatrix,system::AbstractSystem)
    _bfaceregions=bfaceregions(system.grid)
    _bfacenodes=bfacenodes(system.grid)
    for ibface=1:num_bfaces(system.grid)
        ibreg=_bfaceregions[ibface]
        for ispec=1:num_species(system)
            if system.boundary_factors[ispec,ibreg]≈ Dirichlet
                for inode=1:dim_grid(system.grid)
                    U[ispec,_bfacenodes[inode,ibface]]=system.boundary_values[ispec,ibreg]
                end
            end
        end
    end
end

function _initialize!(U::AbstractMatrix,system::AbstractSystem)
    _initialize_dirichlet!(U,system)
    _initialize_inactive_dof!(U,system)
end

function _eval_and_assemble_inactive_species(system::AbstractSystem,U,Uold,F) end

function _eval_and_assemble_inactive_species(system::DenseSystem,U,Uold,F)
    if system.species_homogeneous
        return
    end
    for inode=1:size(system.node_dof,2)
        for ispec=1:size(system.node_dof,1)
            if !isdof(system,ispec,inode)
                F[ispec,inode]+= U[ispec,inode]-Uold[ispec,inode];
                idof=dof(F,ispec,inode)
                system.matrix[idof,idof]+=1.0
            end
        end
    end
end

function _initialize_inactive_dof!(U::AbstractMatrix,system::AbstractSystem) end

function _initialize_inactive_dof!(U::DenseSolutionArray,system::DenseSystem)
    if system.species_homogeneous
        return
    end
    for inode=1:size(system.node_dof,2)
        for ispec=1:size(system.node_dof,1)
            if !isdof(system,ispec,inode)
                U[ispec,inode]=0
            end
        end
    end
end

function Base.show(io::IO,sys::AbstractSystem)
    str=@sprintf("%s(num_species=%d)",typeof(sys),num_species(sys))
    println(io,str)
end

#####################################################
has_generic_operator(sys::AbstractSystem) = sys.physics.generic_operator!=nofunc_generic
has_generic_operator_sparsity(sys::AbstractSystem) =  sys.physics.generic_operator_sparsity!=nofunc_generic_sparsity

##################################################################
"""
$(SIGNATURES)

Number of degrees of freedom for system.
"""
function num_dof(system::AbstractSystem) end

num_dof(system::SparseSystem)= nnz(system.node_dof)

num_dof(system::DenseSystem)= length(system.node_dof)

"""
$(SIGNATURES)

Create a solution vector for system.
If inival is not specified, the entries of the returned vector are undefined.
"""
function unknowns(system::AbstractSystem;inival=undef) end

unknowns(sys::SparseSystem{Tv,Ti,Tm};inival=undef) where {Tv,Ti, Tm}=unknowns(Tv,sys,inival=inival)

unknowns(system::DenseSystem{Tv,Ti,Tm};inival=undef) where {Tv,Ti, Tm} = unknowns(Tv,system,inival=inival)

"""
$(SIGNATURES)

Create a solution vector for system with elements of type `Tu`.
If inival is not specified, the entries of the returned vector are undefined.
"""
function unknowns(Tu::Type, sys::AbstractSystem; inival=undef) end

function unknowns(Tu::Type, system::SparseSystem;inival=undef)
    a0=Array{Tu}(undef,num_dof(system))
    if inival!=undef
        fill!(a0,inival)
    end
    return SparseSolutionArray(SparseMatrixCSC(system.node_dof.m,
                                               system.node_dof.n,
                                               system.node_dof.colptr,
                                               system.node_dof.rowval,
                                               a0
                                               )
    )
end

function unknowns(Tu::Type, sys::DenseSystem; inival=undef)
    a=Array{Tu}(undef,size(sys.node_dof)...)
    if inival!=undef
        fill!(a,inival)
    end
    return a
end

"""
$(SIGNATURES)

Reshape vector to fit as solution to system.
"""
function Base.reshape(v,system::AbstractSystem) end

Base.reshape(v::DenseSolutionArray,system::DenseSystem)=v

Base.reshape(v::SparseSolutionArray,sys::SparseSystem)=v

function Base.reshape(v::AbstractVector, sys::DenseSystem)
    @assert  length(v)==num_dof(sys)
    nspec=num_species(sys)
    reshape(v,Int64(nspec),Int64(length(v)/nspec))
end

function Base.reshape(v::AbstractVector,system::SparseSystem)
    @assert  length(v)==num_dof(system)
    SparseSolutionArray(SparseMatrixCSC(system.node_dof.m,
                                        system.node_dof.n,
                                        system.node_dof.colptr,
                                        system.node_dof.rowval,
                                        Vector(v)
                                        )
                        )
end

"""
$(SIGNATURES)

Calculate norm, paying attention to species distribution over regions
"""
function LinearAlgebra.norm(system::AbstractSystem,u,p) end

function LinearAlgebra.norm(system::DenseSystem,u,p)
    _initialize_inactive_dof!(u,system)
    norm(u,p)
end

LinearAlgebra.norm(system::SparseSystem,u,p)=norm(u.node_dof.nzval,p)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# LQ Factorizations

struct LQ{T,S<:AbstractMatrix} <: Factorization{T}
    factors::S
    τ::Vector{T}
    LQ{T,S}(factors::AbstractMatrix{T}, τ::Vector{T}) where {T,S<:AbstractMatrix} = new(factors, τ)
end

struct LQPackedQ{T,S<:AbstractMatrix} <: AbstractMatrix{T}
    factors::Matrix{T}
    τ::Vector{T}
    LQPackedQ{T,S}(factors::AbstractMatrix{T}, τ::Vector{T}) where {T,S<:AbstractMatrix} = new(factors, τ)
end

LQ(factors::AbstractMatrix{T}, τ::Vector{T}) where {T} = LQ{T,typeof(factors)}(factors, τ)
LQPackedQ(factors::AbstractMatrix{T}, τ::Vector{T}) where {T} = LQPackedQ{T,typeof(factors)}(factors, τ)

"""
    lqfact!(A) -> LQ

Compute the LQ factorization of `A`, using the input
matrix as a workspace. See also [`lq`](@ref).
"""
lqfact!(A::StridedMatrix{<:BlasFloat}) = LQ(LAPACK.gelqf!(A)...)
"""
    lqfact(A) -> LQ

Compute the LQ factorization of `A`. See also [`lq`](@ref).
"""
lqfact(A::StridedMatrix{<:BlasFloat})  = lqfact!(copy(A))
lqfact(x::Number) = lqfact(fill(x,1,1))

"""
    lq(A; full = false) -> L, Q

Perform an LQ factorization of `A` such that `A = L*Q`. The default (`full = false`)
computes a factorization with possibly-rectangular `L` and `Q`, commonly the "thin"
factorization. The LQ factorization is the QR factorization of `Transpose(A)`. If the explicit,
full/square form of `Q` is requested via `full = true`, `L` is not extended with zeros.

!!! note
    While in QR factorization the "thin" factorization is so named due to yielding
    either a square or "tall"/"thin" rectangular factor `Q`, in LQ factorization the
    "thin" factorization somewhat confusingly produces either a square or "short"/"wide"
    rectangular factor `Q`. "Thin" factorizations more broadly are also
    referred to as "reduced" factorizatons.
"""
function lq(A::Union{Number,AbstractMatrix}; full::Bool = false, thin::Union{Bool,Nothing} = nothing)
    # DEPRECATION TODO: remove deprecated thin argument and associated logic after 0.7
    if thin != nothing
        Base.depwarn(string("the `thin` keyword argument in `lq(A; thin = $(thin))` has ",
            "been deprecated in favor of `full`, which has the opposite meaning, ",
            "e.g. `lq(A; full = $(!thin))`."), :lq)
        full::Bool = !thin
    end
    F = lqfact(A)
    L, Q = F.L, F.Q
    return L, !full ? Array(Q) : mul!(Q, Matrix{eltype(Q)}(I, size(Q.factors, 2), size(Q.factors, 2)))
end

copy(A::LQ) = LQ(copy(A.factors), copy(A.τ))

LQ{T}(A::LQ) where {T} = LQ(convert(AbstractMatrix{T}, A.factors), convert(Vector{T}, A.τ))
Factorization{T}(A::LQ{T}) where {T} = A
Factorization{T}(A::LQ) where {T} = LQ{T}(A)
AbstractMatrix(A::LQ) = A.L*A.Q
AbstractArray(A::LQ) = AbstractMatrix(A)
Matrix(A::LQ) = Array(AbstractArray(A))
Array(A::LQ) = Matrix(A)

adjoint(A::LQ{T}) where {T} = QR{T,typeof(A.factors)}(adjoint(A.factors), A.τ)

function getproperty(F::LQ, d::Symbol)
    m, n = size(F)
    if d == :L
        return tril!(getfield(F, :factors)[1:m, 1:min(m,n)])
    elseif d == :Q
        return LQPackedQ(getfield(F, :factors), getfield(F, :τ))
    else
        return getfield(F, d)
    end
end

getindex(A::LQPackedQ, i::Integer, j::Integer) =
    mul!(A, setindex!(zeros(eltype(A), size(A, 2)), 1, j))[i]

function show(io::IO, C::LQ)
    println(io, "$(typeof(C)) with factors L and Q:")
    show(io, C.L)
    println(io)
    show(io, C.Q)
end

LQPackedQ{T}(Q::LQPackedQ) where {T} = LQPackedQ(convert(AbstractMatrix{T}, Q.factors), convert(Vector{T}, Q.τ))
AbstractMatrix{T}(Q::LQPackedQ) where {T} = LQPackedQ{T}(Q)
Matrix(A::LQPackedQ) = LAPACK.orglq!(copy(A.factors),A.τ)
Array(A::LQPackedQ) = Matrix(A)

size(F::LQ, dim::Integer) = size(getfield(F, :factors), dim)
size(F::LQ)               = size(getfield(F, :factors))

# size(Q::LQPackedQ) yields the shape of Q's square form
function size(Q::LQPackedQ)
    n = size(Q.factors, 2)
    return n, n
end
function size(Q::LQPackedQ, dim::Integer)
    if dim < 1
        throw(BoundsError())
    elseif dim <= 2 # && 1 <= dim
        return size(Q.factors, 2)
    else # 2 < dim
        return 1
    end
end

## Multiplication by LQ
mul!(A::LQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} =
    A.L * LAPACK.ormlq!('L', 'N', A.factors, A.τ, B)
mul!(A::LQ{T}, B::QR{T}) where {T<:BlasFloat} =
    A.L * LAPACK.ormlq!('L', 'N', A.factors, A.τ, Matrix(B))
mul!(A::QR{T}, B::LQ{T}) where {T<:BlasFloat} =
    mul!(zeros(eltype(A), size(A)), Matrix(A), Matrix(B))
function *(A::LQ{TA}, B::StridedVecOrMat{TB}) where {TA,TB}
    TAB = promote_type(TA, TB)
    mul!(Factorization{TAB}(A), copy_oftype(B, TAB))
end
function *(A::LQ{TA},B::QR{TB}) where {TA,TB}
    TAB = promote_type(TA, TB)
    mul!(Factorization{TAB}(A), Factorization{TAB}(B))
end
function *(A::QR{TA},B::LQ{TB}) where {TA,TB}
    TAB = promote_type(TA, TB)
    mul!(Factorization{TAB}(A), Factorization{TAB}(B))
end
*(A::Adjoint{<:Any,<:LQ}, B::LQ) = adjoint(A.parent) * B
*(A::LQ, B::Adjoint{<:Any,<:LQ}) = A * adjoint(B.parent)

## Multiplication by Q
### QB
mul!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = LAPACK.ormlq!('L','N',A.factors,A.τ,B)
function (*)(A::LQPackedQ, B::StridedVecOrMat)
    TAB = promote_type(eltype(A), eltype(B))
    mul!(AbstractMatrix{TAB}(A), copy_oftype(B, TAB))
end

### QcB
mul!(adjA::Adjoint{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasReal} =
    (A = adjA.parent; LAPACK.ormlq!('L','T',A.factors,A.τ,B))
mul!(adjA::Adjoint{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} =
    (A = adjA.parent; LAPACK.ormlq!('L','C',A.factors,A.τ,B))
function *(adjA::Adjoint{<:Any,<:LQPackedQ}, B::StridedVecOrMat)
    A = adjA.parent
    TAB = promote_type(eltype(A), eltype(B))
    if size(B,1) == size(A.factors,2)
        mul!(Adjoint(AbstractMatrix{TAB}(A)), copy_oftype(B, TAB))
    elseif size(B,1) == size(A.factors,1)
        mul!(Adjoint(AbstractMatrix{TAB}(A)), [B; zeros(TAB, size(A.factors, 2) - size(A.factors, 1), size(B, 2))])
    else
        throw(DimensionMismatch("first dimension of B, $(size(B,1)), must equal one of the dimensions of A, $(size(A))"))
    end
end

### QBc/QcBc
function *(A::LQPackedQ, adjB::Adjoint{<:Any,<:StridedVecOrMat})
    B = adjB.parent
    TAB = promote_type(eltype(A), eltype(B))
    BB = similar(B, TAB, (size(B, 2), size(B, 1)))
    adjoint!(BB, B)
    return mul!(A, BB)
end
function *(adjA::Adjoint{<:Any,<:LQPackedQ}, adjB::Adjoint{<:Any,<:StridedVecOrMat})
    A, B = adjA.parent, adjB.parent
    TAB = promote_type(eltype(A), eltype(B))
    BB = similar(B, TAB, (size(B, 2), size(B, 1)))
    adjoint!(BB, B)
    return mul!(Adjoint(A), BB)
end

# in-place right-application of LQPackedQs
# these methods require that the applied-to matrix's (A's) number of columns
# match the number of columns (nQ) of the LQPackedQ (Q) (necessary for in-place
# operation, and the underlying LAPACK routine (ormlq) treats the implicit Q
# as its (nQ-by-nQ) square form)
mul!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasFloat} =
    LAPACK.ormlq!('R', 'N', B.factors, B.τ, A)
mul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasReal} =
    (B = adjB.parent; LAPACK.ormlq!('R', 'T', B.factors, B.τ, A))
mul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasComplex} =
    (B = adjB.parent; LAPACK.ormlq!('R', 'C', B.factors, B.τ, A))

# out-of-place right application of LQPackedQs
#
# LQPackedQ's out-of-place multiplication behavior is context dependent. specifically,
# if the inner dimension in the multiplication is the LQPackedQ's second dimension,
# the LQPackedQ behaves like its square form. if the inner dimension in the
# multiplication is the LQPackedQ's first dimension, the LQPackedQ behaves like either
# its square form or its truncated form depending on the shape of the other object
# involved in the multiplication. we treat these cases separately.
#
# (1) the inner dimension in the multiplication is the LQPackedQ's second dimension.
# in this case, the LQPackedQ behaves like its square form.
#
function *(A::StridedVecOrMat, adjQ::Adjoint{<:Any,<:LQPackedQ})
    Q = adjQ.parent
    TR = promote_type(eltype(A), eltype(Q))
    return mul!(copy_oftype(A, TR), Adjoint(AbstractMatrix{TR}(Q)))
end
function *(adjA::Adjoint{<:Any,<:StridedMatrix}, adjQ::Adjoint{<:Any,<:LQPackedQ})
    A, Q = adjA.parent, adjQ.parent
    TR = promote_type(eltype(A), eltype(Q))
    C = adjoint!(similar(A, TR, reverse(size(A))), A)
    return mul!(C, Adjoint(AbstractMatrix{TR}(Q)))
end
#
# (2) the inner dimension in the multiplication is the LQPackedQ's first dimension.
# in this case, the LQPackedQ behaves like either its square form or its
# truncated form depending on the shape of the other object in the multiplication.
#
# these methods: (1) check whether the applied-to matrix's (A's) appropriate dimension
# (columns for A_*, rows for Ac_*) matches the number of columns (nQ) of the LQPackedQ (Q),
# and if so effectively apply Q's square form to A without additional shenanigans; and
# (2) if the preceding dimensions do not match, check whether the appropriate dimension of
# A instead matches the number of rows of the matrix of which Q is a factor (i.e.
# size(Q.factors, 1)), and if so implicitly apply Q's truncated form to A by zero extending
# A as necessary for check (1) to pass (if possible) and then applying Q's square form
#
function *(A::StridedVecOrMat, Q::LQPackedQ)
    TR = promote_type(eltype(A), eltype(Q))
    if size(A, 2) == size(Q.factors, 2)
        C = copy_oftype(A, TR)
    elseif size(A, 2) == size(Q.factors, 1)
        C = zeros(TR, size(A, 1), size(Q.factors, 2))
        copyto!(C, 1, A, 1, length(A))
    else
        _rightappdimmismatch("columns")
    end
    return mul!(C, AbstractMatrix{TR}(Q))
end
function *(adjA::Adjoint{<:Any,<:StridedMatrix}, Q::LQPackedQ)
    A = adjA.parent
    TR = promote_type(eltype(A), eltype(Q))
    if size(A, 1) == size(Q.factors, 2)
        C = adjoint!(similar(A, TR, reverse(size(A))), A)
    elseif size(A, 1) == size(Q.factors, 1)
        C = zeros(TR, size(A, 2), size(Q.factors, 2))
        adjoint!(view(C, :, 1:size(A, 1)), A)
    else
        _rightappdimmismatch("rows")
    end
    return mul!(C, AbstractMatrix{TR}(Q))
end
_rightappdimmismatch(rowsorcols) =
    throw(DimensionMismatch(string("the number of $(rowsorcols) of the matrix on the left ",
        "must match either (1) the number of columns of the (LQPackedQ) matrix on the right ",
        "or (2) the number of rows of that (LQPackedQ) matrix's internal representation ",
        "(the factorization's originating matrix's number of rows)")))

function (\)(A::LQ{TA}, b::StridedVector{Tb}) where {TA,Tb}
    S = promote_type(TA,Tb)
    m = checksquare(A)
    m == length(b) || throw(DimensionMismatch("left hand side has $m rows, but right hand side has length $(length(b))"))
    AA = Factorization{S}(A)
    x = ldiv!(AA, copy_oftype(b, S))
    return x
end
function (\)(A::LQ{TA},B::StridedMatrix{TB}) where {TA,TB}
    S = promote_type(TA,TB)
    m = checksquare(A)
    m == size(B,1) || throw(DimensionMismatch("left hand side has $m rows, but right hand side has $(size(B,1)) rows"))
    AA = Factorization{S}(A)
    X = ldiv!(AA, copy_oftype(B, S))
    return X
end
# With a real lhs and complex rhs with the same precision, we can reinterpret
# the complex rhs as a real rhs with twice the number of columns
function (\)(F::LQ{T}, B::VecOrMat{Complex{T}}) where T<:BlasReal
    c2r = reshape(transpose(reinterpret(T, reshape(B, (1, length(B))))), size(B, 1), 2*size(B, 2))
    x = ldiv!(F, c2r)
    return reshape(collect(reinterpret(Complex{T}, transpose(reshape(x, div(length(x), 2), 2)))),
                           isa(B, AbstractVector) ? (size(F,2),) : (size(F,2), size(B,2)))
end

function ldiv!(A::LQ{T}, B::StridedVecOrMat{T}) where T
    mul!(Adjoint(A.Q), ldiv!(LowerTriangular(A.L),B))
    return B
end

ext: jl

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
using Compat
using Distances
using Ripserer
using SparseArrays
using Random
using StaticArrays
using Suppressor
using Test

using Ripserer:
    distances, births, adjacency_matrix, edges, nv, unsafe_simplex, ChainElement, Chain

include("../testdatasets.jl")
include("interfacetest.jl")

@testset "distances" begin
    for points in (
        [(0, 0), (0, 1), (1, 1), (1, 0)],
        [SVector(0, 0), SVector(0, 1), SVector(1, 1), SVector(1, 0)],
    )
        @test distances(points) ≈ [0 1 √2 1; 1 0 1 √2; √2 1 0 1; 1 √2 1 0]
        @test distances(points, Cityblock()) == [0 1 2 1; 1 0 1 2; 2 1 0 1; 1 2 1 0]
    end
end

@testset "Rips points constructor, sparse=false" begin
    filtration = Rips([(sin(x), cos(x)) for x in range(0, 2π; length=101)[1:(end - 1)]])
    adj = adjacency_matrix(filtration)
    @test all(x -> x > 0, adj[i, j] for i in 1:100 for j in (i + 1):100)
    @test eltype(edges(filtration)) === Simplex{1,Float64,Int}

    filtration = Rips{Int32}([
        (sin(x), cos(x)) for x in range(0.0f0, 2.0f0π; length=101)[1:(end - 1)]
    ])
    adj = adjacency_matrix(filtration)
    @test !issparse(adj)
    @test all(x -> x > 0, adj[i, j] for i in 1:100 for j in (i + 1):100)
    @test eltype(edges(filtration)) === Simplex{1,Float32,Int32}
end

@testset "Rips points constructor, sparse=true" begin
    filtration = Rips(
        [(sin(x), cos(x)) for x in range(0, 2π; length=101)[1:(end - 1)]];
        threshold=0.1,
        sparse=true,
    )
    adj = adjacency_matrix(filtration)
    @test issparse(adj)
    @test maximum(adj) ≤ 0.1
    @test maximum(adj) ≤ threshold(filtration)
end

@testset "Construction does not alter input" begin
    dist = [
        0 9 1 2
        9 0 3 4
        1 3 0 4
        2 4 4 0
    ]
    orig_dist = copy(dist)
    Rips(dist; threshold=1)
    @test dist == orig_dist
    Rips(dist; threshold=1, sparse=true)
    @test dist == orig_dist

    dist = sparse(dist)
    orig_dist = copy(dist)
    Rips(dist; threshold=1)
    @test dist == orig_dist
end

@testset "Warns with duplicate points" begin
    pts = [(1, 0), (1, 1), (0, 1), (0, 0), (0, 0)]
    @test @capture_err(Rips(pts)) ≠ ""
    @test @capture_err(Rips(pts[1:4])) == ""
    @test nv(@suppress Rips(pts)) == 4
end

@testset "Errors" begin
    @testset "Non-square matrices throw an error" begin
        @test_throws DimensionMismatch Rips(zeros(3, 2))
        @test_throws DimensionMismatch Rips(zeros(3, 2); sparse=true)
    end
    @testset "Asymmetric matrices throw an error" begin
        @test_throws ArgumentError Rips([1 1 1; 1 1 1; 1 2 1])
        @test_throws ArgumentError Rips([1 1 1; 2 1 1; 1 2 1]; sparse=true)
    end
    @testset "Edge births must be larger than vertex births" begin
        @test_throws ArgumentError Rips([1 1 1; 1 1 1; 1 1 2])
        @test_throws ArgumentError Rips([1 1 1; 1 2 1; 1 1 1]; sparse=true)
    end
end

@testset "ripserer" begin
    @testset "Dense" begin
        @testset "Icosahedron" begin
            d0, d1, d2 = ripserer(icosahedron; dim_max=2)
            @test d0 == [fill((0.0, 1.0), 11); (0.0, Inf)]
            @test d1 == []
            @test d2 == [(1.0, 2.0)]
        end
        @testset "Cycle with various fields" begin
            d0_2, d1_2, d2_2, d3_2 = ripserer(Rips{Int32}, cycle; dim_max=3)
            d0_7, d1_7, d2_7, d3_7 = ripserer(Rips(cycle); dim_max=3, modulus=7)
            d0_r, d1_r, d2_r, d3_r = ripserer(cycle; dim_max=3, field=Rational{Int})

            @test d0_2 == d0_7 == d0_r == [fill((0, 1), size(cycle, 1) - 1); (0, Inf)]
            @test d1_2 == d1_7 == d1_r == [(1, 6)]
            @test d2_2 == d2_7 == d2_r == fill((6, 7), 5)
            @test d3_2 == d3_7 == d3_r == [(7, 8)]
        end
        @testset "RP2 with various fields" begin
            _, d1_2, d2_2 = ripserer(Rips(projective_plane); dim_max=2)
            _, d1_3, d2_3 = ripserer(Rips, projective_plane; dim_max=2, modulus=3)
            _, d1_331, d2_331 = ripserer(projective_plane; dim_max=2, field=Mod{5})
            _, d1_r, d2_r = ripserer(projective_plane; dim_max=2, field=Rational{Int})
            @test d1_2 == [(1, 2)]
            @test d2_2 == [(1, 2)]
            @test d1_3 == d1_331 == d1_r == []
            @test d2_3 == d2_331 == d1_r == []
        end
        @testset "Icosahedron, threshold=1" begin
            d0, d1, d2 = ripserer(Rips(icosahedron; threshold=1); dim_max=2)
            @test d0 == [fill((0.0, 1.0), 11); (0.0, Inf)]
            @test d1 == []
            @test d2 == [(1.0, Inf)]
        end
        @testset "Icosahedron, threshold=0.5" begin
            d0, d1, d2 = ripserer(icosahedron; dim_max=2, threshold=0.5)
            @test d0 == fill((0.0, Inf), 12)
            @test d1 == []
            @test d2 == []
        end
        @testset "RP2 with various fields, threshold=1" begin
            _, d1_2, d2_2 = ripserer(Rips, projective_plane; dim_max=2, threshold=1)
            _, d1_3, d2_3 = ripserer(projective_plane; dim_max=2, modulus=3, threshold=1)
            _, d1_331, d2_331 = ripserer(
                projective_plane; dim_max=2, field=Mod{5}, threshold=1
            )
            _, d1_r, d2_r = ripserer(
                projective_plane; dim_max=2, field=Rational{Int}, threshold=1
            )
            @test d1_2 == [(1, Inf)]
            @test d2_2 == [(1, Inf)]
            @test d1_3 == d1_331 == d1_r == []
            @test d2_3 == d2_331 == d1_r == []
        end
        @testset "Points as input" begin
            for metric in (Euclidean(), Cityblock())
                pts = torus_points(9)
                @test ripserer(pts; metric=metric) == ripserer(Rips(pts; metric=metric))
            end
        end
        @testset "Cutoff" begin
            d0, d1 = ripserer(rand_dist_matrix(20); cutoff=0.5)
            @test all(persistence.(d0) .> 0.5)
            @test all(persistence.(d1) .> 0.5)
        end
    end

    @testset "Sparse" begin
        @testset "Icosahedron" begin
            d0, d1, d2 = ripserer(sparse(icosahedron); dim_max=2)
            @test d0 == [fill((0.0, 1.0), 11); (0.0, Inf)]
            @test d1 == []
            @test d2 == [(1.0, 2.0)]
        end
        @testset "RP2 with various fields, threshold=1" begin
            _, d1_2, d2_2 = ripserer(sparse(projective_plane); dim_max=2, threshold=1)
            _, d1_3, d2_3 = ripserer(
                Rips(projective_plane; threshold=1, sparse=true); dim_max=2, modulus=3
            )
            _, d1_331, d2_331 = ripserer(
                Rips, sparse(projective_plane); dim_max=2, field=Mod{5}, threshold=1
            )
            _, d1_r, d2_r = ripserer(
                Rips,
                projective_plane;
                sparse=true,
                dim_max=2,
                field=Rational{Int},
                threshold=1,
            )
            @test d1_2 == [(1, Inf)]
            @test d2_2 == [(1, Inf)]
            @test d1_3 == d1_331 == d1_r == []
            @test d2_3 == d2_331 == d1_r == []
        end
        @testset "Equal to Rips" begin
            for thresh in (nothing, 1, 0.5, 0.126)
                data = torus_points(100)
                r_res = ripserer(data; threshold=thresh, dim_max=2)
                s_res_1 = ripserer(Rips, data; threshold=thresh, sparse=true, dim_max=2)

                # Add zeros to diagonal. Adding ones first actually changes the structure of
                # the matrix.
                data2 = sparse(Ripserer.distances(data))
                for i in axes(data2, 1)
                    data2[i, i] = 1
                    data2[i, i] = 0
                end
                s_res_2 = ripserer(data2; threshold=thresh, dim_max=2)

                @test r_res == s_res_1 == s_res_2
            end
        end
    end

    @testset "Representatives" begin
        @testset "Known example" begin
            # This example was generated by getting representatives from ripser.
            _, d1, d2 = ripserer(projective_plane; dim_max=2, reps=true)

            @test simplex.(representative(only(d1))) == [
                Simplex{1}((11, 10), 1),
                Simplex{1}((10, 7), 1),
                Simplex{1}((10, 6), 1),
                Simplex{1}((8, 1), 1),
                Simplex{1}((7, 3), 1),
                Simplex{1}((7, 1), 1),
                Simplex{1}((6, 2), 1),
                Simplex{1}((5, 1), 1),
                Simplex{1}((2, 1), 1),
            ]
            @test coefficient.(representative(only(d1))) == fill(Mod{2}(1), 9)
            @test simplex.(representative(only(d2))) == [Simplex{2}((6, 2, 1), 1)]
            @test coefficient.(representative(only(d2))) == [Mod{2}(1)]
        end
        @testset "Types" begin
            d0, d1, d2, d3 = ripserer(cycle; dim_max=3, reps=true)
            @test d1[1].representative isa Chain

            d0, d1, d2, d3 = ripserer(cycle; dim_max=3, reps=true, field=Rational{Int})
            @test d3[1].representative isa Chain
        end
        @testset "Infinite interval very low threshold" begin
            _, d1 = ripserer(cycle; dim_max=1, reps=true, threshold=1, field=Rational{Int})
            rep = representative(only(d1))
            @test simplex(only(rep)) == birth_simplex(only(d1))
            @test rep isa Chain{Rational{Int},Simplex{1,Int,Int}}
        end
        @testset "Infinite interval higher threshold" begin
            _, d1 = ripserer(cycle; dim_max=1, reps=true, threshold=3, field=Rational{Int})
            rep = representative(only(d1))
            @test !isempty(rep)
            @test rep isa Chain{Rational{Int},Simplex{1,Int,Int}}
        end
        @testset "Critical simplices" begin
            result = ripserer(torus(100); reps=true, threshold=0.5)
            for diag in result
                @test birth.(diag) == birth.(birth_simplex.(diag))
                finite = filter(isfinite, diag)
                @test death.(finite) == birth.(death_simplex.(finite))
                infinite = filter(!isfinite, diag)
                @test all(isnothing, death_simplex.(infinite))
            end
        end
    end

    @testset "Diagram metadata" begin
        filtration = Rips(cycle)
        d0, d1, d2, d3 = ripserer(filtration; dim_max=3, reps=true, field=Rational{Int})
        @test d0.dim == 0
        @test d1.dim == 1
        @test d2.dim == 2
        @test d3.dim == 3
        thresh = Float64(threshold(Rips(cycle)))
        @test d0.threshold ≡ thresh
        @test d1.threshold ≡ thresh
        @test d2.threshold ≡ thresh
        @test d3.threshold ≡ thresh
        field = Rational{Int}
        @test d0.field ≡ field
        @test d1.field ≡ field
        @test d2.field ≡ field
        @test d3.field ≡ field
        @test d0.filtration == filtration
        @test d1.filtration == filtration
        @test d2.filtration == filtration
        @test d3.filtration == filtration
    end

    @testset "Zero-dimensional sublevel set persistence" begin
        @testset "with sparse matrix" begin
            data = [1, 0, 1, 2, 3, 4, 3, 2, 3, 2, 1, 2]

            # Create distance matrix from data, where neighboring points are connected by edges
            # and the edge weights are equal to the max of both vertex births.
            n = length(data)
            dists = spzeros(n, n)
            for i in 1:n
                dists[i, i] = data[i]
            end
            for i in 1:(n - 1)
                j = i + 1
                dists[i, j] = dists[j, i] = max(dists[i, i], dists[j, j])
            end
            # 0-dimensional persistence should find values of minima and maxima of our data.
            d0 = ripserer(dists; dim_max=0)[1]
            @test d0 == [(2, 3), (1, 4), (0, Inf)]
        end

        @testset "with Rips" begin
            data = [1, 0, 1, 2, 3, 4, 3, 2, 3, 2, 1, 2]

            n = length(data)
            dists = zeros(Int, (n, n))
            for i in 1:n, j in 1:n
                if abs(i - j) ≤ 1
                    dists[i, j] = max(data[i], data[j])
                else
                    dists[i, j] = 5
                end
            end
            d0 = ripserer(dists; dim_max=0)[1]
            @test d0 == [(2, 3), (1, 4), (0, Inf)]
        end
    end

    @testset "Homology and explicit cohomology" begin
        @testset "Representative cycle" begin
            res_hom = ripserer(cycle; alg=:homology, reps=true, dim_max=3)
            @test vertices.(simplex.(representative(res_hom[2][1]))) ==
                sort!(vcat([(i + 1, i) for i in 1:17], [(18, 1)]))
        end
        @testset "Infinite intervals" begin
            @test_broken ripserer(
                Rips(cycle; threshold=2); alg=:homology, implicit=true
            )[2][1] == (1.0, Inf)
            @test_broken ripserer(
                Rips, cycle; alg=:homology, threshold=2, implicit=false
            )[2][1] == (1.0, Inf)
            @test ripserer(cycle; alg=:involuted, threshold=2)[2][1] == (1.0, Inf)
        end
    end

    @testset "Interface" begin
        test_filtration(Rips, cycle; modulus=3, dim_max=2)
    end

    @testset "Errors" begin
        @testset "Overflow checking" begin
            @test_throws OverflowError ripserer(Rips{Int16}(ones(1000, 1000)))
        end
        @testset "Unsupported algirithms" begin
            @test_throws ArgumentError ripserer(ones(5, 5); alg=:something)
        end
        @testset "Int or Float64 field type" begin
            @test_throws ErrorException ripserer(ones(5, 5); field=Int)
            @test_throws ErrorException ripserer(ones(5, 5); field=Float64)
            @test_throws ErrorException ripserer(ones(5, 5); field=UInt8)
        end
        @testset "Explicit cohomology reperesentatives unsupported" begin
            @test_throws ErrorException ripserer(cycle; implicit=false, reps=true)
        end
    end
end


text:
# ----------------------------- HELPER FUNCTIONS ----------------------------- #

function _getxyz(plm::PairwiseListMatrix{T,D,TV}) where {T,D,TV}
    names = getlabels(plm)
    names, names, Matrix(plm)
end

function _getxyz(nplm::NamedArray{T,2,PairwiseListMatrix{T,D,TV},DN}) where {T,D,TV,DN}
    names = getlabels(nplm)
    names, names, Matrix(nplm.array)
end

# ---------------------------------- MATRIX ---------------------------------- #

@recipe function plot(plm::PairwiseListMatrix{T,D,TV}) where {T,D,TV}
    seriestype  :=  :heatmap
    yflip --> true
    ratio --> :equal
    _getxyz(plm)
end

@recipe function plot(nplm::NamedArray{T,2,PairwiseListMatrix{T,D,TV},DN}) where {T,D,TV,DN}
    seriestype  :=  :heatmap
    yflip --> true
    ratio --> :equal
    _getxyz(nplm)
end


text:
# This file was generated, do not modify it. # hide
using UrlDownload
using DataFrames
using PrettyPrinting
using PyPlot
using MLJ
ioff() # hide
MLJ.color_off(); # hide

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
function rna(dna)
    rna_str = ""

    # Your code goes here

    return rna_str
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
"""
	bm3d_wie(img::Array{Float64}, imgBasic::Array{Float64}, sigma::AbstractFloat, config::bm3d_config)

2nd step of BM3DDenoise

img: input noisy image
imgBasic: denoised image from first step of BM3D (hard thresholding)
sigma: known or assumed standard deviation of noise
config: bm3d_config
"""
function bm3d_wie(
	img::Array{Float64},
	imgBasic::Array{Float64},
	sigma::AbstractFloat,
	config::bm3d_config,
)

	# parameters
	patchSize = config.wie_patchSize
	searchStride = config.wie_searchStride
	searchWin = config.wie_searchWin
	nMatch = config.wie_nMatch

	# block matching step
	config.show_info && @info "2st get_reference_pixels"
	if config.show_info
		refIndex = @time get_reference_pixels(
			[size(img, 1); size(img, 2)],
			patchSize,
			searchStride,
		)
	else
		refIndex = get_reference_pixels(
			[size(img, 1); size(img, 2)],
			patchSize,
			searchStride,
		)
	end
	config.show_info && @info "2st get_reference_pixels end"

	config.show_info && @info "2st match_patches"
	if config.show_info
		matchTable = @time match_patches(
			imgBasic,
			refIndex,
			patchSize,
			searchWin,
			nMatch,
			config.wie_distFunction,
		)
	else
		matchTable = match_patches(
			imgBasic,
			refIndex,
			patchSize,
			searchWin,
			nMatch,
			config.wie_distFunction,
		)
	end
	config.show_info && @info "2st match_patches end"

	# Don't use similar(), because it need initialize to 0.0
	Wout = zeros(Float64, size(img)...)
	imgOut = zeros(Float64, size(img)...)

	# 3D filtering
	config.show_info && @info "2st 3D filtering"
	if config.show_info
		@time wie_3D_filtering!(
			Wout,
			imgOut,
			img,
			imgBasic,
			matchTable,
			refIndex,
			patchSize,
			nMatch,
			sigma,
			config,
		)
	else
		wie_3D_filtering!(
			Wout,
			imgOut,
			img,
			imgBasic,
			matchTable,
			refIndex,
			patchSize,
			nMatch,
			sigma,
			config,
		)
	end
	config.show_info && @info "2st 3D filtering end"

	return imgOut ./ Wout

end

"""
3D filtering
"""
function wie_3D_filtering!(
	Wout::AbstractArray{<:AbstractFloat,2},
	imgOut::AbstractArray{<:AbstractFloat,2},
	img::AbstractArray{<:AbstractFloat,2},
	imgBasic::AbstractArray{<:AbstractFloat,2},
	matchTable::Array{CartesianIndex{2},3},
	refIndex::Array{CartesianIndex{2},2},
	patchSize::CartesianIndex{2},
	nMatch::Int,
	sigma::AbstractFloat,
	config::bm3d_config,
)
	# Each reference block is processed to reduce memory usage
	I_end, J_end = size(refIndex)

	# Preventing conflicts in group_to_image!
	imgLockPool = Array{ReentrantLock}(undef, I_end, J_end)
	for i in 1:length(imgLockPool)
		imgLockPool[i] = ReentrantLock()
	end

	@views @inbounds Threads.@threads for J in 1:J_end
		# Preventing conflicts in parallel computing
		G3D = zeros(Float64, nMatch + 1, patchSize[1], patchSize[2])
		G3Dbasic = similar(G3D)
		WC = similar(G3Dbasic)
		for I in 1:I_end
			# Compute 3D group spectrum
			form_group!(
				G3D,
				img,
				matchTable,
				refIndex,
				patchSize,
				CartesianIndex(I, J),
				config.wie_group_transform!,
			)
			form_group!(
				G3Dbasic,
				imgBasic,
				matchTable,
				refIndex,
				patchSize,
				CartesianIndex(I, J),
				config.wie_group_transform!,
			)

			# Wiener filtering of 3D groups, using basic estimate as target spectrum
			WC .= G3Dbasic .^ 2 ./ (G3Dbasic .^ 2 .+ sigma^2)
			G3D .*= WC

			# Weight
			T = norm(WC, 1)
			if T > zero(T)
				W = 1.0 / (T * sigma^2)
				G3D .*= W
			else
				W = 1.0
			end

			invert_group!(
				imgOut,
				G3D,
				matchTable,
				refIndex,
				patchSize,
				CartesianIndex(I, J),
				config.wie_group_itransform!,
				imgLockPool,
			)
			group_to_image!(
				Wout,
				W,
				matchTable,
				refIndex,
				patchSize,
				CartesianIndex(I, J),
				imgLockPool,
			)
		end
	end
end

function wie_3D_filtering!(
	Wout::Array{<:AbstractFloat,3},
	imgOut::Array{<:AbstractFloat,3},
	img::Array{<:AbstractFloat,3},
	imgBasic::Array{<:AbstractFloat,3},
	matchTable::Array{CartesianIndex{2},3},
	refIndex::Array{CartesianIndex{2},2},
	patchSize::CartesianIndex{2},
	nMatch::Int,
	sigma::AbstractFloat,
	config::bm3d_config,
)
	@views @inbounds for i in 1:size(img, 3)
		wie_3D_filtering!(
			Wout[:, :, i],
			imgOut[:, :, i],
			img[:, :, i],
			imgBasic[:, :, i],
			matchTable,
			refIndex,
			patchSize,
			nMatch,
			sigma,
			config,
		)
	end
end


text:
# open Demo/ViewController/ViewController.xcworkspace
# Run

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

using Swifter

vc = initial("http://localhost:8080")

@testset "View Controller" begin

    @test 320 == @query vc.view.frame.size.width

    @query vc.label.text = "Hello world"
    @test "Hello world" == @query vc.label.text
    @test "Hello world" == @query vc.view.subviews[0].text

    @query vc.label.alpha = 0.5
    @test 0.5 == @query vc.label.alpha
    @test 0.5 == @query vc.view.subviews[0].alpha

    @query vc.label.alpha = 1
    @test 1 == @query vc.label.alpha
    @test 1 == @query vc.view.subviews[0].alpha

    @query vc.label.userInteractionEnabled = false
    @test false == @query vc.label.userInteractionEnabled
    @test false == @query vc.view.subviews[0].userInteractionEnabled

    @query vc.label.userInteractionEnabled = true
    @test true == @query vc.label.userInteractionEnabled
    @test true == @query vc.view.subviews[0].userInteractionEnabled

    @test 3 == @query vc.view.subviews.count
end

@testset "Application" begin

    @test 320 == @query UIApplication.sharedApplication().delegate.window.frame.size.width
    @query UIApplication.sharedApplication().delegate
    @query UIApplication.sharedApplication().delegate.window
    @query UIApplication.sharedApplication().keyWindow
    @query UIApplication.sharedApplication().windows
    @query UIApplication.sharedApplication().backgroundTimeRemaining

end

@testset "Bundle" begin

    @test "1" == @query NSBundle.mainBundle.infoDictionary["CFBundleVersion"]
    @test "6.0" == @query NSBundle.mainBundle.infoDictionary["CFBundleInfoDictionaryVersion"]
    @test AbstractString["iPhoneSimulator"] == @query NSBundle.mainBundle.infoDictionary["CFBundleSupportedPlatforms"]
    @test "BuildMachineOSBuild" == @query NSBundle.mainBundle.infoDictionary.keys.sort.first
    @test "UISupportedInterfaceOrientations" == @query NSBundle.mainBundle.infoDictionary.keys.sort.last

    @query NSBundle.mainBundle.bundleIdentifier
    @test "com.factor.ViewController" == @query NSBundle.mainBundle.infoDictionary["CFBundleIdentifier"]
    @query NSBundle.mainBundle.infoDictionary.keys.sort
    @query NSBundle.mainBundle.infoDictionary.keys[0]
    @query NSBundle.mainBundle.infoDictionary.keys[1]
    @query NSBundle.mainBundle.infoDictionary.values

    @test "CFBundleSupportedPlatforms" in valueof(@query NSBundle.mainBundle.infoDictionary.keys)
    @test "CFBundleSupportedPlatforms" in keys(valueof(@query NSBundle.mainBundle.infoDictionary))

end

max_stars_repo_name: wookay/Console

max_issues_repo_name: wookay/Console

max_forks_repo_name: wookay/Console


text:
using HTTPDownloads

using Test
using Sockets

using HTTP
import MbedTLS
import JSON

# Find a free port on `network_interface`
function find_free_port(network_interface)
    # listen on port 0 => kernel chooses a free port. See, for example,
    # https://stackoverflow.com/questions/44875422/how-to-pick-a-free-port-for-a-subprocess
    server = listen(network_interface, 0)
    _, free_port = getsockname(server)
    close(server)
    # The kernel can reuse free_port here after $some_time_delay, but apprently
    # this is large enough for Selenium to have used this technique for ten
    # years...
    return free_port
end

# Return SSL config for running a server under self-signed certs
function test_server_ssl_config()
    MbedTLS.SSLConfig(joinpath(@__DIR__, "resources", "server.pem"),
                      joinpath(@__DIR__, "resources", "server.key"))
end

# Return SSL config for running client
function test_client_ssl_config()
    MbedTLS.SSLConfig()

    conf = MbedTLS.SSLConfig()
    MbedTLS.config_defaults!(conf)

    entropy = MbedTLS.Entropy()
    rng = MbedTLS.CtrDrbg()
    MbedTLS.seed!(rng, entropy)
    MbedTLS.rng!(conf, rng)

    MbedTLS.authmode!(conf, MbedTLS.MBEDTLS_SSL_VERIFY_REQUIRED)

    # Trust the test root CA that we self-signed
    root_ca_cert = MbedTLS.crt_parse_file(joinpath(@__DIR__, "resources", "root_ca.crt"))
    MbedTLS.ca_chain!(conf, root_ca_cert)

    return conf
end

function jsonbody(res::HTTP.Response)
    JSON.parse(String(copy(res.body)))
end

#-------------------------------------------------------------------------------

HTTPDownloads.set_downloads_backend()

bingoserver = "https://httpbingo.julialang.org"

@testset "Basic functionality" begin
    res = HTTP.get("$bingoserver/get?a=1&b=2")
    @test res isa HTTP.Response
    body = jsonbody(res)
    @test body["args"]["a"] == ["1"]
    @test body["args"]["b"] == ["2"]

    # Errors mapped to HTTP StatusError
    @test_throws HTTP.ExceptionRequest.StatusError HTTP.get("$bingoserver/status/400")
    @test HTTP.status(HTTP.get("$bingoserver/status/400", status_exception=false)) == 400
end

@testset "Input body data" begin
    postbody_roundtrip(body) = jsonbody(HTTP.post("$bingoserver/post", body=body))["data"]

    data = "hi\x00\x01"
    @test postbody_roundtrip(data) == data
    @test postbody_roundtrip(Vector{UInt8}(codeunits(data))) == data

    # Streams
    @test postbody_roundtrip(IOBuffer(data)) == data
    @test postbody_roundtrip(IOBuffer(codeunits(data))) == data
    @test begin
        io = Base.BufferStream()
        write(io, data)
        close(io) # Should be closewrite()? But that doesn't exist until Julia-1.8
        postbody_roundtrip(io)
    end == data
end

@testset "Output body data" begin
    @test begin
        # Check explicitly provided streams work.
        # Downloads.jl is arguably more correct here than HTTP itself. See:
        # https://github.com/JuliaWeb/HTTP.jl/issues/543
        stream = IOBuffer()
        HTTP.get("$bingoserver/get?a=1", response_stream=stream)
        JSON.parse(String(take!(stream)))["args"]["a"]
    end == ["1"]
    # Right now, HTTP.open() falls through to HTTP.jl internals rather than
    # using Downloads.jl. Check that this works:
    @test begin
        body = nothing
        HTTP.open("GET", "$bingoserver/get?xx=100") do http
            body = String(read(http))
        end
        JSON.parse(body)["args"]["xx"]
    end == ["100"]
end

@testset "Timeouts" begin
    @test HTTP.get("$bingoserver/delay/2", readtimeout=10) isa HTTP.Response
    @test_throws HTTP.TimeoutRequest.ReadTimeoutError HTTP.get("$bingoserver/delay/5", readtimeout=1)

    # Test connection timeout by attempting to connect to an unroutable IP
    # address. Discussion at
    # https://stackoverflow.com/questions/100841/artificially-create-a-connection-timeout-error/37465639
    @test_throws HTTP.TimeoutRequest.ReadTimeoutError HTTP.get("http://10.255.255.1", connect_timeout=1)
end

@testset "Redirects" begin
    # Test that redirects can be disabled
    @test begin
        response = HTTP.get("$bingoserver/redirect/1", redirect=false)
        haskey(Dict(HTTP.headers(response)), "location")
    end
    @test begin
        response = HTTP.get("$bingoserver/redirect/1", redirect=true)
        !haskey(Dict(HTTP.headers(response)), "location")
    end
end

@testset "SSL Config" begin
    test_port = find_free_port(Sockets.localhost)
    server = listen(Sockets.localhost, test_port)

    # Run a local https server with self-signed certificate
    tsk = @async try
        HTTP.listen(Sockets.localhost, test_port; server=server,
                    sslconfig=test_server_ssl_config()) do http
            HTTP.setstatus(http, 200)
            HTTP.startwrite(http)
            HTTP.write(http, "hello, world")
        end
    finally exc
        close(server)
        @error "Error running server" exception=exc,catch_backtrace()
    end

    # Test that we can access this server if we trust our test CA
    response = HTTP.get("https://localhost:$test_port",
                        sslconfig=test_client_ssl_config(),
                        require_ssl_verification=true,
                        retries=0)
    @test String(response.body) == "hello, world"
    @test response.status == 200

    response = HTTP.get("https://localhost:$test_port",
                        require_ssl_verification=false,
                        retries=0)
    @test String(response.body) == "hello, world"
    @test response.status == 200

    close(server)
end


text:
# Weighted row-replication
function replicate{T<:Real}(data::Matrix{T}, w::Vector{T})
    num_samples, num_signals = size(data)
    length(w) == num_samples ||
        throw(DimensionMismatch("Inconsistent array lengths."))
    replicated = zeros(int(sum(w)), num_signals)
    j = 1
    for i = 1:num_samples
        for k = j:(j+int(w[i])-1)
            replicated[k,:] = data[i,:]
        end
        j += int(w[i])
    end
    replicated
end

# Recombine over weights
function recombine{T<:Real}(v::Vector{T}, w::Vector{T})
    num_weights = length(w)
    recombined = zeros(num_weights)
    j = 1
    for i = 1:num_weights
        recombined[i] = v[j] * w[i]
        j += int(w[i])
    end
    recombined
end


text:
using DataDeps

register(DataDep(
  "thyroid-disease-sick-euthyroid",
  "http://archive.ics.uci.edu/ml/datasets/Thyroid+Disease",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/thyroid-disease/sick-euthyroid.data",
  "6b93c0a10c6a1760d47864eb9f4592750cbfac389d9a64b9feec17f2690374c4",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:26,
        categoric_indices=[3:15; 17; 19; 21; 23; 25],
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
using Documenter, Logging, SkipLists, Test
using Base.Threads: nthreads

@info "Running tests with $(nthreads()) threads"

include("list_test_utils.jl")

include("test_node_nonconcurrent.jl")
include("test_list_nonconcurrent.jl")


text:
using GPMLj, Test, Random, PyCall

Random.seed!(123)

@testset "MeanFunctions" begin

    @testset "Additive" begin
        identity(X) = X
        mf = gpflow.mean_functions.Additive(identity, identity)
        @test typeof(mf)<:gpflow.AbstractMeanFunction
        @test typeof(mf)<:gpflow.mean_functions.Additive
        @test typeof(mf.o)<:PyObject
        temp = mf.o
        @test instantiate!(mf) == temp
    end

    @testset "Constant" begin
        mf = gpflow.mean_functions.Constant()
        @test typeof(mf)<:gpflow.AbstractMeanFunction
        @test typeof(mf)<:gpflow.mean_functions.Constant
        @test typeof(mf.o)<:PyObject
        temp = mf.o
        @test instantiate!(mf) == temp
    end

    @testset "Identity" begin
        mf = gpflow.mean_functions.Identity()
        @test typeof(mf)<:gpflow.AbstractMeanFunction
        @test typeof(mf)<:gpflow.mean_functions.Identity
        @test typeof(mf.o)<:PyObject
        temp = mf.o
        @test instantiate!(mf) == temp
    end

    @testset "Linear" begin
        mf = gpflow.mean_functions.Linear()
        @test typeof(mf)<:gpflow.AbstractMeanFunction
        @test typeof(mf)<:gpflow.mean_functions.Linear
        @test typeof(mf.o)<:PyObject
        temp = mf.o
        @test instantiate!(mf) == temp
    end

    @testset "MeanFunction" begin
        mf = gpflow.mean_functions.MeanFunction()
        @test typeof(mf)<:gpflow.AbstractMeanFunction
        @test typeof(mf)<:gpflow.mean_functions.MeanFunction
        @test typeof(mf.o)<:PyObject
        temp = mf.o
        @test instantiate!(mf) == temp
    end

    @testset "SwitchedMeanFunction" begin
        mf1 = gpflow.mean_functions.Linear()
        mf2 = gpflow.mean_functions.Linear()
        mf3 = gpflow.mean_functions.Linear()
        mf = gpflow.mean_functions.SwitchedMeanFunction([mf1, mf2, mf3])
        @test typeof(mf)<:gpflow.AbstractMeanFunction
        @test typeof(mf)<:gpflow.mean_functions.SwitchedMeanFunction
        @test typeof(mf.o)<:PyObject
        temp = mf.o
        @test instantiate!(mf) == temp
    end

end #module


text:
#=
test_search:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-09-21
=#

module TestSearch
    # packages
    using Test

    # external modules
    include("RabinKarpStringSearch/rabin_karp_string_search.jl")
    include("BreadthFirstSearch/breadth_first_search.jl")
    include("DepthFirstSearch/depth_first_search.jl")
    include("HillClimbingSearch/hill_climbing_search.jl")
    include("BestFirstSearch/best_first_search.jl")

    # methods
    function test()
        @testset "Search" begin
            @testset "RabinKarpSearch" begin
                @test RabinKarpSearch.rabin_karp_search("Hello sunshine", "sun") == 7
                @test RabinKarpSearch.rabin_karp_search("GEEKS FOR GEEKS", "GEEKS") == 1
                @test RabinKarpSearch.rabin_karp_search("aaaaaaaaaaaa", "bbb") == -1
            end
            @testset "BreadthFirstSearch" begin
                @test_nowarn BFS.main()
            end
            @testset "DepthFirstSearch" begin
                @test_nowarn DFS.main()
            end
            @testset "HillClimbingSearch" begin
                @test_nowarn HCS.main()
            end
            @testset "BestFirstSearch" begin
                @test_nowarn BestFirstSearch.main()
            end
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .TestSearch
    TestSearch.test()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
function dudt!(du, u, p, t)
    du .= 0.0
    J, ll, lr, lpdm, dgl, dgr, γ, uq = p

    nm = uq.nm
    nq = uq.nq

    ncell = size(u, 1)
    nsp = size(u, 2)

    u_ran = zeros(ncell, nsp, 3, nq)
    for i = 1:ncell
        for j = 1:nsp, k = 1:3
            u_ran[i, j, k, :] .= chaos_ran(u[i, j, k, :], uq)
        end

        tmp = @view u_ran[i, :, :, :]
        positive_limiter(tmp, γ, ps.wp ./ 2, uq.op.quad.weights, ll, lr, 1.0)
    end

    f = zeros(ncell, nsp, 3, nm+1)
    for i = 1:ncell, j = 1:nsp
        _f = zeros(3, nq)
        for k = 1:nq
            _f[:, k] .= euler_flux(u_ran[i, j, :, k], γ)[1] ./ J[i]
        end

        for k = 1:3
            f[i, j, k, :] .= ran_chaos(_f[k, :], uq)
        end
    end

    f_face = zeros(ncell, 3, nm+1, 2)
    for i = 1:ncell, j = 1:3, k = 1:nm+1
        f_face[i, j, k, 1] = dot(f[i, :, j, k], lr)
        f_face[i, j, k, 2] = dot(f[i, :, j, k], ll)
    end

    u_face = zeros(ncell, 3, nq, 2)
    for i = 1:ncell, j = 1:3, k = 1:nq
        u_face[i, j, k, 1] = dot(u_ran[i, :, j, k], lr)
        u_face[i, j, k, 2] = dot(u_ran[i, :, j, k], ll)
    end

    for i = 1:ncell, k = 1:nq
        primL = conserve_prim(u_face[i, :, k, 2], γ)
        primR = conserve_prim(u_face[i, :, k, 1], γ)
        if min(primL[end], primR[end]) < 0
            @info "negative interpolated temperature: ($(1/primL[end]), $(1/primR[end]))"

            #=u_mean = [sum(u_ran[i, :, idx, k] .* ps.wp ./ 2) for idx in axes(u_ran, 3)]
            for pp1 = 1:nsp
                u_ran[i, pp1, :, k] .= u_mean
            end=#
            tmp = @view u_ran[i, :, :, :]
            positive_limiter(tmp, γ, ps.wp ./ 2, uq.op.quad.weights, ll, lr, 0.99)

            for j = 1:3
                u_face[i, j, k, 1] = dot(u_ran[i, :, j, k], lr)
                u_face[i, j, k, 2] = dot(u_ran[i, :, j, k], ll)
            end
        end
    end

    fq_interaction = zeros(ncell + 1, 3, nq)
    for i = 2:ncell, j = 1:nq
        #=primL = conserve_prim(u_face[i-1, :, j, 1], γ)
        primR = conserve_prim(u_face[i, :, j, 1], γ)
        if primL[end] < 0
            primL[end] = 1 / 1e-6
            u_face[i-1, :, j, 1] .= prim_conserve(primL, γ)
        end
        if primR[end] < 0
            primR[end] = 1 / 1e-6
            u_face[i, :, j, 1] .= prim_conserve(primR, γ)
        end=#

        fw = @view fq_interaction[i, :, j]
        flux_hll!(fw, u_face[i-1, :, j, 1], u_face[i, :, j, 2], γ, 1.0)
    end

    f_interaction = zeros(ncell + 1, 3, nm+1)
    for i = 2:ncell, j = 1:3
        f_interaction[i, j, :] .= ran_chaos(fq_interaction[i, j, :], uq)
    end

    rhs1 = zero(u)
    for i = 1:ncell, ppp1 = 1:nsp, k = 1:3, l = 1:nm+1
        rhs1[i, ppp1, k, l] = dot(f[i, :, k, l], lpdm[ppp1, :])
    end

    idx = 2:ncell-1 # ending points are Dirichlet
    for i in idx, ppp1 = 1:nsp, k = 1:3, l = 1:nm+1
        du[i, ppp1, k, l] =
            -(
                rhs1[i, ppp1, k, l] +
                (f_interaction[i, k, l] / J[i] - f_face[i, k, l, 2]) * dgl[ppp1] +
                (f_interaction[i+1, k, l] / J[i] - f_face[i, k, l, 1]) * dgr[ppp1]
            )
    end
end

max_stars_repo_name: CSMMLab/FRSG

max_issues_repo_name: CSMMLab/FRSG

max_forks_repo_name: CSMMLab/FRSG


text:
using Test
# Toy dataset.
# Format: each row is an example.
# The last column is the label name, a string.
# The first two columns are features.
# Feel free to play with it by adding more features & examples.
# Interesting note: I've written this so the 2nd and 5th examples
# have the same features, but different labels - so we can see how the
# tree handles this case.

"""
Find the unique values for a column in a dataset.
Could make the result Union{Set{Real}, Set{AbstractString}}???
"""
function unique_vals(rows::Data, col::Int)::Set{Feature}
    Set([row[col] for row in rows])
end

training_data = Data([
    ["Green", 3, "Apple"],
    ["Yellow", 3, "Apple"],
    ["Red", 1, "Grape"],
    ["Red", 1, "Grape"],
    ["Yellow", 3, "Lemon"],
])

"""
    Column labels.
    These are used only to print the tree.
"""
const header = ["color", "diameter", "label"]
Base.show(io::IO, question::Question) = CART.show(io, question, header)
#######
# Demo:
# Let"s look at some example to understand how Gini Impurity works.
#
# First, we"ll look at a dataset with no mixing.
@testset "Gini" begin
    no_mixing = Data([["Apple"], ["Apple"]])
    # this will return 0
    @test gini(no_mixing) == 0.0
    # Now, we"ll look at dataset with a 50:50 apples:oranges ratio
    some_mixing = Data([["Apple"], ["Orange"]])
    # this will return 0.5 - meaning, there"s a 50% chance of misclassifying
    # a random example we draw from the dataset.
    @test gini(some_mixing) == 0.5

    # Now, we'll look at a dataset with many different labels
    lots_of_mixing = Data([["Apple"],
                           ["Orange"],
                           ["Grape"],
                           ["Grapefruit"],
                           ["Blueberry"]])
    # This will return 0.8
    @test gini(lots_of_mixing) ≈ 0.8
end
#######

#######
# Demo:
# Calculate the uncertainy of our training data.
@testset "Information gain" begin
    current_uncertainty = gini(training_data)
    @test current_uncertainty ≈ 0.6399999999999999

    # How much information do we gain by partioning on "Green"?
    true_rows, false_rows = partition(training_data, Question(1, "Green"))
    gain = info_gain(true_rows, false_rows, current_uncertainty)
    @test gain ≈ 0.1399999999999999

    # What about if we partioned on "Red" instead?
    true_rows, false_rows = partition(training_data, Question(1,"Red"))
    gain = info_gain(true_rows, false_rows, current_uncertainty)
    @test gain ≈ 0.37333333333333324
end

@testset "Best split" begin
end

@testset "Print tree" begin
    my_tree = build_tree(training_data)
    print_tree(my_tree)
end


text:
struct LatentTimeSeriesModel{R,E,N,D,R1,E1,D1,T}
    rnn::R1
    enc::E1
    node::N
    dec::D1
    p1::T
    p2::T
    p3::T
    p4::T

    function LatentTimeSeriesModel(rnn, enc, node, dec)
        p1, re1 = Flux.destructure(rnn)
        p2, re2 = Flux.destructure(enc)
        p3 = node.p
        p4, re3 = Flux.destructure(dec)
        return new{
            typeof(rnn),
            typeof(enc),
            typeof(node),
            typeof(dec),
            typeof(re1),
            typeof(re2),
            typeof(re3),
            typeof(p3),
        }(
            re1,
            re2,
            node,
            re3,
            p1,
            p2,
            p3,
            p4,
        )
    end
end

Flux.trainable(m::LatentTimeSeriesModel) = (m.p1, m.p2, m.p3, m.p4)

function (m::LatentTimeSeriesModel{R,E,N,D})(
    x::AbstractArray{T,3},
    p1 = m.p1,
    p2 = m.p2,
    p3 = m.p3,
    p4 = m.p4;
    node_kwargs...,
) where {R,E,N,D,T}
    rnn = m.rnn(p1)::R
    out = rnn(x)

    enc = m.enc(p2)::E
    out = enc(out)

    latent_dim = size(out, 1) ÷ 2
    μ₀ = out[1:latent_dim, :]
    logσ² = out[latent_dim+1:end, :]

    sample = CUDA.randn(size(μ₀, 1), size(μ₀, 2))::CuArray{Float32,2}
    z₀ = sample .* exp.(logσ² / 2) .+ μ₀

    res, nfe, sv = m.node(z₀, p3; node_kwargs...)

    dec = m.dec(p4)::D

    res = reshape(res, size(res, 1), :)
    result = dec(res)
    result = reshape(result, size(result, 1), :, size(x, 3))

    return result, μ₀, logσ², nfe, sv
end


text:
# Julia wrapper for header: /usr/include/libavutil/hwcontext_videotoolbox.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_map_videotoolbox_format_to_pixfmt,
    av_map_videotoolbox_format_from_pixfmt

function av_map_videotoolbox_format_to_pixfmt(cv_fmt::Integer)
    ccall((:av_map_videotoolbox_format_to_pixfmt, libavutil), AVPixelFormat, (UInt32,), cv_fmt)
end

function av_map_videotoolbox_format_from_pixfmt(pix_fmt::AVPixelFormat)
    ccall((:av_map_videotoolbox_format_from_pixfmt, libavutil), UInt32, (AVPixelFormat,), pix_fmt)
end


text:
export GPC

abstract type AbstractGP end

# A collection of GPs (GPC == "GP Collection"). Used to keep track of GPs.
mutable struct GPC
    n::Int
    GPC() = new(0)
end

@nograd GPC

next_index(gpc::GPC) = gpc.n + 1

#
# Projecting AbstractGPs onto a finite dimensional marginal
#

import Base: rand, length
import Distributions: logpdf, ContinuousMultivariateDistribution

export mean, std, cov, marginals, rand, logpdf, elbo, dtc

"""
    FiniteGP{Tf<:AbstractGP, Tx<:AV, TΣy}

The finite-dimensional projection of the AbstractGP `f` at `x`.
"""
struct FiniteGP{Tf<:AbstractGP, Tx<:AV, TΣy} <: ContinuousMultivariateDistribution
    f::Tf
    x::Tx
    Σy::TΣy
end
FiniteGP(f::AbstractGP, x::AV, σ²::AV{<:Real}) = FiniteGP(f, x, Diagonal(σ²))
FiniteGP(f::AbstractGP, x::AV, σ²::Real) = FiniteGP(f, x, Fill(σ², length(x)))
FiniteGP(f::AbstractGP, x::AV) = FiniteGP(f, x, 1e-18)

length(f::FiniteGP) = length(f.x)

"""
    mean(fx::FiniteGP)

Compute the mean vector of `fx`.

```jldoctest
julia> f = GP(Matern52(), GPC());

julia> x = randn(11);

julia> mean(f(x)) == zeros(11)
true
```
"""
mean(fx::FiniteGP) = mean_vector(fx.f, fx.x)

"""
    cov(f::FiniteGP)

Compute the covariance matrix of `fx`.

## Noise-free observations

```jldoctest cov_finitegp
julia> f = GP(Matern52(), GPC());

julia> x = randn(11);

julia> # Noise-free

julia> cov(f(x)) == Stheno.pw(Matern52(), x)
true
```

## Isotropic observation noise

```jldoctest cov_finitegp
julia> cov(f(x, 0.1)) == Stheno.pw(Matern52(), x) + 0.1 * I
true
```

## Independent anisotropic observation noise

```jldoctest cov_finitegp
julia> s = rand(11);

julia> cov(f(x, s)) == Stheno.pw(Matern52(), x) + Diagonal(s)
true
```

## Correlated observation noise

```jldoctest cov_finitegp
julia> A = randn(11, 11); S = A'A;

julia> cov(f(x, S)) == Stheno.pw(Matern52(), x) + S
true
```
"""
cov(f::FiniteGP) = cov(f.f, f.x) + f.Σy

"""
    cov(fx::FiniteGP, gx::FiniteGP)

Compute the cross-covariance matrix between `fx` and `gx`.

```jldoctest
julia> f = GP(Matern32(), GPC());

julia> x1 = randn(11);

julia> x2 = randn(13);

julia> cov(f(x1), f(x2)) == pw(Matern32(), x1, x2)
true
```
"""
cov(fx::FiniteGP, gx::FiniteGP) = cov(fx.f, gx.f, fx.x, gx.x)

"""
    marginals(f::FiniteGP)

Compute a vector of Normal distributions representing the marginals of `f` efficiently.
In particular, the off-diagonal elements of `cov(f(x))` are never computed.

```jldoctest
julia> f = GP(Matern32(), GPC());

julia> x = randn(11);

julia> fs = marginals(f(x));

julia> mean.(fs) == mean(f(x))
true

julia> std.(fs) == sqrt.(diag(cov(f(x))))
true
```
"""
marginals(f::FiniteGP) = Normal.(mean(f), sqrt.(cov_diag(f.f, f.x) .+ diag(f.Σy)))

"""
    rand(rng::AbstractRNG, f::FiniteGP, N::Int=1)

Obtain `N` independent samples from the marginals `f` using `rng`. Single-sample methods
produce a `length(f)` vector. Multi-sample methods produce a `length(f)` x `N` `Matrix`.

```jldoctest
julia> f = GP(Matern32(), GPC());

julia> x = randn(11);

julia> rand(f(x)) isa Vector{Float64}
true

julia> rand(MersenneTwister(123456), f(x)) isa Vector{Float64}
true

julia> rand(f(x), 3) isa Matrix{Float64}
true

julia> rand(MersenneTwister(123456), f(x), 3) isa Matrix{Float64}
true
```
"""
function rand(rng::AbstractRNG, f::FiniteGP, N::Int)
    μ, C = mean(f), cholesky(Symmetric(cov(f)))
    return μ .+ C.U' * randn(rng, promote_type(eltype(μ), eltype(C)), length(μ), N)
end
rand(f::FiniteGP, N::Int) = rand(Random.GLOBAL_RNG, f, N)
rand(rng::AbstractRNG, f::FiniteGP) = vec(rand(rng, f, 1))
rand(f::FiniteGP) = vec(rand(f, 1))

"""
    logpdf(f::FiniteGP, y::AbstractVecOrMat{<:Real})

The logpdf of `y` under `f` if is `y isa AbstractVector`. logpdf of each column of `y` if
`y isa Matrix`.

```jldoctest
julia> f = GP(Matern32(), GPC());

julia> x = randn(11);

julia> y = rand(f(x));

julia> logpdf(f(x), y) isa Real
true

julia> Y = rand(f(x), 3);

julia> logpdf(f(x), Y) isa AbstractVector{<:Real}
true
```
"""
logpdf(f::FiniteGP, y::AbstractVector{<:Real}) = first(logpdf(f, reshape(y, :, 1)))

function logpdf(f::FiniteGP, Y::AbstractMatrix{<:Real})
    μ, C = mean(f), cholesky(Symmetric(cov(f)))
    T = promote_type(eltype(μ), eltype(C), eltype(Y))
    return -((size(Y, 1) * T(log(2π)) + logdet(C)) .+ diag_Xt_invA_X(C, Y .- μ)) ./ 2
end

"""
   elbo(f::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)

The saturated Titsias Evidence LOwer Bound (ELBO) [1]. `y` are observations of `f`, and `u`
are pseudo-points.

```jldoctest
julia> f = GP(Matern52(), GPC());

julia> x = randn(1000);

julia> z = range(-5.0, 5.0; length=13);

julia> y = rand(f(x, 0.1));

julia> elbo(f(x, 0.1), y, f(z)) < logpdf(f(x, 0.1), y)
true
```

[1] - M. K. Titsias. "Variational learning of inducing variables in sparse Gaussian
processes". In: Proceedings of the Twelfth International Conference on Artificial
Intelligence and Statistics. 2009.
"""
function elbo(f::FiniteGP, y::AV{<:Real}, u::FiniteGP)
    _dtc, chol_Σy, A = _compute_intermediates(f, y, u)
    return _dtc - (tr_Cf_invΣy(f, f.Σy, chol_Σy) - sum(abs2, A)) / 2
end

"""
    dtc(f::FiniteGP, y::AV{<:Real}, u::FiniteGP)

The Deterministic Training Conditional (DTC) [1]. `y` are observations of `f`, and `u`
are pseudo-points.

```jldoctest
julia> f = GP(Matern52(), GPC());

julia> x = randn(1000);

julia> z = range(-5.0, 5.0; length=256);

julia> y = rand(f(x, 0.1));

julia> isapprox(dtc(f(x, 0.1), y, f(z)), logpdf(f(x, 0.1), y); atol=1e-3, rtol=1e-3)
true
```

[1] - M. Seeger, C. K. I. Williams and N. D. Lawrence. "Fast Forward Selection to Speed Up
Sparse Gaussian Process Regression". In: Proceedings of the Ninth International Workshop on
Artificial Intelligence and Statistics. 2003
"""
dtc(f::FiniteGP, y::AV{<:Real}, u::FiniteGP) = first(_compute_intermediates(f, y, u))

# Factor out computations common to the `elbo` and `dtc`.
function _compute_intermediates(f::FiniteGP, y::AV{<:Real}, u::FiniteGP)
    consistency_check(f, y, u)
    chol_Σy = cholesky(f.Σy)

    A = cholesky(Symmetric(cov(u))).U' \ (chol_Σy.U' \ cov(f, u))'
    Λ_ε = cholesky(Symmetric(A * A' + I))
    δ = chol_Σy.U' \ (y - mean(f))

    tmp = logdet(chol_Σy) + logdet(Λ_ε) + sum(abs2, δ) - sum(abs2, Λ_ε.U' \ (A * δ))
    _dtc = -(length(y) * typeof(tmp)(log(2π)) + tmp) / 2
    return _dtc, chol_Σy, A
end

function consistency_check(f, y, u)
    @assert length(f) == size(y, 1)
end
Zygote.@nograd consistency_check

import Base: \
\(A::AbstractMatrix, B::Diagonal) = A \ Matrix(B)

\(A::Union{LowerTriangular, UpperTriangular}, B::Diagonal) = A \ Matrix(B)
\(A::Adjoint{<:Any, <:Union{LowerTriangular, UpperTriangular}}, B::Diagonal) = A \ Matrix(B)

# Compute tr(Cf / Σy) efficiently for different types of Σy. For dense Σy you obviously need
# to compute the entirety of Cf, which is bad, but for particular structured Σy one requires
# only a subset of the elements. Σy isa UniformScaling is version usually considered.
function tr_Cf_invΣy(f::FiniteGP, Σy::UniformScaling, chol_Σy::Cholesky)
    return sum(cov_diag(f.f, f.x)) / Σy.λ
end
function tr_Cf_invΣy(f::FiniteGP, Σy::Diagonal, chol_Σy::Cholesky)
    return sum(cov_diag(f.f, f.x) ./ diag(Σy))
end
tr_Cf_invΣy(f::FiniteGP, Σy::Matrix, chol_Σy::Cholesky) = tr(chol_Σy \ cov(f))
# function tr_Cf_invΣy(f::FiniteGP, Σy::BlockDiagonal, chol_Σy::Cholesky)
#     C = cholesky(Symmetric(_get_kernel_block_diag(f, cumulsizes(Σy, 1))))
#     return tr_At_A(chol_Σy.U' \ C.U')
# end

# function _get_kernel_block_diag(f::FiniteGP, cs)
#     k = kernel(f.f)
#     ids = map(n->cs[n]:cs[n+1]-1, 1:length(cs)-1)
#     xs = map(id->f.x[id], ids)
#     Σs = map(x->pw(k, x), xs)
#     return block_diagonal(Σs)
# end

# function _get_kernel_block_diag(f::FiniteGP{<:GP{<:BlockMean, <:BlockKernel}, <:BlockData}, cs)
#     k = kernel(f.f)
#     ids = map(n->cs[n]:cs[n+1]-1, 1:length(cs)-1)
#     @assert _test_block_consistency(ids, f)
#     xs = blocks(f.x)
#     Σs = map(n->pw(k.ks[n], xs[n]), 1:length(xs))
#     return block_diagonal(Σs)
# end

_test_block_consistency(ids, f) = length.(ids) == length.(blocks(f.x))
Zygote.@nograd _test_block_consistency

# """
#     elbo(f::FiniteGP, y::AV{<:Real}, u::FiniteGP, mε::AV{<:Real}, Λε::AM{<:Real})

# The unsaturated Titsias-ELBO.
# """
# function elbo(f::FiniteGP, y::AV{<:Real}, u::FiniteGP, mε::AV{<:Real}, Λε::AM{<:Real})
#     @assert length(u.x) == length(mε)
#     @assert size(Λε) == (length(mε), length(mε))
#     # do stuff.
# end

import Base: |, merge
export ←, |, Obs

"""
    Observation

Represents fixing a paricular (finite) GP to have a particular (vector) value.
"""
struct Observation{Tf<:FiniteGP, Ty<:Vector}
    f::Tf
    y::Ty
end

const Obs = Observation

←(f, y) = Observation(f, y)
get_f(c::Observation) = c.f
get_y(c::Observation) = c.y


text:
function visualize{T, D}(style::Style, img::Texture{T, D, 2}, data::Dict{Symbol, Any})
  kernel = data[:kernel]
  w, h  = img.dims
  texparams = [
     (GL_TEXTURE_MIN_FILTER, GL_NEAREST),
    (GL_TEXTURE_MAG_FILTER, GL_NEAREST),
    (GL_TEXTURE_WRAP_S,  GL_CLAMP_TO_EDGE),
    (GL_TEXTURE_WRAP_T,  GL_CLAMP_TO_EDGE)
  ]

  v, uv, indexes = genquad(0f0, 0f0, w, h)
  if typeof(kernel) <: Real
    filterkernel = float32(kernel)
  elseif eltype(kernel) <: Union(AbstractArray, Real)
    filterkernel = Texture(kernel, parameters=texparams)
  end

  data = @compat Dict(
    :vertex           => GLBuffer(v, 2),
    :index            => indexbuffer(indexes),
    :uv               => GLBuffer(uv, 2),
    :image            => img,
    :normrange        => data[:normrange],
    :filterkernel     => filterkernel,
    :projectionview   => data[:screen].orthographiccam.projectionview,
    :model            => data[:model]
  )

  fragdatalocation = [(0, "fragment_color"),(1, "fragment_groupid")]
  textureshader    = TemplateProgram(joinpath(shaderdir, "uv_vert.vert"), joinpath(shaderdir, "texture.frag"), attributes=data, fragdatalocation=fragdatalocation)

  obj = RenderObject(data, textureshader)

  prerender!(obj, glDisable, GL_DEPTH_TEST, enabletransparency, glDisable, GL_CULL_FACE)
  postrender!(obj, render, obj.vertexarray)
  obj
end


text:
export  premrestotalsd
using Distributions

"Constructs line of business & calculates statistical parameters"
function NLLob(lob::DataFrame,
               β_raw_vec,
               rf,
               corr_prem_res,
               df_lobs_prem_risk)
  prem = Array(Float64, 2)
  name = lob[1, :name]
  ## calculate written net premium and earned net premium
  re_prop_q = [lob[1, :re_prop_q_py],
               lob[1, :re_prop_q_cy],
               lob[1, :re_prop_q_ny]]
  upr = [lob[1, :upr_py], lob[1, :upr_cy], lob[1, :upr_ny]]
  ## upr_gross[PY] not defined, set == Inf:
  upr_gross =    upr ./ (1 .- [1, re_prop_q[PY], re_prop_q[CY]])
  prem_gross_w =
    [lob[1, :prem_gross_w_py], lob[1, :prem_gross_w_cy]]
  prem_w = (1 .- [re_prop_q[PY], re_prop_q[CY]]) .* prem_gross_w
  prem_gross_cy =
    prem_gross_w[CY] + upr_gross[CY] - upr_gross[NY]
  for t in [PY, CY]
    prem[t] = prem_w[t] + upr[t] - upr[t + 1]
  end
  ## calculate statistical parameters
  vol_prem = max(prem[PY], prem[CY])
  np_fac =
    df_lobs_prem_risk[df_lobs_prem_risk[:abbrev] .== name,
                      :re_np_prem][1,1]
  vol_prem_vc =
    df_lobs_prem_risk[df_lobs_prem_risk[:abbrev] .== name,
                      :prem_gross_vc][1,1]
  vol_prem_sd = vol_prem * vol_prem_vc * np_fac
  β =  β_raw_vec / sum(β_raw_vec)
  vol_res =
    lob[1, :res_undisc_py] * (β ⋅ cumprod(1 ./ (1 .+ rf)))
  vol_res_vc =
    df_lobs_prem_risk[df_lobs_prem_risk[:abbrev] .== name,
                      :res_net_vc][1,1]
  vol_res_sd = vol_res * vol_res_vc
  vol_prem_res_sd =
    sqrt(([vol_prem_sd, vol_res_sd]' *
          corr_prem_res *
            [vol_prem_sd, vol_res_sd])[1,1])
  return NLLob(name,
               lob[1, :index],
               prem_gross_w, prem_w,
               prem_gross_cy, prem,
               upr_gross, upr,
               re_prop_q, β,
               vol_prem,  vol_prem_vc,
               vol_res,
               vol_prem_res_sd)
end

"""
`premrestotalsd(nllobs::Vector{NLLob}, corr_lob::Matrix{Float64})`

Calculates the standard deviation of the total premium reserve
"""
function premrestotalsd(nllobs::Vector{NLLob},
                        corr_lob::Matrix{Float64})
  indices = Array(Int, 0)
  prem_res_sd = Array(Float64, 0)
  for i in 1:length(nllobs)
    push!(indices, nllobs[i].index)
    push!(prem_res_sd, nllobs[i].vol_prem_res_sd)
  end
  prem_res_total_sd =
    sqrt(prem_res_sd' * corr_lob[indices, indices] * prem_res_sd)
  return prem_res_total_sd
end


text:
using Test
using EliminateGraphs

@testset "Vertex set" begin
    @test Vertex(3) ∪ NeighborCover(4) isa UnionOf{Vertex, NeighborCover}
end


text:
include("IOM/src/models/EMOM/Topography.jl")

using Formatting
using NCDatasets
using ArgParse
using JSON

function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--zdomain-file"
            help = "NetCDF file containing z_w information"
            arg_type = String
            required = true

        "--topo-file"
            help = "Topography file."
            arg_type = String
            default = ""

        "--forcing-file"
            help = "Annual forcing file. It should contain: TAUX, TAUY, SWFLX, NSWFLX, VSFLX"
            arg_type = String
            default = ""

        "--forcing-file-varnames"
            help = "Variables that should be checked (comma separated). Default: TAUX,TAUY,SWFLX,NSWFLX,VSFLX,TEMP,SALT,QFLX_TEMP,QFLX_SALT. Unexisted variables will be skipped."
            arg_type = String
            default = "TAUX,TAUY,SWFLX,NSWFLX,VSFLX,TEMP,SALT,QFLX_TEMP,QFLX_SALT"

        "--forcing-file-xdim"
            help = "Varname of x dimension in the forcing file"
            arg_type = String
            default = "Nx"

        "--forcing-file-ydim"
            help = "Varname of y dimension in the forcing file"
            arg_type = String
            default = "Ny"

        "--forcing-file-zdim"
            help = "Varname of z dimension in the forcing file"
            arg_type = String
            default = "Nz"

        "--forcing-file-tdim"
            help = "Varname of time dimension in the forcing file"
            arg_type = String
            default = "time"

    end

    return parse_args(s)
end

parsed = parse_commandline()

print(json(parsed,4))

println("This program is used to see if")
println("1. Surface mask in the `domain-file` matches the surface mask in the topography created by `topo-file` and `zdomain-file`.")
println("2. If `forcing-file` is given, the program checks if 3D mask (T-grid) implied by `topo-file` is consistent with `forcing-file` (i.e. if mask[k,i,j] == 1, then forcing[k,i,j] should not be missing value).")

Dataset(parsed["domain-file"], "r") do ds
    global Nx = ds.dim["ni"]
    global Ny = ds.dim["nj"]
    global mask_sT = ds["mask"][:]

    if length(size(mask_sT)) == 2
        mask_sT = reshape(mask_sT, 1, size(mask_sT)...)
    end
end

Dataset(parsed["zdomain-file"], "r") do ds
    global z_w = ds["z_w"][:]
    global Nz = length(z_w) - 1
end

Dataset(parsed["topo-file"], "r") do ds
    global Nz_bot = ds["Nz_bot"][:]
end

topo = Topography(
    Nz_bot, Nx, Ny, z_w
)

flag_check1_pass = true
flag_check2_pass = true

print("### Check 1 : mask in domain-file should match topo.mask_T[1, :, :]. Answer: ")
if all(mask_sT .== topo.mask_T[1:1, :, :])
    println("Yes. It is consistent")
    flag_check1_pass = true
else
    println("*** Error: No, they are not. Please check. ***")
end

if parsed["forcing-file"] != ""
    println("`forcing-file` is given.")
    print("### Check 2 : If topo.mask_T[k,i,j] == 1 then forcing[k,i,j] should be finite. Answer: ")

    mask_sT_idx = mask_sT .== 1.0
    mask_T_idx  = topo.mask_T  .== 1.0

    Dataset(parsed["forcing-file"], "r") do ds

        Nt = ds.dim[parsed["forcing-file-tdim"]]

        for varname in split(parsed["forcing-file-varnames"], ",")
            if haskey(ds, varname)
                print(format("Check variable `{:s}`... ", varname))
                v = ds[varname]

                sdim = 0

                if dimnames(v) == (parsed["forcing-file-xdim"], parsed["forcing-file-ydim"], parsed["forcing-file-zdim"],  parsed["forcing-file-tdim"])
                    sdim = 3
                elseif dimnames(v) == (parsed["forcing-file-xdim"], parsed["forcing-file-ydim"], parsed["forcing-file-tdim"])
                    sdim = 2
                else
                    println("*** Error : Unrecognized dimensions: ", dimnames(v), ". Please check. ***")
                    flag_check2_pass = false
                    continue
                end

                if sdim == 3
                    mask_idx_used = mask_T_idx
                    mask_used = topo.mask_T
                    spatial_slice = (Colon(), Colon(), Colon())
                    reshape_dim   = (Nx, Ny, Nz) # reshape before permute
                elseif sdim == 2
                    mask_idx_used = mask_sT_idx
                    mask_used = mask_sT
                    spatial_slice = (Colon(), Colon())
                    reshape_dim   = (Nx, Ny, 1) # reshape before permute
                else
                    throw(ErrorException("Unknown situation. sdim = " * string(sdim)))
                end

                mask_used_cnt = sum(mask_used)

                for t = 1:Nt
                    var  = permutedims(reshape(nomissing(ds[varname][spatial_slice..., t],  NaN), reshape_dim...), [3, 1, 2])
                    var_isfinite_cnt = sum( isfinite.(var) )

                    if all(isfinite.(var[mask_idx_used]))
                        println(format("Pass"))
                        if var_isfinite_cnt != mask_used_cnt
                            println("Some of the data might not be used. (Forcing data has it but will not be used).")
                        else
                            println("All the data match the mask completely.")
                        end
                    else
                        println(format("*** Error: Variable `{:s}` is having missing values where mask is not zero. Please check. ***", varname))
                        println("var_isfinite_cnt = ", var_isfinite_cnt, "; mask_used_cnt = ", mask_used_cnt)
                        flag_check2_pass = false
                    end
                end

            end
        end
    end
end

println("### Summary: ")

println("Check 1: ", flag_check1_pass)

if parsed["forcing-file"] != ""
    println("Check 2: ", flag_check2_pass)
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
"""
    idct2!(x, fft_type::FFTType = Radix2Type())

Inverse discrete cosine transform of type II.

X = IDCT2(Y) inverts the DCT2 transform, returning the
original vector if Y was obtained using Y = DCT2(X). It is based on
the staggered-grid definition
                          N
   X(j) = 2/N*( Y(1)/2 +  ∑ Y(k) cos (π*(k-1)*(j-1/2)/N) )
                         k=2

X = IDCT2(Y,N) pads or truncates the vector Y to length N
before transforming.

If Y is a matrix, the IDCT2 operation is applied to
each column.

                          N
   X(j) = 2/N*( Y(1)/2 +  ∑ Y(k) cos (π*(k-1)*(j-1/2)/N) )
                         k=2
"""
function idct2!(x, fft_type=Radix2Type())
    x .= idct2(x, fft_type)
end
function idct2(a, fft_type::FFTType=Radix2Type())
    FT = eltype(a)
    if min(size(a)...)==1
        if size(a,2)>1
            do_trans = true;
        else
            do_trans = false;
        end
        a = a[:];
    else
        do_trans = false;
    end
    n = size(a,1);
    m = size(a,2);
    # Pad or truncate a if necessary
    if size(a,1)<n
      aa = zeros(n,m);
      aa[1:size(a,1),:] .= a;
    else
      aa = a[1:n,:];
    end
    y = Array{Complex{FT},2}(undef, 2*n, m)
    fill!(y, 0)
    y[1:n,:] .= aa;
    i = Complex(0,1);
    e = 0.5 * exp.(-i*FT(0.5) * π * (0:n-1) ./ n);
    for l = 1:m
        y[1:n,l] .= y[1:n,l] .* e;
    end
    y[n+2:2*n,:] .= reverse(conj.(y[2:n,:]); dims=1);
    yy = fft(y, 1);
    b = (2/n) * yy[1:n,:];
    if isreal(a)
        b = real.(b)
    end
    if do_trans
        b = b'
    end
    return b
end


text:
function ci_tapir(b::Cthulhu.Bookmark)
    ci, = code_typed(b)
    return Core.Compiler.lower_tapir(b.mi, ci)
end

function ircode_tapir(b::Cthulhu.Bookmark)
    ci, = code_typed(b)
    return Core.Compiler.lower_tapir_to_ircode(b.mi, ci)
end


text:
using DSGE
using BenchmarkTools, JLD

path = dirname(@__FILE__)
filepath = "$path/../../reference/BondLabor/jacobian.jld"

# Benchmark current function
m = BondLabor()
steadystate!(m)
trial = @benchmark DSGE.jacobian($m) gcsample = true

# # Optionally over-write the existing reference trial
# write_ref_trial(trial, "jacobian", filepath = filepath)

# Compute time differential
print_all_benchmarks(trial, filepath, "jacobian")


text:
# ---
# title: 1534. Count Good Triplets
# id: problem1534
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/count-good-triplets/description/>
# hidden: true
# ---
#
# Given an array of integers `arr`, and three integers `a`, `b` and `c`. You
# need to find the number of good triplets.
#
# A triplet `(arr[i], arr[j], arr[k])` is **good** if the following conditions
# are true:
#
#   * `0 <= i < j < k < arr.length`
#   * `|arr[i] - arr[j]| <= a`
#   * `|arr[j] - arr[k]| <= b`
#   * `|arr[i] - arr[k]| <= c`
#
# Where `|x|` denotes the absolute value of `x`.
#
# Return _the number of good triplets_.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [3,0,1,1,9,7], a = 7, b = 2, c = 3
#     Output: 4
#     Explanation:  There are 4 good triplets: [(3,0,1), (3,0,1), (3,1,1), (0,1,1)].
#
#
# **Example 2:**
#
#
#
#     Input: arr = [1,1,2,2,3], a = 0, b = 0, c = 1
#     Output: 0
#     Explanation: No triplet satisfies all conditions.
#
#
#
#
# **Constraints:**
#
#   * `3 <= arr.length <= 100`
#   * `0 <= arr[i] <= 1000`
#   * `0 <= a, b, c <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end