text:
# This file contains functions for the Profit Efficiency DEA model
"""
    ProfitDEAModel
An data structure representing a profit DEA model.
"""
struct ProfitDEAModel <: AbstractProfitDEAModel
    n::Int64
    m::Int64
    s::Int64
    Gx::Symbol
    Gy::Symbol
    monetary::Bool
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    lambda::SparseMatrixCSC{Float64, Int64}
    techeff::Vector
    alloceff::Vector
    normalization::Vector
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    deaprofit(X, Y, W, P; Gx, Gy)
Compute profit efficiency using data envelopment analysis model for
inputs `X`, outputs `Y`, price of inputs `W`, and price of outputs `P`.

# Direction specification:

The directions `Gx` and `Gy` can be one of the following symbols.
- `:Zeros`: use zeros.
- `:Ones`: use ones.
- `:Observed`: use observed values.
- `:Mean`: use column means.
- `:Monetary`: use direction so that profit inefficiency is expressed in monetary values.

Alternatively, a vector or matrix with the desired directions can be supplied.

# Optional Arguments
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [1 1; 1 1; 0.75 1.5; 0.5 2; 0.5 2; 2 2; 2.75 3.5; 1.375 1.75];

julia> Y = [1 11; 5 3; 5 5; 2 9; 4 5; 4 2; 3 3; 4.5 3.5];

julia> P = [2 1; 2 1; 2 1; 2 1; 2 1; 2 1; 2 1; 2 1];

julia> W = [2 1; 2 1; 2 1; 2 1; 2 1; 2 1; 2 1; 2 1];

julia> deaprofit(X, Y, W, P, Gx = :Monetary, Gy = :Monetary)
Profit DEA Model
DMUs = 8; Inputs = 2; Outputs = 2
Returns to Scale = VRS
Gx = Monetary; Gy = Monetary
─────────────────────────────────────
   Profit     Technical    Allocative
─────────────────────────────────────
1     2.0   0.0           2.0
2     2.0  -5.41234e-16   2.0
3     0.0   0.0           0.0
4     2.0   0.0           2.0
5     2.0   0.0           2.0
6     8.0   6.0           2.0
7    12.0  12.0          -1.77636e-15
8     4.0   3.0           1.0
─────────────────────────────────────
```
"""
function deaprofit(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    W::Union{Matrix,Vector}, P::Union{Matrix,Vector};
    Gx::Union{Symbol,Matrix,Vector}, Gy::Union{Symbol,Matrix,Vector},
    monetary::Bool = false,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::ProfitDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    nw, mw = size(W, 1), size(W, 2)
    np, sp = size(P, 1), size(P, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nw != nx
        throw(DimensionMismatch("number of rows in W and X ($nw, $nx) are not equal"));
    end
    if np != ny
        throw(DimensionMismatch("number of rows in P and Y ($np, $ny) are not equal"));
    end
    if mw != m
        throw(DimensionMismatch("number of columns in W and X ($mw, $m) are not equal"));
    end
    if sp != s
        throw(DimensionMismatch("number of columns in P and Y ($sp, $s) are not equal"));
    end

    # Build or get user directions
    if typeof(Gx) == Symbol
        Gxsym = Gx

        if Gx == :Zeros
            Gx = zeros(size(X))
        elseif Gx == :Ones
            Gx = ones(size(X))
        elseif Gx == :Observed
            Gx = X
        elseif Gx == :Mean
            Gx = repeat(mean(X, dims = 1), size(X, 1))
        elseif Gx == :Monetary
            GxGydollar = 1 ./ (sum(P, dims = 2) + sum(W, dims = 2));
            Gx = repeat(GxGydollar, 1, m);
        else
            throw(ArgumentError("Invalid `Gx`"));
        end

    else
        Gxsym = :Custom
    end

    if typeof(Gy) == Symbol
        Gysym = Gy

        if Gy == :Zeros
            Gy = zeros(size(Y))
        elseif Gy == :Ones
            Gy = ones(size(Y))
        elseif Gy == :Observed
            Gy = Y
        elseif Gy == :Mean
            Gy = repeat(mean(Y, dims = 1), size(Y, 1))
        elseif Gy == :Monetary
            GxGydollar = 1 ./ (sum(P, dims = 2) + sum(W, dims = 2));
            Gy = repeat(GxGydollar, 1, s);
        else
            throw(ArgumentError("Invalid `Gy`"));
        end

    else
        Gysym = :Custom
    end

    if (size(Gx, 1) != size(X, 1)) | (size(Gx, 2) != size(X, 2))
        throw(DimensionMismatch("size of Gx and X ($(size(Gx)), $(size(X))) are not equal"));
    end
    if (size(Gy, 1) != size(Y, 1)) | (size(Gy, 2) != size(Y, 2))
        throw(DimensionMismatch("size of Gy and Y ($(size(Gy)), $(size(Y))) are not equal"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    # Get maximum profit targets and lambdas
    n = nx

    Xtarget, Ytarget, plambdaeff = deamaxprofit(X, Y, W, P, optimizer = optimizer)

    # Profit, technical and allocative efficiency
    maxprofit = sum(P .* Ytarget, dims = 2) .- sum(W .* Xtarget, dims = 2)

    pefficiency  = vec(maxprofit .- ( sum(P .* Y, dims = 2) .- sum(W .* X, dims = 2)))
    normalization = vec(sum(P .* Gy, dims = 2) .+ sum(W .* Gx, dims = 2))
    techefficiency = efficiency(deaddf(X, Y, Gx = Gx, Gy = Gy, rts = :VRS, slack = false, optimizer = optimizer))

    if monetary
        techefficiency = techefficiency .* normalization
    else
        pefficiency = pefficiency ./ normalization
    end

    allocefficiency = pefficiency - techefficiency

    return ProfitDEAModel(n, m, s, Gxsym, Gysym, monetary, names, pefficiency, plambdaeff, techefficiency, allocefficiency, normalization, Xtarget, Ytarget)

end

function Base.show(io::IO, x::ProfitDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    dmunames = names(x)

    eff = efficiency(x)
    techeff = efficiency(x, :Technical)
    alloceff = efficiency(x, :Allocative)

    if !compact
        print(io, "Profit DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Returns to Scale = VRS")
        print(io, "\n")
        print(io, "Gx = ", string(x.Gx), "; Gy = ", string(x.Gy))
        print(io, "\n")
        show(io, CoefTable(hcat(eff, techeff, alloceff), ["Profit", "Technical", "Allocative"], dmunames))
    end

end


text:
checkbandwidths(N, M, l::AbstractVector{Int}, u::AbstractVector{Int}) =
    M == 1 || (length(u) == M && length(l) == M) ||
    throw(DimensionMismatch("For a matrix of $(N)×$(M) blocks, $(M) lower and upper column bandwidths are required"))

#### Routines for BandedSizes
function bb_blockstarts(b_size, l::AbstractVector{Int}, u::AbstractVector{Int})
    N,M = nblocks(b_size)
    L,U = maximum(l), maximum(u)
    b_start = BandedMatrix{Int}(undef, (N, M), (L, U))
    -L > U && return b_start

    checkbandwidths(N, M, l, u)

    ind_shift = 0
    for J = 1:M
        KR = max(1,J-u[J]):min(J+l[J],N)
        if !isempty(KR)
            b_start[KR,J] .= ind_shift .+ view(cumulsizes(b_size,1),KR) .- cumulsizes(b_size,1)[KR[1]] .+ 1

            num_rows = cumulsizes(b_size,1,KR[end]+1)-cumulsizes(b_size,1,KR[1])
            num_cols = blocksize(b_size, 2, J)
            ind_shift += num_rows*num_cols
        end
    end

    b_start
end

bb_blockstarts(b_size, l::Integer, u::Integer) = bb_blockstarts(b_size, Fill(l, nblocks(b_size,2)), Fill(u, nblocks(b_size,2)))

function bb_blockstrides(b_size, l::AbstractVector{Int}, u::AbstractVector{Int})
    N, M = nblocks(b_size)
    L,U = maximum(l), maximum(u)
    checkbandwidths(N, M, l, u)
    b_strides = Vector{Int}(undef, M)
    for J=1:M
        KR = max(1,J-u[J]):min(J+l[J],N)
        if !isempty(KR)
            b_strides[J] = cumulsizes(b_size,1,KR[end]+1)-cumulsizes(b_size,1,KR[1])
        else
            b_strides[J] = 0
        end
    end
    b_strides
end

bb_blockstrides(b_size, l::Integer, u::Integer) = bb_blockstrides(b_size, Fill(l, nblocks(b_size,1)), Fill(u, nblocks(b_size,2)))

struct BlockSkylineSizes{BS<:AbstractBlockSizes{2}, LL<:AbstractVector{Int}, UU<:AbstractVector{Int}, BStarts, BStrides} <: AbstractBlockSizes{2}
    block_sizes::BS
    block_starts::BStarts # gives index where the blocks start, usually a BandedMatrix{Int}
    block_strides::BStrides # gives stride to next block for J-th column, usually a Vector{Int}
    l::LL
    u::UU
end

BlockSkylineSizes(b_size::BlockSizes{2}, l::AbstractVector{Int}, u::AbstractVector{Int}) =
    BlockSkylineSizes(b_size, bb_blockstarts(b_size, l, u), bb_blockstrides(b_size, l, u), l, u)

BlockSkylineSizes(rows::AbstractVector{Int}, cols::AbstractVector{Int}, l::AbstractVector{Int}, u::AbstractVector{Int}) =
    BlockSkylineSizes(BlockSizes(rows,cols), l, u)

const BlockBandedSizes = BlockSkylineSizes{DefaultBlockSizes{2}, Fill{Int,1,Tuple{OneTo{Int}}}, Fill{Int,1,Tuple{OneTo{Int}}},
                                            BandedMatrix{Int,Matrix{Int},OneTo{Int}}, Vector{Int}}

BlockBandedSizes(b_size::BlockSizes{2}, l::Int, u::Int) =
    BlockSkylineSizes(b_size, Fill(l, nblocks(b_size,2)), Fill(u, nblocks(b_size,2)))
BlockBandedSizes(rows::AbstractVector{Int}, cols::AbstractVector{Int}, l::Int, u::Int) =
    BlockSkylineSizes(rows, cols, Fill(l, length(cols)), Fill(u, length(cols)))

colblockbandwidths(bs::BlockSkylineSizes) = (bs.l, bs.u)

blockstart(block_sizes::BlockSkylineSizes, K, J) = block_sizes.block_starts[K,J]
blockstride(block_sizes::BlockSkylineSizes, J) = block_sizes.block_strides[J]

blockstart(A::AbstractMatrix, K, J) = blockstart(blocksizes(A),K,J)
blockstride(A::AbstractMatrix, J) = blockstride(blocksizes(A),J)

blockbandwidths(bs::BlockSkylineSizes) = maximum.(colblockbandwidths(bs))
blockbandwidth(bs::BlockSkylineSizes, i::Int) = blockbandwidths(bs)[i]

==(A::BlockSkylineSizes, B::BlockSkylineSizes) = A.block_sizes == B.block_sizes && A.block_starts == B.block_starts &&
    A.l == B.l && A.u == B.u

cumulsizes(B::BlockSkylineSizes) = cumulsizes(B.block_sizes)

colrange(B::BlockSkylineSizes, J::Integer) = max(1, J-B.u[J]):min(nblocks(B.block_sizes,1), J+B.l[J])

function bb_numentries(B::BlockSkylineSizes)
    b_size = B.block_sizes
    numentries = 0
    N = nblocks(b_size,1)
    for J = 1:nblocks(b_size,2),
        KR = colrange(B, J)
        num_rows = cumulsizes(b_size,1,KR[end]+1) - cumulsizes(b_size,1,KR[1])
        num_cols = blocksize(b_size, 2, J)
        numentries += num_rows*num_cols
    end
    numentries
end

function _BandedBlockMatrix end

#  A block matrix where only the bands are nonzero
#   isomorphic to BandedMatrix{Matrix{T}}
struct BlockSkylineMatrix{T, DATA<:AbstractVector{T}, BS<:AbstractBlockSizes{2}} <: AbstractBlockBandedMatrix{T}
    data::DATA
    block_sizes::BS

    global function _BlockSkylineMatrix(data::DATA, block_sizes::BS) where {T,DATA<:AbstractVector{T}, BS<:AbstractBlockSizes{2}}
        new{T,DATA,BS}(data, block_sizes)
    end
end

const BlockBandedMatrix{T} = BlockSkylineMatrix{T, Vector{T}, BlockBandedSizes}

# Auxiliary outer constructors
@inline _BlockBandedMatrix(data::AbstractVector, bs::BlockBandedSizes) =
    _BlockSkylineMatrix(data, bs)

@inline _BlockBandedMatrix(data::AbstractVector, (kr,jr)::NTuple{2, AbstractVector{Int}}, (l,u)::NTuple{2, Int}) =
    _BlockBandedMatrix(data, BlockBandedSizes(kr,jr, l,u))

@inline BlockSkylineMatrix{T}(::UndefInitializer, block_sizes::BlockSkylineSizes) where T =
    _BlockSkylineMatrix(Vector{T}(undef, bb_numentries(block_sizes)), block_sizes)

@inline BlockBandedMatrix{T}(::UndefInitializer, block_sizes::BlockBandedSizes) where T =
    _BlockSkylineMatrix(Vector{T}(undef, bb_numentries(block_sizes)), block_sizes)

@inline BlockBandedMatrix{T}(::UndefInitializer, block_sizes::BlockSizes, (l,u)::NTuple{2, Int}) where T =
    BlockSkylineMatrix{T}(undef, BlockBandedSizes(block_sizes,l,u))

@inline BlockSkylineMatrix{T}(::UndefInitializer, block_sizes::BlockSizes, (l,u)::NTuple{2, AbstractVector{Int}}) where T =
    BlockSkylineMatrix{T}(undef, BlockSkylineSizes(block_sizes,l,u))

"""
    BlockSkylineMatrix{T,LL,UU}(undef, (rows, cols), (l::LL, u::UU))

returns an undef `sum(rows)`×`sum(cols)` block-banded matrix `A`
of type `T` with block-bandwidths `(l,u)` and where `A[Block(K,J)]`
is a `Matrix{T}` of size `rows[K]`×`cols[J]`.

`(l,u)` may be integers for constant bandwidths or integer vectors of
lengths `rows` and `cols`, respectively, for ragged bands.
"""
BlockSkylineMatrix

@inline BlockBandedMatrix{T}(::UndefInitializer, dims::NTuple{2, AbstractVector{Int}}, lu::NTuple{2, Int}) where T =
    BlockSkylineMatrix{T}(undef, BlockBandedSizes(dims..., lu...))

@inline BlockSkylineMatrix{T}(::UndefInitializer, dims::NTuple{2, AbstractVector{Int}}, lu::NTuple{2, AbstractVector{Int}}) where T =
    BlockSkylineMatrix{T}(undef, BlockSkylineSizes(dims..., lu...))

function BlockSkylineMatrix{T}(A::AbstractMatrix, block_sizes::BlockSkylineSizes) where T
    ret = BlockSkylineMatrix(Zeros{T}(size(A)), block_sizes)
    for J = Block.(1:nblocks(ret, 2)), K = blockcolrange(ret, Int(J))
        kr, jr = globalrange(block_sizes, (Int(K), Int(J)))
        view(ret, K, J) .= view(A, kr, jr)
    end
    ret
end

function BlockSkylineMatrix{T}(A::AbstractBlockBandedMatrix, block_sizes::BlockSkylineSizes) where T
    ret = BlockSkylineMatrix(Zeros{T}(size(A)), block_sizes)
    block_sizes == blocksizes(A) || throw(ArgumentError())
    for J = Block.(1:nblocks(ret, 2)), K = blockcolrange(ret, Int(J))
        view(ret, K, J) .= view(A, K, J)
    end
    ret
end

BlockBandedMatrix{T}(A::AbstractMatrix, block_sizes::BlockBandedSizes) where T =
    BlockSkylineMatrix{T}(A, block_sizes)

##
# Special cases
##

function BlockSkylineMatrix{T}(Z::Zeros, block_sizes::BlockSkylineSizes) where T
    if size(Z) ≠ size(block_sizes)
        throw(DimensionMismatch("Size of input $(size(Z)) must be consistent with $(size(block_sizes))"))
    end
    _BlockSkylineMatrix(zeros(T, bb_numentries(block_sizes)), block_sizes)
end

function BlockSkylineMatrix{T}(E::Eye, block_sizes::BlockSkylineSizes) where T
    if size(E) ≠ size(block_sizes)
        throw(DimensionMismatch("Size of input $(size(E)) must be consistent with $(sum.(dims))"))
    end
    ret = BlockSkylineMatrix(Zeros{T}(size(E)), block_sizes)
    ret[diagind(ret)] .= one(T)
    ret
end

function BlockSkylineMatrix{T}(A::UniformScaling, block_sizes::BlockSkylineSizes) where T
    ret = BlockSkylineMatrix(Zeros{T}(size(block_sizes)), block_sizes)
    ret[diagind(ret)] .= convert(T, A.λ)
    ret
end

BlockSkylineMatrix(A::Union{AbstractMatrix,UniformScaling},
                        block_sizes::BlockSkylineSizes) = BlockSkylineMatrix{eltype(A)}(A, block_sizes)
BlockBandedMatrix(A::Union{AbstractMatrix,UniformScaling},
                        block_sizes::BlockBandedSizes) = BlockBandedMatrix{eltype(A)}(A, block_sizes)

BlockSkylineMatrix{T}(A::Union{AbstractMatrix,UniformScaling},
                           dims::NTuple{2,AbstractVector{Int}}, lu::NTuple{2,AbstractVector{Int}}) where T =
                               BlockSkylineMatrix{T}(A, BlockSkylineSizes(dims..., lu...))

BlockBandedMatrix{T}(A::Union{AbstractMatrix,UniformScaling},
                           dims::NTuple{2,AbstractVector{Int}}, lu::NTuple{2,Int}) where T =
                               BlockSkylineMatrix{T}(A, BlockBandedSizes(dims..., lu...))

BlockSkylineMatrix(A::Union{AbstractMatrix,UniformScaling},
                        dims::NTuple{2, AbstractVector{Int}},
                        lu::NTuple{2,AbstractVector{Int}}) = BlockSkylineMatrix{eltype(A)}(A, dims, lu)
BlockBandedMatrix(A::Union{AbstractMatrix,UniformScaling},
                        dims::NTuple{2, AbstractVector{Int}},
                        lu::NTuple{2,Int}) = BlockBandedMatrix{eltype(A)}(A, dims, lu)

BlockBandedMatrix(A::BlockBandedMatrix, lu::NTuple{2,Int}) = BlockBandedMatrix(A, BlockBandedSizes(blocksizes(A).block_sizes, lu...))

function convert(::Type{BlockSkylineMatrix}, A::AbstractMatrix)
    @assert isblockbanded(A)
    block_sizes = convert(BlockSkylineSizes, A.block_sizes)

    ret = BlockSkylineMatrix{eltype(A)}(undef, block_sizes)
    for J = Block.(1:nblocks(ret, 2)), K = blockcolrange(ret, Int(J))
        view(ret, K, J) .= view(A, K, J)
    end
    ret
end

function convert(::Type{BlockBandedMatrix}, A::AbstractMatrix)
    @assert isblockbanded(A)
    convert(BlockSkylineMatrix, A)
end

BlockSkylineMatrix(A::AbstractMatrix) = convert(BlockSkylineMatrix, A)
BlockBandedMatrix(A::AbstractMatrix) = convert(BlockBandedMatrix, A)

similar(A::BlockSkylineMatrix, T::Type=eltype(A), bs::BlockSkylineSizes=blocksizes(A)) =
    BlockSkylineMatrix{T}(undef, bs)

################################
# BlockSkylineMatrix Interface #
################################

MemoryLayout(::Type{<:BlockSkylineMatrix}) = BlockBandedColumnMajor()
colblockbandwidths(A::BlockSkylineMatrix) = (A.block_sizes.l, A.block_sizes.u)
blockbandwidths(A::BlockSkylineMatrix) = maximum.(colblockbandwidths(A))
BroadcastStyle(::Type{<:BlockSkylineMatrix}) = BlockSkylineStyle()
BroadcastStyle(::Type{<:BlockBandedMatrix}) = BlockBandedStyle()

################################
# AbstractBlockArray Interface #
################################

@inline blocksizes(block_array::BlockSkylineMatrix) = block_array.block_sizes

zeroblock(A::BlockSkylineMatrix, K::Int, J::Int) =
    Matrix(Zeros{eltype(A)}(blocksize(A, (K, J))))

@inline function getblock(A::BlockSkylineMatrix, K::Int, J::Int)
    @boundscheck blockcheckbounds(A, K, J)
    if -A.block_sizes.l[J] ≤ J - K ≤ A.block_sizes.u[J]
        convert(Matrix, view(A, Block(K, J)))
    else
        zeroblock(A, K, J)
    end
end

# @inline function Base.getindex(block_arr::BlockArray{T,N}, blockindex::BlockIndex{N}) where {T,N}
#     @boundscheck checkbounds(block_arr.blocks, blockindex.I...)
#     @inbounds block = block_arr.blocks[blockindex.I...]
#     @boundscheck checkbounds(block, blockindex.α...)
#     @inbounds v = block[blockindex.α...]
#     return v
# end

###########################
# AbstractArray Interface #
###########################

# @inline function Base.similar(block_array::BlockSkylineMatrix{T}, ::Type{T2}) where {T,N,T2}
#     BlockArray(similar(block_array.blocks, Array{T2, N}), copy(block_array.block_sizes))
# end

Base.size(arr::BlockSkylineMatrix) =
    @inbounds return (cumulsizes(arr.block_sizes.block_sizes,1)[end] - 1, cumulsizes(arr.block_sizes.block_sizes,2)[end] - 1)

@inline function getindex(A::BlockSkylineMatrix, i::Int, j::Int)
    @boundscheck checkbounds(A, i, j)
    bi = global2blockindex(A.block_sizes, (i, j))
    @inbounds v = view(A, Block(bi.I))[bi.α...]
    return v
end

@inline function setindex!(A::BlockSkylineMatrix{T}, v, i::Int, j::Int) where T
    @boundscheck checkbounds(A, i, j)
    bi = global2blockindex(A.block_sizes, (i, j))
    V = view(A, Block(bi.I))
    @inbounds V[bi.α...] = convert(T, v)::T
    return v
end

## structured matrix methods ##
function Base.replace_in_print_matrix(A::BlockSkylineMatrix, i::Integer, j::Integer, s::AbstractString)
    bi = global2blockindex(A.block_sizes, (i, j))
    I,J = bi.I
    -A.block_sizes.l[J] ≤ J-I ≤ A.block_sizes.u[J] ? s : Base.replace_with_centered_mark(s)
end

############
# Indexing #
############

# function _check_setblock!(block_arr::BlockArray{T, N}, v, block::NTuple{N, Int}) where {T,N}
#     for i in 1:N
#         if size(v, i) != blocksize(block_arr.block_sizes, (i, block[i]))
#             throw(DimensionMismatch(string("tried to assign $(size(v)) array to ", blocksize(block_arr, block), " block")))
#         end
#     end
# end
#
#
@inline function setblock!(A::BlockSkylineMatrix, v, K::Int, J::Int)
    @boundscheck blockcheckbounds(A, K, J)

    V = view(A, Block(K), Block(J))
    V .= v
    return A
end

# @propagate_inbounds function Base.setindex!(block_array::BlockArray{T, N}, v, block_index::BlockIndex{N}) where {T,N}
#     getblock(block_array, block_index.I...)[block_index.α...] = v
# end

########
# Misc #
########

# @generated function Base.Array(block_array::BlockArray{T, N, R}) where {T,N,R}
#     # TODO: This will fail for empty block array
#     return quote
#         block_sizes = block_array.block_sizes
#         arr = similar(block_array.blocks[1], size(block_array)...)
#         @nloops $N i i->(1:nblocks(block_sizes, i)) begin
#             block_index = @ntuple $N i
#             indices = globalrange(block_sizes, block_index)
#             arr[indices...] = getblock(block_array, block_index...)
#         end
#
#         return arr
#     end
# end
#
# @generated function Base.copyto!(block_array::BlockArray{T, N, R}, arr::R) where {T,N,R <: AbstractArray}
#     return quote
#         block_sizes = block_array.block_sizes
#
#         @nloops $N i i->(1:nblocks(block_sizes, i)) begin
#             block_index = @ntuple $N i
#             indices = globalrange(block_sizes, block_index)
#             copyto!(getblock(block_array, block_index...), arr[indices...])
#         end
#
#         return block_array
#     end
# end
#
# function Base.fill!(block_array::BlockArray, v)
#     for block in block_array.blocks
#         fill!(block, v)
#     end
# end

##################
# BlockBandedBlock
#
#   views of the blocks satisfy the Matrix interface, and are memory-compatible
#   with StridedMatrix.
##################

const BlockBandedBlock{T,LL,UU} = SubArray{T,2,BlockSkylineMatrix{T,LL,UU},Tuple{BlockSlice1,BlockSlice1},false}

# gives the columns of parent(V).data that encode the block
blocks(V::BlockBandedBlock)::Tuple{Int,Int} = first(first(parentindices(V)).block.n),
first(last(parentindices(V)).block.n)

######################################
# Matrix interface  for Blocks #
######################################

MemoryLayout(::Type{<:BlockBandedBlock}) = ColumnMajor()

function Base.unsafe_convert(::Type{Ptr{T}}, V::BlockBandedBlock{T}) where T
    A = parent(V)
    K,J = blocks(V)
    Base.unsafe_convert(Ptr{T}, A.data) + sizeof(T)*(blockstart(A,K,J)-1)
end

strides(V::BlockBandedBlock) = (1,parent(V).block_sizes.block_strides[blocks(V)[2]])

@propagate_inbounds function getindex(V::BlockBandedBlock, k::Int, j::Int)
    @boundscheck checkbounds(V, k, j)
    A = parent(V)
    K,J = blocks(V)
    if -A.block_sizes.l[J] ≤ J-K ≤ A.block_sizes.u[J]
        b_start = blockstart(A,K,J)
        b_start == 0 && return zero(eltype(V))
        b_stride = blockstride(A,J)
        A.data[b_start + k-1 + (j-1)*b_stride ]
    else
        zero(eltype(V))
    end
end

@propagate_inbounds function setindex!(V::BlockBandedBlock, v, k::Int, j::Int)
    @boundscheck checkbounds(V, k, j)
    A = parent(V)
    K,J = blocks(V)
    if -A.block_sizes.l[J] ≤ J-K ≤ A.block_sizes.u[J]
        b_start = blockstart(A,K,J)
        # TODO: What to do if b_start == 0 ?
        b_stride = A.block_sizes.block_strides[J]
        A.data[b_start + k-1 + (j-1)*b_stride ] = v
    elseif iszero(v) # allow setindex for 0 data
        v
    else
        throw(BandError(A, J-K))
    end
end

ext: jl


text:
using LinearAlgebra, Random, Unicode, DelimitedFiles
using Combinatorics: binomial

@testset "bimatrix_generators.jl" begin

    @testset "blotto_game" begin
        h = 3
        t = 7
        rho = 0.5

        g = @inferred(blotto_game(h, t, rho))

        @testset "test_nb_actions" begin
            nb_actions = (factorial(t+h-1)/
                          (factorial(h-1)*factorial(t)))

            @test g.nums_actions == (nb_actions, nb_actions)
        end

        @testset "test_constant_diagonal" begin
            for i=1:2
                d = diag(g.players[i].payoff_array)
                @test length(unique(d)) == 1
            end
        end

        @testset "test_seed" begin
            g1 = blotto_game(MersenneTwister(32), h, t, rho)
            g2 = blotto_game(MersenneTwister(32), h, t, rho)
            for i in 1:2
                @test g1.players[i].payoff_array == g2.players[i].payoff_array
            end
        end

    end

    @testset "ranking_game" begin
        g = @inferred(ranking_game(100))
        p1_array = g.players[1].payoff_array
        p2_array = g.players[2].payoff_array

        @testset "test_weakly_decreasing_rowise_payoffs" begin
            @test all(p1_array[:, 1:(end-1)] - p1_array[:, 2:end] .>=0)
            @test all(p2_array[:, 1:(end-1)] - p2_array[:, 2:end] .>=0)
        end

        @testset "test_elements_first_row" begin
            @test sum(g[1, 1]) == 1
            possible_elements = [0, 1, 0.5]
            @test all([value in possible_elements for value in p1_array[1, :]])
            @test all([value in possible_elements for value in p2_array[1, :]])
        end

        @testset "test_seed" begin
            g1 = ranking_game(MersenneTwister(32), 20)
            g2 = ranking_game(MersenneTwister(32), 20)
            for i in 1:2
                @test g1.players[i].payoff_array == g2.players[i].payoff_array
            end
        end

    end

    @testset "sgc_game" begin
        k = 2
        n = 4*k - 1
        g = @inferred(sgc_game(k))

        s = "
            0.750 0.750 1.000 0.500 0.500 1.000 0.000 0.500 0.000 0.500 0.000
            0.500 0.000 0.500 0.500 1.000 0.750 0.750 1.000 0.500 0.000 0.500
            0.000 0.500 0.000 0.500 0.000 0.500 1.000 0.500 0.500 1.000 0.750
            0.750 0.000 0.500 0.000 0.500 0.000 0.500 0.000 0.500 0.500 0.000
            0.500 0.000 0.500 0.000 0.750 0.000 0.000 0.750 0.000 0.000 0.000
            0.000 0.500 0.000 0.500 0.000 0.500 0.000 0.000 0.750 0.750 0.000
            0.000 0.000 0.000 0.000 0.500 0.000 0.500 0.000 0.500 0.000 0.000
            0.000 0.000 0.000 0.750 0.000 0.000 0.750 0.500 0.000 0.500 0.000
            0.500 0.000 0.000 0.000 0.000 0.000 0.000 0.750 0.750 0.000"
        payoffs = readdlm(IOBuffer(Unicode.normalize(s, stripcc=true)))
        payoffs = reshape(payoffs, (2, n^2))
        payoff_matrices = [reshape(payoffs[i, :], (n, n)) for i in 1:2]
        payoff_matrices[2] = transpose(payoff_matrices[2])

        for i in 1:2
            @test g.players[i].payoff_array == payoff_matrices[i]
        end
    end

    @testset "tournament_game" begin
        n, k = 6, 4
        m = binomial(n, k)
        g = @inferred(tournament_game(n, k))

        @testset "test_size" begin
            @test g.nums_actions == (n, m)
        end

        @testset "test_payoff_values" begin
            possible_values = [0., 1.]
            for player in g.players
                @test issubset(player.payoff_array, possible_values)
            end

            max_num_dominated_subsets = sum([binomial(i, k) for i in k:(n-1)])
            @test sum(g.players[1].payoff_array) <= max_num_dominated_subsets
            @test all(sum(g.players[2].payoff_array, dims=2) .== k)
        end

        @testset "test_seed" begin
            seed = 0
            g1 = tournament_game(n, k; seed=seed)
            g2 = tournament_game(n, k; seed=seed)

            for i in 1:2
                @test g1.players[i].payoff_array == g2.players[i].payoff_array
            end
        end

        @testset "test_throws_argument_error_too_large_inputs" begin
            n, k = 100, 50
            @test_throws ArgumentError tournament_game(n, k)
        end
    end

    @testset "unit_vector_game" begin
        n = 100
        g = @inferred(unit_vector_game(n))

        @testset "test_size" begin
            @test g.nums_actions == (n, n)
        end

        @testset "test_payoff_values" begin
            @test all(sum(g.players[1].payoff_array, dims=1) .== 1.)
        end

        @testset "test_avoid_pure_nash" begin
            NEs = pure_nash(unit_vector_game(n, avoid_pure_nash=true), tol=0.)
            @test length(NEs) == 0
        end

        @testset "test_seed" begin
            seed = 0
            n = 100
            g1 = unit_vector_game(MersenneTwister(seed), n)
            g2 = unit_vector_game(MersenneTwister(seed), n)
            for i in 1:2
                @test g1.players[i].payoff_array == g2.players[i].payoff_array
            end
        end

        @testset "test_redraw" begin
            seed = 6
            rng = MersenneTwister(seed)
            n = 2
            g = unit_vector_game(rng, n, avoid_pure_nash=true)
            NEs = pure_nash(g, tol=0.)
            @test length(NEs) == 0
        end

        @testset "test_throws_argument_error_avoid_pure_nash_n_1" begin
            n = 1
            @test_throws ArgumentError unit_vector_game(n, avoid_pure_nash=true)
        end
    end

end


text:
using LightGraphs
using SimpleWeightedGraphs
using LinearAlgebra
using GraphPlot
using Colors
using Serialization
using Xtals
using ..ChemistryFeaturization.Utils.GraphBuilding
using ..ChemistryFeaturization.AbstractType: AbstractAtoms

# TO CONSIDER: store ref to featurization rather than the thing itself? Does this matter for any performance we care about?
"""
    AtomGraph

A type representing an atomic structure as a graph (`gr`).

# Fields
- `graph::SimpleWeightedGraph{<:Integer,<:Real}`: the graph representing the structure. See
  [`build_graph`](@ref) for more on generating the weights.
- `elements::Vector{String}`: list of elemental symbols corresponding to each node of the
  graph
- `laplacian::Matrix{<:Real}`: Normalized graph Laplacian matrix, stored to speed up
  convolution operations by avoiding recomputing it every pass.
- `id::String`: Optional, an identifier, e.g. to correspond with tags/labels of an imported
  dataset.
"""
mutable struct AtomGraph <: AbstractAtoms
    graph::SimpleWeightedGraph{<:Integer,<:Real}
    elements::Vector{String}
    laplacian::Matrix{<:Real} # wanted to use LightGraphs.LinAlg.NormalizedGraphLaplacian but seems this doesn't support weighted graphs?
    id::String # or maybe we let it be a number too?
end

# one without features or featurization initialized yet
function AtomGraph(
    graph::SimpleWeightedGraph{A,B},
    elements::Vector{String},
    id = "",
) where {B<:Real,A<:Integer}
    # check that elements is the right length
    num_atoms = size(graph)[1]
    @assert length(elements) == num_atoms "Element list length doesn't match graph size!"

    # this was previously B.(normalized_laplacian(graph)) - won't that potentially give rise to compatibility issues if B is a custom type?
    laplacian = normalized_laplacian(graph)
    AtomGraph(graph, elements, laplacian, id)
end

# initialize directly from adjacency matrix
AtomGraph(adj::Array{R}, elements::Vector{String}, id = "") where {R<:Real} =
    AtomGraph(SimpleWeightedGraph(adj), elements, id)

"""
    AtomGraph(input_file_path, id = splitext(input_file_path)[begin]; output_file_path = nothing, overwrite_file = false, use_voronoi = false, cutoff_radius = 8.0, max_num_nbr = 12, dist_decay_func = inverse_square)

Construct an AtomGraph object from a structure file.

# Required Arguments
- `input_file_path::String`: path to file containing structure (must be readable by ASE.io.read)

# Optional Arguments
- `id::String`: ID associated with structure (e.g. identifier from online database). Defaults to name of input file if undefined.
- `output_file_path = nothing`: If provided, structure will be serialized to file at this location
- `overwrite_file::Bool = false`: whether to overwrite an existing file at `output_file_path`
- `use_voronoi::Bool = false`: Whether to build neighbor lists using Voronoi decompositions
- `cutoff_radius::Real = 8.0`: If not using Voronoi neighbor lists, longest allowable distance to a neighbor, in Angstroms
- `max_num_nbr::Integer = 12`: If not using Voronoi neighbor lists, largest allowable number of neighbors
- `dist_decay_func = inverse_square`: Function by which to assign edge weights according to distance between neighbors

# Note
`max_num_nbr` is a "soft" limit – if multiple neighbors are at the same distance, the full neighbor list may be longer.
"""
function AtomGraph(
    input_file_path::String,
    id::String = splitext(input_file_path)[begin];
    output_file_path::Union{String,Nothing} = nothing,
    overwrite_file::Bool = false,
    use_voronoi::Bool = false,
    cutoff_radius::Real = 8.0,
    max_num_nbr::Integer = 12,
    dist_decay_func::Function = inverse_square,
)

    local ag

    if !isfile(input_file_path)
        @warn "$input_file_path does not exist. Cannot build graph from a non-existent file."
        return missing
    end

    if splitext(input_file_path)[end] == ".jls" # deserialize
        ag = deserialize(input_file_path)
        ag.id = id

    else # try actually building the graph
        try
            adj_mat, elements = build_graph(
                input_file_path,
                use_voronoi = use_voronoi,
                cutoff_radius = cutoff_radius,
                max_num_nbr = max_num_nbr,
                dist_decay_func = dist_decay_func,
            )
            ag = AtomGraph(adj_mat, elements, id)
        catch
            @warn "Unable to build graph for $input_file_path"
            return missing
        end
    end

    to_serialize = !isnothing(output_file_path)
    if to_serialize
        if isfile(output_file_path) && !(overwrite_file)
            @info "Output file already exists, and `overwrite_file` is set to false.\nIf you want to overwrite the existing graph, set `overwrite=true`, or remove the existing file and retry."
        else
            serialize(output_file_path, ag)
        end
    end

    return ag
end

"""
    AtomGraph(crys::Crystal; id="", cutoff_radius = 8.0, max_num_nbr = 12, dist_decay_func = inverse_square)

Construct an AtomGraph object from a Crystal object (defined in Xtals.jl). For now, only supports cutoff-based graph building.

# Required Arguments
- `crys::Crystal`: Crystal from which to build the graph

# Optional Arguments
- `id::String`: ID associated with structure (e.g. identifier from online database). Defaults to the empty string.
- `cutoff_radius::Real = 8.0`: Longest allowable distance to a neighbor, in Angstroms
- `max_num_nbr::Integer = 12`: Largest allowable number of neighbors
- `dist_decay_func = inverse_square`: Function by which to assign edge weights according to distance between neighbors

# Note
`max_num_nbr` is a "soft" limit – if multiple neighbors are at the same distance, the full neighbor list may be longer.
"""
function AtomGraph(
    crys::Crystal;
    id::String = "",
    cutoff_radius::Real = 8.0,
    max_num_nbr::Integer = 12,
    dist_decay_func::Function = inverse_square,
)
    adj_mat, elements = build_graph(
        crys;
        cutoff_radius = cutoff_radius,
        max_num_nbr = max_num_nbr,
        dist_decay_func = dist_decay_func,
    )
    ag = AtomGraph(adj_mat, elements, id)
end

# pretty printing, short version
function Base.show(io::IO, ag::AtomGraph)
    st = "AtomGraph $(ag.id) with $(nv(ag.graph)) nodes, $(ne(ag.graph)) edges"
    print(io, st)
end

# pretty printing, long version
function Base.show(io::IO, ::MIME"text/plain", ag::AtomGraph)
    st = "AtomGraph $(ag.id) with $(nv(ag.graph)) nodes, $(ne(ag.graph)) edges\n\tatoms: $(ag.elements)"
    print(io, st)
end

"""
    normalized_laplacian(graph)

Compute the normalized graph Laplacian matrix of the input graph, defined as

``I - D^{-1/2} A D^{-1/2}``

where ``A`` is the adjacency matrix and ``D`` is the degree matrix.
"""
function normalized_laplacian(g::G) where {G<:LightGraphs.AbstractGraph}
    a = adjacency_matrix(g)
    d = vec(sum(a, dims = 1))
    inv_sqrt_d = diagm(0 => d .^ (-0.5f0))
    laplacian = Float32.(I - inv_sqrt_d * a * inv_sqrt_d)
    !any(isnan, laplacian) || throw(
        ArgumentError(
            "NaN values in graph Laplacian! This is most likely due to atomic separations larger than the specified cutoff distance leading to block zeros in the adjacency matrix...try increasing the cutoff distance or inspecting your structure to ensure the file is correct.",
        ),
    )
    return laplacian
end

normalized_laplacian(ag::AtomGraph) = ag.laplacian

elements(ag::AtomGraph) = ag.elements

# now visualization stuff...

"Get a list of colors to use for graph visualization."
function graph_colors(atno_list, seed_color = colorant"cyan4")
    atom_types = unique(atno_list)
    atom_type_inds = Dict(atom_types[i] => i for i = 1:length(atom_types))
    color_inds = [atom_type_inds[i] for i in atno_list]
    colors = distinguishable_colors(length(atom_types), seed_color)
    return colors[color_inds]
end

"Helper function for sorting because edge ordering isn't preserved when converting to SimpleGraph."
function lt_edge(
    e1::SimpleWeightedGraphs.SimpleWeightedEdge{<:Integer,<:Real},
    e2::SimpleWeightedGraphs.SimpleWeightedEdge{<:Integer,<:Real},
)
    if e1.src < e2.src
        return true
    elseif e1.dst < e2.dst
        return true
    else
        return false
    end
end

"Compute edge widths (proportional to weights on graph) for graph visualization."
function graph_edgewidths(ag::AtomGraph)
    edgewidths = []
    edges_sorted = sort([e for e in edges(ag.graph)], lt = lt_edge)
    for e in edges_sorted
        append!(edgewidths, e.weight)
    end
    return edgewidths
end

"Visualize a given graph."
function visualize(ag::AtomGraph)
    # gplot doesn't work on weighted graphs
    sg = SimpleGraph(adjacency_matrix(ag.graph))
    plt = gplot(
        sg,
        nodefillc = graph_colors(ag.elements),
        nodelabel = ag.elements,
        edgelinewidth = graph_edgewidths(ag),
    )
    display(plt)
end


text:
"Create AtomicSolution for SDE."
function Solutions.AtomicSolution(equation::AbstractEquationSDE{DT,TT}) where {DT,TT}
    AtomicSolutionSDE(equation.t₀, equation.q₀[begin], zeros(DT,equation.m), zeros(DT,equation.m))
end

"Create AtomicSolution for SDE."
function Solutions.AtomicSolution(solution::SolutionSDE{AT,TT}) where {DT, TT, AT <: AbstractArray{DT}}
    AtomicSolutionSDE(get_initial_conditions(solution, 1)..., zeros(DT,solution.nm), zeros(DT,solution.nm))
end

"Create AtomicSolution for partitioned SDE."
function Solutions.AtomicSolution(equation::AbstractEquationPSDE{DT,TT}) where {DT,TT}
    AtomicSolutionPSDE(equation.t₀, equation.q₀[begin], equation.p₀[begin], zeros(DT,equation.m), zeros(DT,equation.m))
end

"Create AtomicSolution for partitioned SDE."
function Solutions.AtomicSolution(solution::SolutionPSDE{AT,TT}) where {DT, TT, AT <: AbstractArray{DT}}
    AtomicSolutionPSDE(get_initial_conditions(solution, 1)..., zeros(DT,solution.nm), zeros(DT,solution.nm))
end


text:
# --- PUBLIC METHODS ----------------------------------------------------------- #
function get_reaction_for_rn_number(rn_number_array::Array{String,1})::Some

    try

        # initialize -
        reaction_object_array = Array{KEGGReaction,1}()
        for rn_number in rn_number_array
            rxn_object = get_reaction_for_rn_number(rn_number) |> check
            if (isnothing(rxn_object) == false)
                push!(reaction_object_array, rxn_object)
            end
        end

        # return -
        return Some(reaction_object_array)
    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function get_reaction_for_rn_number(rn_number::String)::Some

    try

        # initialize -
        kegg_reaction = KEGGReaction()
        kegg_reaction.reaction_number = rn_number

        # formulae the URL string -
        url_string = "$(_KEGG_GET_URL)/$(rn_number)"

        # make the call -
        http_body = http_get_call_with_url(url_string) |> check
        if (isempty(http_body) == true)
            return Some(nothing)
        end

        # split around the \n -
        record_array = string.(split(http_body, "\n"))

        # grab: equation section -
        equation_section = extract_db_file_section(record_array, "EQUATION")
        if (isnothing(equation_section) == true)
            return Some(nothing)
        end
        rxn_string = string(split(equation_section, repeat(" ", 3))[2] |> lstrip)
        kegg_reaction.reaction_markup = rxn_string

        # reaction components -
        rxn_component_array = string.(split(rxn_string, " <=> "))
        kegg_reaction.reaction_forward = rxn_component_array[1]
        kegg_reaction.reaction_reverse = rxn_component_array[2]

        # grab: enzyme section -> ec number
        enzyme_section = extract_db_file_section(record_array, "ENZYME")
        if (isnothing(enzyme_section) == true)
            return Some(nothing)
        end
        tmp_ec_number = string(split(enzyme_section, repeat(" ", 4))[2] |> lstrip)
        kegg_reaction.ec_number = "ec:$(tmp_ec_number)"

        # grab: enzyme name -
        name = extract_db_file_section(record_array, "NAME")
        kegg_reaction.enzyme_name = string(split(name, repeat(" ", 5))[2] |> lstrip)

        # build the st dictionary -
        std = extract_stoichiometric_dictionary(rxn_string)
        kegg_reaction.stoichiometric_dictionary = std

        # return -
        return Some(kegg_reaction)

    catch error
        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function get_reactions_for_ec_number(ec_number_array::Array{String,1})::Some

    try

        # no ec numbers ... then no reactions!
        if (isempty(ec_number_array) == true)
            return Some(nothing)
        end

        # initialize -
        reaction_array = Array{KEGGReaction,1}()

        # process -
        for ec_number in ec_number_array

            # do we have an error?
            local_reaction_array = get_reactions_for_ec_number(ec_number) |> check

            # no - then, if this is not nothing, add -
            if (isnothing(local_reaction_array) == false)
                for reaction in local_reaction_array
                    push!(reaction_array, reaction)
                end
            end
        end

        # return -
        return Some(reaction_array)

    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function get_reactions_for_ec_number(ec_number::String)::Some

    try

        # Check: does this contain a ec: prefix? Yes, chop it off -
        local_ec_number = ec_number
        if (occursin("ec:", ec_number) == true)
            local_ec_number = ec_number[4:end]
        end

        # initalize -
        record_array = KEGGReaction[]

        # formulate url -
        linkage_url_string = "$(_KEGG_LINKAGE_REACTION_URL)/$(local_ec_number)"

        # call to get linkage array -
        ec_rn_linkage_array = http_get_call_with_url(linkage_url_string) |> check

        # if the server returned something, it should be in the form ec# rn#
        # we need the rn# to pull down the reaction string -
        if (length(ec_rn_linkage_array) > 1)

            # ec_rn_linkage can have more than one association -
            # split around the \n, and then process each item
            P = reverse(split(chomp(ec_rn_linkage_array), '\n'))
            while (!isempty(P))

                # pop -
                local_record = pop!(P)

                # The second fragment is the rn# -
                # split around the \t
                rn_number = string(split(local_record, '\t')[2])

                # ok, run a second query to pull the reaction string down -
                reaction_string = chomp((http_get_call_with_url("$(_KEGG_LIST_URL)/$(rn_number)")) |> check)

                # if we have something in reaction string, create a string record -
                if (length(reaction_string) != 0 &&
                    contains(reaction_string, ";") == true)

                    # create a new Reaction wrapper -
                    reaction_wrapper = KEGGReaction()

                    # populate w/easy stuff -
                    reaction_wrapper.ec_number = ec_number
                    reaction_wrapper.reaction_number = rn_number

                    # parse the body string -
                    tmp_enzyme_name = string(split(reaction_string, ";")[1])
                    reaction_wrapper.enzyme_name = string(split(tmp_enzyme_name, "\t")[2])

                    # split into forward and reverse strings -
                    tmp_full_reaction_string = lstrip(chomp(string(split(reaction_string, ";")[2])))
                    reaction_wrapper.reaction_markup = tmp_full_reaction_string

                    tmp_string_frag = split(tmp_full_reaction_string, " <=> ")
                    if (length(tmp_string_frag) == 2)

                        # forward -
                        reaction_wrapper.reaction_forward = lowercase(string(tmp_string_frag[1]))

                        # reverse -
                        reaction_wrapper.reaction_reverse = lowercase(string(tmp_string_frag[2]))
                    else
                        reaction_wrapper.reaction_forward = missing
                        reaction_wrapper.reaction_reverse = missing
                    end

                    # for now ...
                    tmp_reaction_string = "$(reaction_wrapper.reaction_forward) <=> $(reaction_wrapper.reaction_reverse)"
                    reaction_wrapper.stoichiometric_dictionary = extract_stoichiometric_dictionary(tmp_reaction_string)

                    # cache -
                    push!(record_array, reaction_wrapper)
                end
            end
        end

        # check - zero entries?
        if (length(record_array) == 0)
            return Some(nothing)
        end

        # return list of records -
        return Some(record_array)
    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end
# ------------------------------------------------------------------------------ #


text:
using JudiLing
using Documenter
using DocumenterLaTeX

makedocs(;
    modules=[JudiLing],
    authors="Xuefeng Luo",
    repo="https://github.com/MegamindHenry/JudiLing.jl/blob/{commit}{path}#L{line}",
    sitename="JudiLing.jl",
    format=LaTeX(platform = "none"),
    pages=[
        "Home" => "index.md",
        "Manual" => Any[
            "Make Cue Matrix" => "man/make_cue_matrix.md",
            "Make Semantic Matrix" => "man/make_semantic_matrix.md",
            "Cholesky" => "man/cholesky.md",
            "Make Adjacency Matrix" => "man/make_adjacency_matrix.md",
            "Make Yt Matrix" => "man/make_yt_matrix.md",
            "Find Paths" => "man/find_path.md",
            "Evaluation" => "man/eval.md",
            "Output" => "man/output.md",
            "Test Combo" => "man/test_combo.md",
            "Utils" => "man/utils.md"
        ]
    ],
)


text:
using PortHamiltonian

eigval, eigvec = eigen(ph);
num_freq = frequencies(eigval)
exact_freq = pi/2*collect(1:2:2*N);
err = exact_freq - num_freq;

@test abs(err[2]) < 1e-13

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
# Julia 0.5
using HDF5, JLD, DataFrames, Gadfly, Colors

for (fname, param) in (("perf_density", :density),
                       ("perf_size", :size))
    !isdir(fname) && mkdir(fname)
    df = readtable(fname*".csv", header=false, names=[:time, param, :test, :module])
    bdf = by(df, [:test,:module,param]) do df
        t = df[:time]
        DataFrame(min=minimum(t), median=median(t), mean=mean(t), p95=quantile(t,.95), max=maximum(t))
    end
    for d in groupby(bdf, :test)
        p = plot(d, x=param, y=:median, ymin=:min, ymax=:p95, color=:module,
                    Geom.line, Geom.ribbon, Scale.x_log10, Scale.y_log10,
                    Guide.Title(string(d[:test][1])), Guide.YLabel("time (s)"),
                    Theme(background_color=colorant"white"))
        draw(PNG(joinpath(fname, string(d[:test][1], ".png")), 4inch, 3inch), p)
    end
    # Relative plots
    rdf = by(bdf, [:test,param]) do df
        t = df[:min]
        mods = df[:module]
        baset = t[mods .== "Base"]
        testt = t[mods .!= "Base"]
        DataFrame(min=testt./baset)
    end
    for d in groupby(rdf, :test)
        p = plot(d, x=param, y=:min, Geom.line, Scale.x_log10,
                    Guide.Title(string(d[:test][1])), Guide.YLabel("min relative to base"),
                    Theme(background_color=colorant"white"))
        draw(PNG(joinpath(fname, string(d[:test][1], "_rel.png")), 4inch, 3inch), p)
    end
end


text:
# 3D Toy problem to test the solver

using luteos

P = P1() # Polynomial order of solution

mesh   = Mesh3D( "cube", P, N = 30 )#17)
master = Master3D( P )

mat = Material( E = 1, ν = 0.0, dim = mesh.dim )

# source function
function funcS( p::Array{Float64,2} )
  return fill(1.0, size(p,1), 3)
end

function funcB( p::Array{Float64,2} )
  return fill(0.0, size(p,1), 3)
end
function funcBNH( p::Array{Float64,2} )
  return fill(1.0, size(p,1), 3)
end

# prob = Problem( funcS, [funcB, funcB, funcB, funcB, funcB, funcB], [1,1,1,1,1,1],
#                 name="Example 1", bcnorm=false )
prob = Elas( "Example 1", mat, funcS,  [1,1,1,1,1,1], false,[funcB, funcB, funcB, funcB, funcB, funcBNH] )

@time (uhath, uh, σh ) = hdgSolve( master, mesh, prob, method = :GMRES )

# write solution
writeTecplot( "elas3D.dat", prob, mesh, uh, σh )

@printf("No. elements: %i", size(σh,3))


text:
__precompile__()
module GaussianRandomFields

# dependencies
using SpecialFunctions

# using PyPlot   # *****

using FastGaussQuadrature

using Arpack

using FFTW

using LinearAlgebra

using Statistics

# import statements
import Base: show, ndims, -

# import PyPlot: plot, surf, contour, contourf, plot_trisurf, tricontourf

# export statements
export IsotropicCovarianceStructure, AnisotropicCovarianceStructure, CovarianceFunction, ndims, apply # from covariance_functions.jl

export SeparableCovarianceFunction # from separable.jl

export Matern # from matern.jl

export Exponential # from exponential.jl

export SquaredExponential, Gaussian # from squaredexponential.jl

export AnisotropicExponential # from anisotropic_exponential.jl

export GaussianRandomField, sample, randdim # from gaussian_random_fields.jl

export Cholesky # from cholesky.jl

export Spectral # from spectral.jl

export KarhunenLoeve, rel_error # from karhunen_loeve.jl

export CirculantEmbedding # from circulant_embedding.jl

export GaussLegendre, EOLE, Simpson, Midpoint, Trapezoidal # from quadrature.jl

# export plot, surf, contour, contourf, tricontourf, plot_trisurf, plot_eigenvalues, plot_eigenfunction, plot_covariance_matrix # from plots.jl

export star, Lshape # from data/

# include statements
include("covariance_functions/covariance_functions.jl")

include("covariance_functions/matern.jl")

include("covariance_functions/exponential.jl")

include("covariance_functions/squared_exponential.jl")

include("covariance_functions/anisotropic_exponential.jl")

include("gaussian_random_fields.jl")

include("fem.jl")

include("generators/cholesky.jl")

include("generators/spectral.jl")

include("generators/quadrature.jl")

include("generators/karhunen_loeve.jl")

include("generators/circulant_embedding.jl")

include("analytic.jl")

include("separable.jl")

include("../data/star.jl")

include("../data/Lshape.jl")

# include("plot.jl")

end # module


text:
ConductanceStamp(G,pa::Int,pb::Int,N::Int) = (sparse([pa,pa,pb,pb],[pa,pb,pa,pb],[G,-G,-G,G],N,N),spzeros(N))
TransConductanceStamp(G,oa::Int,ob::Int,ia::Int,ib::Int,N::Int) = (sparse([oa,oa,ob,ob],[ia,ib,ia,ib],[G,-G,-G,G],N,N),spzeros(N))
VoltageSourceStamp(V,p::Int,m::Int,C::Int,N::Int) = (sparse([p,m,C,C],[C,C,p,m],[1,-1,-1,1],N,N),sparse([C],[1],[-V],N,1))
CurrentSourceStamp(I,p::Int,m::Int,N::Int) = (spzeros(N,N),sparse([p,m],[1,1],[-I,I],N,1))

function SymStamp(T::CircuitComponent{W},N::Int) where W <: TransistorModel
    name = T.name
    gmSym = Symbol("gm_{$name}")
    gdsSym = (Symbol("gds_{$name}"))
    IdSym = (Symbol("Id_{$name}"))
    syms = @variables  $gmSym $gdsSym $IdSym
    d = T.connections["d"]+1
    g = T.connections["g"]+1
    s = T.connections["s"]+1
    M1,I1 = TransConductanceStamp(syms[1],d,s,g,s,N)
    M2,I2 = ConductanceStamp(syms[2],d,s,N)
    #M3,I3 = ConductanceStamp(1/syms[1],g,s,N)
    M3,I3 = CurrentSourceStamp(syms[3],d,s,N)
    M4,I4 = CurrentSourceStamp(syms[1],d,s,N)
    addSymParameter.(Ref(T),syms)
    return M1+M2+M3,I1+I2+I3,syms
end

function MT_SymStamp(T::CircuitComponent{W},N::Int) where W <: TransistorModel
    name = T.name
    #gmSym = Symbol("gm_{$name}")
    #gdsSym = (Symbol("gds_{$name}"))
    IdSym = (Symbol("Id_{$name}"))
    syms = @variables  $gmSym $gdsSym $IdSym
    d = T.connections["d"]+1
    g = T.connections["g"]+1
    s = T.connections["s"]+1
    #M1,I1 = TransConductanceStamp(syms[1],d,s,g,s,N)
    #M2,I2 = ConductanceStamp(syms[2],d,s,N)
    #M3,I3 = ConductanceStamp(1/syms[1],g,s,N)
    M3,I3 = CurrentSourceStamp(syms[3],d,s,N)
    M4,I4 = CurrentSourceStamp(syms[1],d,s,N)
    addSymParameter.(Ref(T),syms)
    return M3,I3,syms
end

function SymStamp(T::CircuitComponent{W},N::Int,p::Num) where W <: TransistorModel
    name = T.name
    gmSym = Symbol("gm_{$name}")
    gdsSym = Symbol("gds_{$name}")
    CgsSym = Symbol("Cgs_{$name}")
    CgdSym = Symbol("Cgd_{$name}")
    #CdsSym = Symbol("Cds_{$name}")
    syms = @variables  $gmSym $gdsSym $CgsSym $CgdSym
    d = T.connections["d"]+1
    g = T.connections["g"]+1
    s = T.connections["s"]+1
    M1,I1 = TransConductanceStamp(syms[1],d,s,g,s,N)
    M2,I2 = ConductanceStamp(syms[2],d,s,N)
    M3,I3 = ConductanceStamp(1/syms[1],g,s,N)
    M4,I4 = ConductanceStamp(p*syms[3],g,s,N)
    M5,I5 = ConductanceStamp(p*syms[4],g,d,N)
    addSymParameter.(Ref(T),syms)
    #M6,I6 = ConductanceStamp(p*syms[5],d,s,N)
    return M1+M2+M3+M4+M5,I1+I2+I3+I4+I5,syms
end

function SymStamp(T::CircuitComponent{Resistor},N::Int,s=0)
    name = T.name
    G = Symbol("$name")
    syms = @variables $G
    p = T.connections["p"]+1
    n = T.connections["n"]+1
    M2,I2 = ConductanceStamp(1/syms[1],p,n,N)
    addSymParameter.(Ref(T),syms)
    return M2,I2,syms
end

function SymStamp(T::CircuitComponent{VoltageSource},N::Int,s=0)
    name = T.name
    G = Symbol("$name")
    syms = @variables $G
    p = T.connections["p"]+1
    n = T.connections["n"]+1
    j = T.connections["j"]+1
    addSymParameter.(Ref(T),syms)
    M2,I2 = VoltageSourceStamp(syms[1],p,n,j,N)
    return M2,I2,syms
end

function SymStamp(T::CircuitComponent{CurrentSource},N::Int,s=0)
    name = T.name
    G = Symbol("$name")
    syms = @variables $G
    p = T.connections["p"]+1
    n = T.connections["n"]+1
    addSymParameter.(Ref(T),syms)
    M2,I2 = CurrentSourceStamp(syms[1],p,n,N)
    return M2,I2,syms
end

function SymStamp(T::CircuitComponent{TransConductance},N::Int,s=0)
    name = T.name
    G = Symbol("$name")
    syms = @variables $G
    ip = T.connections["ip"]+1
    in = T.connections["im"]+1
    op = T.connections["op"]+1
    on = T.connections["om"]+1
    addSymParameter.(Ref(T),syms)
    M2,I2 = TransConductanceStamp(syms[1],op,on,ip,in,N)
    return M2,I2,syms
end


text:
using Plots
using DelimitedFiles
using LinearAlgebra
using ProgressMeter
using LsqFit
using KernelDensity
using Statistics
using LaTeXStrings
using StatsBase
using Measures
using JLD2
using PyCall
stats = pyimport("scipy.stats")
include("../Analysis/AlphaShapes/AlphaShapes.jl")
include("../Analysis/analysis-utils.jl")
# this is usually an "overnighter"
using ArgParse

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "-T"
            arg_type = Float64
			default = 10.			# seconds of simulated data
    end
    return parse_args(s)
end

args = parse_commandline()
@info args

model_seconds = args["T"]

@info "Simulating $model_seconds seconds of model data each iteration"

"""
A wrapper to run the model and collect the trajectories
"""
function Sim(N,alpha,tau,v,mu,mur=1.,dt=1.0/900.0)
    cmd = `../CUDAABP --initial-packing-fraction 0.5 -N $N -T $model_seconds -a 1 -b 1 -alpha $alpha -tau $tau -dt $dt -mu $mu -mur $mur -v $v -Dr 2.34 --random-seed $(Int(floor(time()))) -silent 1`
    run(cmd)
    T = readdlm("trajectories.txt",',')
    #h = Int(floor(size(T,1)/2))
    d_data = T[:,1:2]
    return d_data
end

d = Sim(10,1,1,1,30,1.);

 len = Int(ceil(size(d,1) / 10))

runs = 3
Ps = [200]

n = 15

alpha = collect(linspace(-2.0,0.0,n))
tau = collect(linspace(0.0,30.0,n))

tau = tau[2:end]

DR = 2.34
v0 = 13.19*2.0
mu = 31.624014716770823

bins_density = Float64.(logspace(1e-4,6.0,100))
bins_time = 100

Data = Dict()
prog = Progress(runs*n*n*length(Ps))
for a in alpha
    for b in tau
        for p in 1:length(Ps)
            #@info Ps[p]
            pdfs_model = zeros(runs,length(bins_density))
            Ts = zeros(runs,len,Ps[p],2)
            for r in 1:runs
                D = Sim(Ps[p],a,b,v0,mu,1.);
                for i in 1:size(Ts,2)
                    for j in 1:size(Ts,3)
                        # extract the trajectories
                        Ts[r,i,j,:] = D[(i-1)*Ps[p]+j,:]
                    end
                end
                if b == 0.0
                    # density is not recorded in simulation since the pre-factor is 0.0
                    d = zeros(size(Ts,2),size(Ts,3))
                    for i in 1:size(Ts,2)
                        d[i,:] = Density(Ts[r,i,:,1:2]./2.0)
                    end
                    d = flatten(d)
                else
                    # density data from CUDA code
                    d = readdlm("trajectories.txt",',')[:,4]
                end
                # extract second half of data
                h = Int(floor(size(d,1)/2.))
                kde = stats.gaussian_kde(d[h:end])
                pdfs_model[r,:] = kde(bins_density)
                next!(prog)
            end
            # save to dict
            Data["$(Ps[p])-$a-$b"] = pdfs_model
	    Data["Trajectory-$(Ps[p])-$a-$b"] = Ts
        end
    end
end
# write to this dir
JLD2.@save "alpha-tau-pd-15.jld" Data

max_stars_repo_path: Model Analysis/get-model-data-pd.jl

max_stars_repo_name: harveydevereux/CUDA-Whirligigs

max_issues_repo_path: Model Analysis/get-model-data-pd.jl

max_issues_repo_name: harveydevereux/CUDA-Whirligigs

max_forks_repo_path: Model Analysis/get-model-data-pd.jl

max_forks_repo_name: harveydevereux/CUDA-Whirligigs


text:
module CovidData

import DataFrames, CSV, ZipFile, WorldBankData, HTTP, Dates, GZip
using LinearAlgebra: dot
#using Statistics

export covidjhudata, statedata

"""
    covidjhudata()

Downloads most recent JHU CSSE data on covid cases, deaths, and recoveries.

Returns a DataFrame
"""
function covidjhudata(;redownload=false)

  url = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"
  ids = Symbol.(["Province/State","Country/Region", "Lat","Long"])
  dfs=map(name->begin
          res = HTTP.get(replace(url, "confirmed" => name))
          df = CSV.read(res.body)
          df = DataFrames.stack(df, DataFrames.Not(ids), variable_name=:Date, value_name = Symbol(name))
          end, ["confirmed", "deaths", "recovered"])

  df = DataFrames.join(dfs..., on=[ids..., :Date], kind=:outer)
  df[!,:Date] .= Dates.Date.(String.(df[!,:Date]), "mm/dd/yy") .+ Dates.Year(2000)
  for col in names(df)
    df[!,col]=replace(df[!,col], NaN=>missing)
  end
  DataFrames.rename!(df, [Symbol("Province/State") => :Province,
                          Symbol("Country/Region") => :Country])

  # merge country populations from world bank
  ambig = []
  df.iso2c = Vector{Union{Missing, String}}(undef, DataFrames.nrow(df))
  df.iso2c .= missing
  for c in unique(df.Country)
    res = WorldBankData.search_wdi("countries","name",Regex(c))
    if DataFrames.nrow(res) != 1
      push!(ambig, (c, res))
    else
      df[df.Country.==c, :iso2c] .= res.iso2c[1]
    end
  end
  # hand checked corrections
  df[df.Country.=="China", :iso2c] .= "CN"
  df[(df.Country.=="China") .& (df.Province.=="Hong Kong"), :iso2c] .= "HK"
  df[df.Country.=="Congo (Brazzaville)", :iso2c] .= "CG"
  df[df.Country.=="Congo (Kinshasa)", :iso2c] .= "CD"
  df[df.Country.=="Czechia", :iso2c] .= "CZ"
  df[df.Country.=="Guinea", :iso2c] .= "GN"
  df[df.Country.=="Korea, South", :iso2c] .= "KR"
  df[df.Country.=="Kyrgyzstan", :iso2c] .= "KG"
  df[df.Country.=="Niger", :iso2c] .= "CZ"
  df[df.Country.=="Nigeria", :iso2c] .= "NG"
  df[df.Country.=="Saint Lucia", :iso2c] .= "LC"
  df[df.Country.=="Saint Vincent and the Grenadines", :iso2c] .= "VC"
  df[df.Country.=="Slovakia", :iso2c] .= "SK"
  df[df.Country.=="South Africa", :iso2c] .= "ZA"
  df[df.Country.=="Sudan", :iso2c] .= "SD"
  df[df.Country.=="Dominica", :iso2c] .= "DM"
  df[df.Country.=="US", :iso2c] .= "US"

  csvfile=normpath(joinpath(dirname(Base.find_package("CovidData")),"..","data","cpop.csv"))
  if (redownload || !isfile(csvfile))
    mkpath(normpath(joinpath(dirname(Base.find_package("CovidData")),"..","data")))
    tmp = DataFrames.DataFrame(iso2c=unique(skipmissing(df.iso2c)))
    tmp.cpop = Vector{Union{Missing,Float64}}(undef, DataFrames.nrow(tmp))
    tmp.cpop .= missing
    for c in tmp.iso2c
      println(c)
      tmp.cpop[c.==tmp.iso2c] .=
        try
          wb = WorldBankData.wdi("SP.POP.TOTL",c, 2010)
          sort!(wb,[:year])
          r = DataFrames.nrow(wb)
          while ismissing(wb.SP_POP_TOTL[r]) && r>1
            r -= 1
          end
          wb.SP_POP_TOTL[r]
        catch
          missing
        end
    end
    CSV.write(csvfile, tmp)
  end
  tmp = CSV.read(csvfile)
  df = join(df, tmp, on=[:iso2c], kind=:left)

  # Add province populations for Canada
  zipfile=normpath(joinpath(dirname(Base.find_package("CovidData")),"..","data","17100009-eng.zip"))
  if !isfile(zipfile)
    download("https://www150.statcan.gc.ca/n1/tbl/csv/17100009-eng.zip", zipfile)
  end
  csvfile=normpath(joinpath(dirname(Base.find_package("CovidData")),"..","data","17100009.csv"))
  if !isfile(csvfile)
    datadir = normpath(joinpath(dirname(Base.find_package("CovidData")),"..","data"))
    run(`unzip $zipfile -d $datadir`)
  end
  cnd=CSV.read(csvfile)
  lastdate=maximum(cnd.REF_DATE)
  cnd=cnd[cnd.REF_DATE.==lastdate,[:GEO,:VALUE]]
  DataFrames.rename!(cnd, [:GEO => :Province, :VALUE => :ppop])
  df = join(df, cnd, on=[:Province], kind=:left)

  # province populations for China
  chinapop=Dict("Guangdong" => 111690000, "Shandong" => 100060000, "Henan" => 95590000,
               "Sichuan" => 83020000, "Jiangsu" => 80290000, "Hebei" => 75200000 ,
               "Hunan" => 68600000 , "Anhui" => 62550000 ,  "Hubei" => 59020000 ,
               "Zhejiang" => 56570000, "Guangxi" => 48850000, "Yunnan" =>  48010000 ,
               "Jiangxi" => 46220000 , "Liaoning" => 43690000 , "Fujian" => 39110000,
               "Shaanxi" => 38350000, "Heilongjiang" => 37890000,"Shanxi" => 36820000 ,
               "Guizhou" => 35550000,  "Chongqing" => 30750000, "Jilin" => 27170000,
               "Gansu" => 26260000 , "Inner Mongolia" => 25290000, "Xinjiang" => 24450000,
               "Shanghai" => 24180000, "Taiwan " => 23562318, "Beijing" => 21710000 ,
               "Tianjin" => 15570000 ,  "Hainan" => 9170000 , "Hong Kong Hong Kong" => 7335384 ,
               "Ningxia" => 6820000 , "Qinghai" => 5980000 ,  "Tibet" => 3370000 ,
               "Macau" => 644900  )
  for k in keys(chinapop)
    df.ppop[.!ismissing.(df.Province) .& (df.Province.==k)] .= chinapop[k]
  end
  cp = maximum(df.cpop[df.Country.=="China"])
  hk = (df.Country.=="China") .& (df.Province.=="Hong Kong")
  hkp = unique(df.cpop[hk])
  df.cpop[hk] .= cp
  df.ppop[hk] .= hkp

  # province populations for Australia
  auspop = Dict("New South Wales" => 7317500, "Victoria" => 5640900,
                "Queensland" => 4599400, "Western Australia" =>2366900,
                "South Australia" => 1659800, "Tasmania" =>511000,
                "Australian Capital Territory" => 366900,
                "Northern Territory" => 231200)
  for k in keys(auspop)
    df.ppop[.!ismissing.(df.Province) .& (df.Province.==k)] .= auspop[k]
  end
  return(df)
end

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

"""
    statedata(;filename="covidstates.Rda", policies=:dates, fillmissingmobility=false)

Loads data on US state case numbers, recoveries, deaths, policies, and
movement. This dataset combines
- JHU data on cases, deaths, & recoveries
- NYT data on cases & deaths
- Google mobility reports
- Policy data from Raifman et al
- Testing and hospitalizations from Covid Tracking Project
- Operations of businesses and hourly workers from homebase

The downloading and combination of datasets is done in R.

Options:

- `policies in [:dates, :indicators]` The downloaded policy date
  contains the date each policy went into effect. Set to `:indicators`
  if you want time-varying indicator variables for whether each policy
  was in place.

- `fillmissingmobility = false` should missing values from Google
  mobility reports be filled forward and backward in time? If `true`,
  then prior to the start of data reporting, we set google mobility
  reports to 0, and after the most recent date report, we set them to
  the last non-missing value.

Returns a DataFrame
"""
function statedata(;filename="covidstates.csv", policies=:dates, fillmissingmobility=false)
  fullpath = normpath(joinpath(dirname(Base.find_package("CovidData")),"..","data",filename))
  mkpath(normpath(joinpath(dirname(Base.find_package("CovidData")),"..","data")))
  if !isfile(fullpath) || Dates.days(Dates.now() - Dates.unix2datetime(mtime(fullpath)))>=1
    @info "$filename is more than 1 day old, attempting to recreate"
    rfile = normpath(joinpath(dirname(Base.find_package("CovidData")),"..","R","cases_and_policies.R"))
    try
      rwd =  normpath(joinpath(dirname(Base.find_package("CovidData")),"..","R"))
      cmd = `R -e "setwd(\"$rwd\"); source(\"$rfile\")"`
      run(cmd)
    catch err
      @warn "R encountered the following error:\n"*
      err.msg*
      "\n It is likely that some R packages are missing. Try running\n"*
      "$fullpath\n in R to debug.\n\n"
      "Using old data that was last modified on $(Dates.unix2datetime(mtime(fullpath)))"
    end
  end
  dat = CSV.read(fullpath, missingstring="NA")

  if policies==:indicators
    @info "Creating time-varying policy indicators from policy dates."
    pvars = [Symbol("Stay.at.home..shelter.in.place")  Symbol("End.relax.stay.at.home.shelter.in.place");
             Symbol("State.of.emergency") missing;
             Symbol("Date.closed.K.12.schools") missing;
             Symbol("Closed.gyms") Symbol("Repened.gyms");
             Symbol("Closed.movie.theaters") Symbol("Reopened.movie.theaters");
             Symbol("Closed.day.cares") missing;
             Symbol("Date.banned.visitors.to.nursing.homes") missing;
             Symbol("Closed.non.essential.businesses") Symbol("Reopen.businesses");
             Symbol("Closed.restaurants.except.take.out") Symbol("Reopen.restaurants")]
    for (r, p) in enumerate(pvars[:,1])
      newp = falses(size(dat,1))
      for st in unique(dat.state)
        ss = dat.state .== st
        day = unique(skipmissing(dat[ss,p]))
        if (length(day) != 1)
          newp[ss] .= false
        else
          newp[ss] .= (dat[ss, :date] .>= day)
        end
        if !ismissing(pvars[r,2])
          day = unique(skipmissing(dat[ss,pvars[r,2]]))
          if (length(day) == 1)
            newp[ss][dat[ss,:date] .>= day] .= false
          end
        end
      end
      dat[!,p] = newp
    end

  elseif policies != :dates
    @error "invalid option for policies"
  end
  dat = sort(dat, [:state, :date])
  if fillmissingmobility

    @info "Filling missing google mobility variables with 0's before\n"*
    "reporting began, and their last known value for the past few days."

    mvars = [:retail_and_recreation_percent_change_from_baseline,
             :grocery_and_pharmacy_percent_change_from_baseline,
             :parks_percent_change_from_baseline ,
             :transit_stations_percent_change_from_baseline,
             :workplaces_percent_change_from_baseline,
             :residential_percent_change_from_baseline,
             :percentchangebusinesses]
    for v in mvars
      i = ismissing.(dat[!,v]) .& (dat.date .<= Dates.Date("2020-03-01"))
      dat[i, v] .= 0
      for gdf in DataFrames.groupby(dat, :state)
        l = findlast(.!ismissing.(gdf[!,v]))
        lastval = (l === nothing) ? missing : gdf[l,v]
        gdf[ismissing.(gdf[:,v]),v] .= lastval
      end
    end
  end

  # UNACAST social distancing measures
  filename = "sds-v3-full-state.csv.gz"
  fullpath = normpath(joinpath(dirname(Base.find_package("CovidData")),"..","data",filename))
  if isfile(fullpath)
    if (Dates.days(Dates.now() - Dates.unix2datetime(mtime(fullpath)))>=1)
      @info "Unacast data is more than 1 day old, attempting to redownload"
      rfile = normpath(joinpath(dirname(Base.find_package("CovidData")),"..","R","unacast.R"))
      try
        rwd =  normpath(joinpath(dirname(Base.find_package("CovidData")),"..","R"))
        cmd = `R -e "setwd(\"$rwd\"); source(\"$rfile\")"`
        run(cmd)
      catch err
        @warn "R encountered the following error:\n"*
        err.msg*
        "\n It is likely that some R packages are missing. Try running\n"*
        "$fullpath\n in R to debug.\n\n"
        "Using old unacast data that was last modified on $(Dates.unix2datetime(mtime(fullpath)))"
      end
    end
    unacast = GZip.open(fullpath) do io
      CSV.read(io)
    end
    dat = DataFrames.join(dat, unacast, on=[:fips=>:state_fips, :state=>:state_name,
                                            :ST=>:state_code, :date=>:date])
  else
    @warn "Unacast social distancing scoreboard data not found\n"*
    "If you want Unacast social distancing variables included in data\n."*
    "Download $filename and save as $fullpath \n"
  end

  return(dat)
end

end # module


text:
# SPDX-License-Identifier: X11
# 2020-12-26
# Resistors in Parallel (200pt)

using Printf: @printf;

readline()
A = parse.(Int, (split ∘ readline)())

@printf("%.6f\n", (inv ∘ sum ∘ map)(inv, A))

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
abstract type AbstractPerson end
abstract type AbstractStudent <: AbstractPerson end
abstract type AbstractTeacher <: AbstractPerson end

struct Person <: AbstractPerson
  name::String
end

struct Student <: AbstractStudent
  name::String
  grade::Int
  hobby::String
end

struct MusicStudent <: AbstractStudent
  grade::Int
end

struct Teacher <: AbstractTeacher
  name::String
  grade::Int
end

get_name(x::AbstractPerson) = x.name

stud = Student("name", 5, "swim")
mus = MusicStudent(7)

get_name(mus)

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency