text:
##
# Represent a banded matrix
# [ a_11 a_12
#   a_21 a_22 a_23
#   a_31 a_32 a_33 a_34
#        a_42 a_43 a_44  ]
# ordering the data like  (cobbndsmns first)
#       [ *      a_12   a_23    a_34
#         a_11   a_22   a_33    a_44
#         a_21   a_32   a_43    *
#         a_31   a_42   *       *       ]
###

function _BandedMatrix end

mutable struct BandedMatrix{T, CONTAINER, RAXIS} <: AbstractBandedMatrix{T}
    data::CONTAINER  # l+u+1 x n (# of columns)
    raxis::RAXIS # axis for rows (col axis comes from data)
    l::Int # lower bandwidth ≥0
    u::Int # upper bandwidth ≥0
    global function _BandedMatrix(data::AbstractMatrix{T}, raxis::AbstractUnitRange, l, u) where {T}
        if size(data,1) ≠ l+u+1  && !(size(data,1) == 0 && -l > u)
           error("Data matrix must have number rows equal to number of bands")
        else
            new{T, typeof(data), typeof(raxis)}(data, raxis, l, u)
        end
    end
end

_BandedMatrix(data::AbstractMatrix, m::Integer, l, u) = _BandedMatrix(data, Base.OneTo(m), l, u)

const DefaultBandedMatrix{T} = BandedMatrix{T,Matrix{T},OneTo{Int}}

bandedcolumns(_) = BandedColumns{UnknownLayout}()
bandedcolumns(::ML) where ML<:AbstractStridedLayout = BandedColumns{ML}()
bandedcolumns(::ML) where ML<:AbstractFillLayout = BandedColumns{ML}()

MemoryLayout(A::Type{BandedMatrix{T,Cont,Axes}}) where {T,Cont,Axes} = bandedcolumns(MemoryLayout(Cont))

## Constructors

"""
    BandedMatrix{T}(undef, (n, m), (l, u))

returns an uninitialized `n`×`m` banded matrix of type `T` with bandwidths `(l,u)`.
"""
BandedMatrix{T, C, OneTo{Int}}(::UndefInitializer, (n,m)::NTuple{2,Integer}, (a,b)::NTuple{2,Integer}) where {T<:BlasFloat, C<:AbstractMatrix{T}} =
    _BandedMatrix(C(undef,max(0,b+a+1),m), n, a, b)
BandedMatrix{T, C, OneTo{Int}}(::UndefInitializer, (n,m)::NTuple{2,Integer}, (a,b)::NTuple{2,Integer}) where {T, C<:AbstractMatrix{T}} =
    _BandedMatrix(C(undef,max(0,b+a+1),m),n,a,b)
BandedMatrix{T, C, OneTo{Int}}(::UndefInitializer, (n,m)::NTuple{2,Integer}, (a,b)::NTuple{2,Integer})  where {T<:Number, C<:AbstractMatrix{T}} =
    _BandedMatrix(fill!(similar(C, max(0,b+a+1),m), zero(T)),n,a,b)

BandedMatrix{T, C}(::UndefInitializer, (n,m)::NTuple{2,Integer}, (a,b)::NTuple{2,Integer}) where {T<:BlasFloat, C<:AbstractMatrix{T}} =
    _BandedMatrix(C(undef,max(0,b+a+1),m), n, a, b)
BandedMatrix{T, C}(::UndefInitializer, (n,m)::NTuple{2,Integer}, (a,b)::NTuple{2,Integer}) where {T, C<:AbstractMatrix{T}} =
    _BandedMatrix(C(undef,max(0,b+a+1),m),n,a,b)
BandedMatrix{T, C}(::UndefInitializer, (n,m)::NTuple{2,Integer}, (a,b)::NTuple{2,Integer})  where {T<:Number, C<:AbstractMatrix{T}} =
    _BandedMatrix(fill!(similar(C, max(0,b+a+1),m), zero(T)),n,a,b)

BandedMatrix{T}(::UndefInitializer, nm::NTuple{2,Integer}, ab::NTuple{2,Integer}) where T =
BandedMatrix{T, Matrix{T}}(undef,nm,ab)
BandedMatrix{T}(::UndefInitializer, nm::NTuple{2,OneTo{Int}}, ab::NTuple{2,Integer}) where T =
    BandedMatrix{T}(undef, length.(nm), ab)

@deprecate BandedMatrix{T, C}(::UndefInitializer, n::Integer, ::Colon, a::Integer, b::Integer)  where {T, C<:AbstractMatrix{T}} BandedMatrix{T,C}(undef,n,n+b,a,b)
@deprecate BandedMatrix{T, C}(::UndefInitializer, n::Integer, m::Integer, a::Integer, b::Integer) where {T, C<:AbstractMatrix{T}} BandedMatrix{T, C}(undef, (n,m), (a,b))
@deprecate BandedMatrix{T}(::UndefInitializer, n::Integer, m::Integer, a::Integer, b::Integer)  where {T} BandedMatrix{T}(undef,(n,m),(a,b))
@deprecate BandedMatrix{T}(::UndefInitializer, n::Integer, ::Colon, a::Integer, b::Integer)  where {T} BandedMatrix{T}(undef,(n,n+b),(a,b))

BandedMatrix{V}(M::BandedMatrix) where {V} = _BandedMatrix(AbstractMatrix{V}(M.data), M.raxis, M.l, M.u)
BandedMatrix(M::BandedMatrix{V}) where {V} = _BandedMatrix(AbstractMatrix{V}(M.data), M.raxis, M.l, M.u)

convert(::Type{BandedMatrix{V}}, M::BandedMatrix{V}) where {V} = M
convert(::Type{BandedMatrix{V}}, M::BandedMatrix) where {V} =
        _BandedMatrix(convert(AbstractMatrix{V}, M.data), M.raxis, M.l, M.u)
convert(::Type{BandedMatrix}, M::BandedMatrix) = M

convert(BM::Type{BandedMatrix{T,C}}, M::BandedMatrix{T,C}) where {T,C<:AbstractMatrix{T}} = M
convert(BM::Type{BandedMatrix{T,C,AXIS}}, M::BandedMatrix{T,C,AXIS}) where {T,C<:AbstractMatrix{T},AXIS} = M
convert(BM::Type{BandedMatrix{T,C,OneTo{Int}}}, M::BandedMatrix{T,C,OneTo{Int}}) where {T,C<:AbstractMatrix{T}} = M

function convert(::Type{BandedMatrix{<:, C}}, M::BandedMatrix) where C
    M.data isa C && return M
    _BandedMatrix(convert(C, M.data), M.raxis, M.l, M.u)
end
function convert(BM::Type{BandedMatrix{T, C}}, M::BandedMatrix) where {T, C <: AbstractMatrix{T}}
    M.data isa C && return M
    _BandedMatrix(convert(C, M.data), M.raxis, M.l, M.u)
end
convert(BM::Type{BandedMatrix{T, C, AXIS}}, M::BandedMatrix) where {T, C <: AbstractMatrix{T}, AXIS} =
    _BandedMatrix(convert(C, M.data), convert(AXIS, M.raxis), M.l, M.u)
convert(BM::Type{BandedMatrix{T, C, OneTo{Int}}}, M::BandedMatrix) where {T, C <: AbstractMatrix{T}} =
    _BandedMatrix(convert(C, M.data), convert(OneTo{Int}, M.raxis), M.l, M.u)

for MAT in (:AbstractBandedMatrix, :AbstractMatrix, :AbstractArray)
    @eval begin
        convert(::Type{$MAT{T}}, M::BandedMatrix) where {T} = convert(BandedMatrix{T}, M)
        convert(::Type{$MAT}, M::BandedMatrix) = M
        $MAT{T}(M::BandedMatrix) where {T} = BandedMatrix{T}(M)
        $MAT(M::BandedMatrix{T}) where {T} = BandedMatrix{T}(M)
    end
end

function convert(::Type{BandedMatrix{<:,C,OneTo{Int}}}, M::AbstractMatrix) where {C}
    Container = typeof(convert(C, similar(M, 0, 0)))
    T = eltype(Container)
    ret = BandedMatrix{T, Container}(undef, size(M), bandwidths(M))
    for k=1:size(M,1),j=1:size(M,2)
        ret[k,j] = convert(T, M[k,j])
    end
    ret
end

function convert(::Type{BandedMatrix{T,C,OneTo{Int}}}, M::AbstractMatrix) where {T, C}
    Container = typeof(convert(C, similar(M, T, 0, 0)))
    ret = BandedMatrix{T, Container}(undef, size(M), bandwidths(M))
    for k=1:size(M,1),j=1:size(M,2)
        ret[k,j] = convert(T, M[k,j])
    end
    ret
end

function convert(::Type{BandedMatrix{T,C,OneTo{Int}}}, M::AbstractMatrix) where {T, C<:AbstractMatrix{T}}
    Container = typeof(convert(C, similar(M, T, 0, 0)))
    ret = BandedMatrix{T, Container}(undef, size(M), bandwidths(M))
    for k=1:size(M,1),j=1:size(M,2)
        ret[k,j] = convert(T, M[k,j])
    end
    ret
end

convert(::Type{BandedMatrix{<:, C}}, M::AbstractMatrix) where {C} = convert(BandedMatrix{<:,C,OneTo{Int}}, M)
convert(::Type{BandedMatrix{T,C}}, M::AbstractMatrix) where {T, C} = convert(BandedMatrix{T,C,OneTo{Int}}, M)

convert(::Type{BandedMatrix{T}}, M::AbstractMatrix) where {T} =
    convert(BandedMatrix{T, typeof(similar(M, T, 0, 0))}, M)

convert(::Type{BandedMatrix}, M::AbstractMatrix) = convert(BandedMatrix{eltype(M)}, M)
convert(::Type{DefaultBandedMatrix}, M::AbstractMatrix{T}) where T = convert(DefaultBandedMatrix{T}, M)

copy(B::BandedMatrix) = _BandedMatrix(copy(B.data), B.raxis, B.l, B.u)

promote_rule(::Type{BandedMatrix{T1, C1}}, ::Type{BandedMatrix{T2, C2}}) where {T1,C1, T2,C2} =
    BandedMatrix{promote_type(T1,T2), promote_type(C1, C2)}

for (op,bop) in ((:(rand),:brand),)
    @eval begin
        $bop(::Type{T},n::Integer,m::Integer,a::Integer,b::Integer) where {T} =
            _BandedMatrix($op(T,max(0,b+a+1),m),n,a,b)
        $bop(::Type{T},n::Integer,a::Integer,b::Integer) where {T} = $bop(T,n,n,a,b)
        $bop(::Type{T},n::Integer,::Colon,a::Integer,b::Integer) where {T} = $bop(T,n,n+b,a,b)
        $bop(::Type{T},::Colon,m::Integer,a::Integer,b::Integer) where {T} = $bop(T,m+a,m,a,b)
        $bop(n::Integer,m::Integer,a::Integer,b::Integer) = $bop(Float64,n,m,a,b)
        $bop(n::Integer,a::Integer,b::Integer) = $bop(n,n,a,b)

        $bop(::Type{T},n::Integer,m::Integer,a) where {T} = $bop(T,n,m,-a[1],a[end])
        $bop(::Type{T},n::Number,::Colon,a) where {T} = $bop(T,n,:,-a[1],a[end])
        $bop(::Type{T},::Colon,m::Integer,a) where {T} = $bop(T,:,m,-a[1],a[end])
        $bop(::Type{T},n::Integer,a) where {T} = $bop(T,n,-a[1],a[end])
        $bop(n::Integer,m::Integer,a) = $bop(Float64,n,m,-a[1],a[end])
        $bop(n::Integer,a) = $bop(n,-a[1],a[end])

        $bop(B::AbstractMatrix) =
            $bop(eltype(B),size(B,1),size(B,2),bandwidth(B,1),bandwidth(B,2))
    end
end

"""
    brand(T,n,m,l,u)

Creates an `n×m` banded matrix  with random numbers in the bandwidth of type `T` with bandwidths `(l,u)`
"""
brand

## Conversions from AbstractArrays, we include FillArrays in case `zeros` is ever faster
BandedMatrix{T}(A::AbstractMatrix, bnds::NTuple{2,Integer}) where T =
    BandedMatrix{T, Matrix{T}}(A, bnds)
function BandedMatrix{T, C}(A::AbstractMatrix, bnds::NTuple{2,Integer}) where {T, C <: AbstractMatrix{T}}
    (n,m) = size(A)
    (l,u) = bnds
    ret = BandedMatrix{T, C}(undef, (n, m), (l, u))
    @inbounds for j = 1:m, k = max(1,j-u):min(n,j+l)
        inbands_setindex!(ret, A[k,j], k, j)
    end
    ret
end

BandedMatrix(A::AbstractMatrix{T}, bnds::NTuple{2,Integer}) where T =
    BandedMatrix{T}(A, bnds)

BandedMatrix{V}(Z::Zeros{T,2}, bnds::NTuple{2,Integer}) where {T,V} =
    _BandedMatrix(zeros(V,max(0,sum(bnds)+1),size(Z,2)),size(Z,1),bnds...)

BandedMatrix(E::Eye{T}, bnds::NTuple{2,Integer}) where T = BandedMatrix{T}(E, bnds)
function BandedMatrix{T}(E::Eye, bnds::NTuple{2,Integer}) where T
    ret=BandedMatrix(Zeros{T}(E), bnds)
    ret[band(0)] .= one(T)
    ret
end

BandedMatrix{T,C}(A::AbstractMatrix) where {T, C<:AbstractMatrix{T}} = BandedMatrix{T,C}(A, bandwidths(A))
BandedMatrix{T}(A::AbstractMatrix) where T = BandedMatrix{T}(A, bandwidths(A))

# use bandeddata if possible
_BandedMatrix(::BandedColumns, A::AbstractMatrix) = _BandedMatrix(copy(bandeddata(A)), size(A,1), bandwidths(A)...)
_BandedMatrix(_, A::AbstractMatrix) = BandedMatrix(A, bandwidths(A))
BandedMatrix(A::AbstractMatrix) = _BandedMatrix(MemoryLayout(typeof(A)), A)

"""
    BandedMatrix{T}(kv::Pair, (m,n), (l,u))

Construct a m × n BandedMatrix with bandwidths (l,u)
from `Pair`s of diagonals and vectors.
Vector `kv.second` will be placed on the `kv.first` diagonal.
"""
function BandedMatrix{T}(kv::Tuple{Vararg{Pair{<:Integer,<:AbstractVector}}},
                         mn::NTuple{2,Integer},
                         lu::NTuple{2,Integer}) where T
    m,n = mn
    l,u = lu
    data = zeros(T, u+l+1, n)
    for (k,v) in kv
        p = length(v)
        if k ≤ 0
            data[u-k+1,1:p] = v
        else
            data[u-k+1,k+1:k+p] = v
        end
    end

    return _BandedMatrix(data, m, l, u)
end

BandedMatrix(kv::Tuple{Vararg{Pair{<:Integer,<:AbstractVector}}},
                         mn::NTuple{2,Integer},
                         lu::NTuple{2,Integer}) =
    BandedMatrix{promote_type(map(x -> eltype(x.second), kv)...)}(kv, mn, lu)

function BandedMatrix{T}(kv::Tuple{Vararg{Pair{<:Integer,<:AbstractVector}}},
                         nm::NTuple{2,Integer}) where T
    u = mapreduce(x -> x.first, max, kv)
    l = -mapreduce(x -> x.first, min, kv)
    BandedMatrix{T}(kv, nm, (l,u))
end

BandedMatrix(kv::Tuple{Vararg{Pair{<:Integer,<:AbstractVector}}}, nm::NTuple{2,Integer}) =
    BandedMatrix{promote_type(map(x -> eltype(x.second), kv)...)}(kv, nm)

function BandedMatrix{T}(kv::Pair{<:Integer,<:AbstractVector}...) where T
    n = mapreduce(x -> length(x.second) + abs(x.first), max, kv)
    u = mapreduce(x -> x.first, max, kv)
    l = -mapreduce(x -> x.first, min, kv)
    BandedMatrix{T}(kv, (n,n), (l,u))
end

"""
    BandedMatrix(kv::Pair{<:Integer,<:AbstractVector}...)

Construct a square matrix from `Pair`s of diagonals and vectors.
Vector `kv.second` will be placed on the `kv.first` diagonal.
"""
BandedMatrix(kv::Pair{<:Integer,<:AbstractVector}...) =
    BandedMatrix{promote_type(map(x -> eltype(x.second), kv)...)}(kv...)

"""
    similar(bm::BandedMatrix, [T::Type], [n::Integer, m::Integer, [l::Integer, u::Integer]])
    similar(bm::BandedSubBandedMatrix, [T::Type], [n::Integer, m::Integer, [l::Integer, u::Integer]])

Creates a banded matrix similar in type to the input. Where the eltype `T`, the
sizes `n`, and `m`, the band limits `l`, and `u` are not provided, they default
to the values in the input banded matrix.
"""
function similar(bm::BandedMatrix, T::Type=eltype(bm),
                      n::Integer=size(bm, 1), m::Integer=size(bm, 2),
                      l::Integer=bandwidth(bm, 1), u::Integer=bandwidth(bm, 2))
    data = similar(bm.data, T, max(0, u+l+1), m)
    _BandedMatrix(data, n, l, u)
end

similar(bm::AbstractBandedMatrix, n::Integer, m::Integer) = similar(bm, eltype(bm), n, m)
similar(bm::AbstractBandedMatrix, n::Integer, m::Integer, l::Integer, u::Integer) =
    similar(bm, eltype(bm), m, n, l, u)
similar(bm::AbstractBandedMatrix, nm::Tuple{<:Integer,<:Integer}) = similar(bm, nm...)

## Abstract Array Interface

axes(A::BandedMatrix) = (A.raxis, axes(A.data,2))
size(A::BandedMatrix) = (length(A.raxis), size(A.data,2))

## banded matrix interface
bandeddata(A::BandedMatrix) = A.data
bandwidths(A::BandedMatrix) = (A.l , A.u)

IndexStyle(::Type{BandedMatrix{T}}) where {T} = IndexCartesian()

# ~~ getindex ~~

# fast method used below
@inline inbands_getindex(data::AbstractMatrix, u::Integer, k::Integer, j::Integer) =
    data[u + k - j + 1, j]

@inline function inbands_getindex(A::BandedMatrix, k::Integer, j::Integer)
    # it takes a bit of time to extract A.data, A.u since julia checks if those fields exist
    # the @inbounds here will suppress those checks
    @inbounds r = inbands_getindex(A.data, A.u, k, j)
    r
end

# work around for Any matrices
_zero(T) = zero(T)
_zero(::Type{Any}) = nothing

# banded get index, used for banded matrices with other data types
@inline function banded_getindex(data::AbstractMatrix, l::Integer, u::Integer, k::Integer, j::Integer)
    if -l ≤ j-k ≤ u
        inbands_getindex(data, u, k, j)
    else
        _zero(eltype(data))
    end
end

# scalar - integer - integer
@inline function getindex(A::BandedMatrix, k::Integer, j::Integer)
    @boundscheck checkbounds(A, k, j)
    @inbounds r = banded_getindex(A.data, A.l, A.u, k, j)
    r
end

# ~ indexing along a band
# we reduce it to converting a View

# type to represent a view of a band
const BandedMatrixBand{T} = SubArray{T, 1, ReshapedArray{T,1,BandedMatrix{T},
                                Tuple{MultiplicativeInverses.SignedMultiplicativeInverse{Int}}}, Tuple{BandSlice}, false}

band(V::BandedMatrixBand) = first(parentindices(V)).band.i

# gives a view of the parent's data matrix
function dataview(V::BandedMatrixBand)
    A = parent(parent(V))
    b = band(V)
    m,n = size(A)
    if b > 0
        view(A.data, A.u - b + 1, b+1:min(n,m+b))
    elseif b == 0
        view(A.data, A.u - b + 1, 1:min(n,m))
    else # b < 0
        view(A.data, A.u - b + 1, 1:min(n,m+b))
    end
end

function convert(::Type{Vector{T}}, V::BandedMatrixBand) where T
    A = parent(parent(V))
    if -A.l ≤ band(V) ≤ A.u
        Vector{T}(dataview(V))
    else
        zeros(T, length(V))
    end
end

convert(::Type{Array{T}}, A::BandedMatrixBand) where T = convert(Vector{T}, A)
convert(::Type{Array}, A::BandedMatrixBand) = convert(Vector{eltype(A)}, A)
convert(::Type{Vector}, A::BandedMatrixBand)= convert(Vector{eltype(A)}, A)

convert(::Type{AbstractArray{T}}, A::BandedMatrixBand{T}) where T = A
convert(::Type{AbstractVector{T}}, A::BandedMatrixBand{T}) where T = A
convert(::Type{AbstractArray}, A::BandedMatrixBand{T}) where T = A
convert(::Type{AbstractVector}, A::BandedMatrixBand{T}) where T = A

convert(::Type{AbstractArray{T}}, A::BandedMatrixBand) where T = convert(Vector{T}, A)
convert(::Type{AbstractVector{T}}, A::BandedMatrixBand) where T = convert(Vector{T}, A)

# ~ indexing along a row

# give range of data matrix corresponding to colrange/rowrange
data_colrange(A::BandedMatrix{T}, i::Integer) where {T} =
    (max(1,A.u+2-i):min(size(A,1)+A.u-i+1,size(A.data,1))) .+
                                ((i-1)*size(A.data,1))

data_rowrange(A::BandedMatrix{T}, i::Integer) where {T} = range((i ≤ 1+A.l ? A.u+i : (i-A.l)*size(A.data,1)) ,
                                size(A.data,1)-1 ,  # step size
                                i+A.u ≤ size(A,2) ? A.l+A.u+1 : size(A,2)-i+A.l+1)

# ~~ setindex! ~~

# ~ Special setindex methods ~

# fast method used below
@inline @propagate_inbounds function inbands_setindex!(data::AbstractMatrix{T}, u::Integer, v, k::Integer, j::Integer) where {T}
    data[u + k - j + 1, j] = convert(T, v)::T
    v
end

# slow fall back method
@inline @propagate_inbounds function inbands_setindex!(A::BandedMatrix, v, k::Integer, j::Integer)
    # it takes a bit of time to extract A.data, A.u since julia checks if those fields exist
    # the @inbounds here will suppress those checks
    @inbounds r = inbands_setindex!(A.data, A.u, v, k, j)
    r
end

@inline @propagate_inbounds function banded_setindex!(data::AbstractMatrix, l::Integer, u::Integer, v, k::Integer, j::Integer)
    if -l ≤ j-k ≤ u
        inbands_setindex!(data, u, v, k, j)
    elseif !iszero(v)  # allow setting outside bands to zero
        throw(BandError(data,j-k))
    else # v == 0
        v
    end
end

# scalar - integer - integer
@inline @propagate_inbounds function setindex!(A::BandedMatrix, v, k::Integer, j::Integer)
    @boundscheck checkbounds(A, k, j)
    @inbounds r = banded_setindex!(A.data, A.l, A.u, v, k ,j)
    r
end

# matrix - colon - colon
@inline @propagate_inbounds function setindex!(A::BandedMatrix{T}, v::AbstractMatrix, kr::Colon, jr::Colon) where {T}
    @boundscheck checkdimensions(size(A), size(v))
    @boundscheck checkbandmatch(A, v, kr, jr)

    for j=1:size(A,2), k=colrange(A,j)
        @inbounds A[k,j] = v[k,j]
    end
    A
end

@propagate_inbounds function setindex!(A::BandedMatrix{T}, v::AbstractVector, ::Colon) where {T}
    A[:, :] = reshape(v,size(A))
end

# ~ indexing along a band

# vector - band - colon
@inline @propagate_inbounds function setindex!(A::BandedMatrix{T}, V::AbstractVector, b::Band) where {T}
    @boundscheck checkband(A, b)
    @boundscheck checkdimensions(diaglength(A, b), V)
    row = A.u - b.i + 1
    data, i = A.data, max(b.i + 1, 1)
    for v in V
        data[row, i] = convert(T, v)::T
        i += 1
    end
    V
end

# ~ indexing along columns
# vector - colon - integer -- A[:, 1] = [1, 2, 3] - not allowed
@inline @propagate_inbounds function setindex!(A::BandedMatrix{T}, V::AbstractVector, kr::Colon, j::Integer) where {T}
    @boundscheck checkbounds(A, kr, j)
    @boundscheck checkdimensions(1:size(A,1), V)
    @boundscheck checkbandmatch(A,V,:,j)

    A.data[data_colrange(A,j)] = V[colrange(A,j)]
    V
end

# vector - BandRange - integer -- A[1, BandRange] = 2
@propagate_inbounds setindex!(A::BandedMatrix, V::AbstractVector, ::BandRangeType, j::Integer) =
    (A[colrange(A, j), j] = V) # call range method

# vector - range - integer -- A[1:3, 1] = [1, 2, 3]
@inline @propagate_inbounds function setindex!(A::BandedMatrix, V::AbstractVector, kr::AbstractRange, j::Integer)
    @boundscheck checkbounds(A, kr, j)
    @boundscheck checkdimensions(kr, V)
    @boundscheck checkbandmatch(A, V, kr, j)

    a = colstart(A, j)
    b = colstop(A, j)

    data, u, i = A.data, A.u, 0
    for v in V
        k = kr[i+=1]
        if a ≤ k ≤ b
            inbands_setindex!(data, u, v, k, j)
        end
    end
    V
end

# ~ indexing along a row

# vector - integer - colon -- A[1, :] = [1, 2, 3] - not allowed
@inline @propagate_inbounds function setindex!(A::BandedMatrix{T}, V::AbstractVector, k::Integer, jr::Colon) where {T}
    @boundscheck if k < 1 || k > size(A,1)
        throw(BoundsError(A, (k, jr)))
    end
    @boundscheck if size(A,2) ≠ length(V)
        throw(DimensionMismatch("tried to assign $(length(V)) vector to $(size(A,1)) destination"))
    end

    for j = 1:rowstart(A,k)-1
        if V[j] ≠ zero(T)
            throw(BandError(A, _firstdiagrow(A, k)))
        end
    end
    for j = rowstop(A,j)+1:size(A,2)
        if V[j] ≠ zero(T)
            throw(BandError(A, _firstdiagrow(A, k)))
        end
    end

    A.data[data_rowrange(A,k)] = V[rowrange(A,k)]
    V
end

# vector - integer - BandRange -- A[1, BandRange] = [1, 2, 3]
@propagate_inbounds setindex!(A::BandedMatrix, V::AbstractVector, k::Integer, ::BandRangeType) =
    (A[k, rowstart(A, k):rowstop(A, k)] = V) # call range method

# vector - integer - range -- A[1, 2:3] = [3, 4]
@inline @propagate_inbounds function setindex!(A::BandedMatrix, V::AbstractVector, k::Integer, jr::AbstractRange)
    @boundscheck checkbounds(A, k, jr)
    @boundscheck checkdimensions(jr, V)
    @boundscheck checkbandmatch(A, V, k, jr)

    a = rowstart(A, k)
    b = rowstop(A, k)

    data, u, i = A.data, A.u, 0
    for v in V
        j = jr[i+=1]
        if a ≤ j ≤ b
            inbands_setindex!(data, u, v, k, j)
        end
    end
    V
end

# ~ indexing over a rectangular block

# matrix - range - range
@inline @propagate_inbounds function setindex!(A::BandedMatrix, V::AbstractMatrix, kr::AbstractRange, jr::AbstractRange)
    @boundscheck checkbounds(A, kr, jr)
    @boundscheck checkdimensions(kr, jr, V)
    @boundscheck checkbandmatch(A, V, kr, jr)

    data, u, l = A.data, A.u, A.l
    jj = 1
    for j in jr
        kk = 1
        for k in kr
            if -l ≤ j - k ≤ u
                # we index V manually in column-major order
                inbands_setindex!(data, u, V[kk, jj], k, j)
                kk += 1
            end
        end
        jj += 1
    end
    V
end

# ~~ end setindex! ~~

function convert(::Type{Matrix}, A::BandedMatrix)
    ret=zeros(eltype(A),size(A,1),size(A,2))
    for j = 1:size(ret,2), k = colrange(ret,j)
        @inbounds ret[k,j] = A[k,j]
    end
    ret
end

function _banded_colval(B::AbstractMatrix)
    data = bandeddata(B)
    j = Vector{Int}(undef, length(data))
    n,m = size(data)
    for κ=1:n, ℓ=1:m
        j[κ+n*(ℓ-1)] = ℓ
    end
    j
end

function _banded_rowval(B::AbstractMatrix)
    data = bandeddata(B)
    i = Vector{Int}(undef, length(data))
    n,m = size(data)
    Bn = size(B,1)
    l,u = bandwidths(B)
    for κ=1:n, ℓ=1:m
        ii=κ+ℓ-u-1
        i[κ+n*(ℓ-1)] = min(max(ii,1),Bn)
    end
    i
end

function _banded_nzval(B::AbstractMatrix)
    data = bandeddata(B)
    i = Vector{Int}(undef, length(data))
    n,m = size(data)
    Bn = size(B,1)
    vb = copy(vec(data))
    l,u = bandwidths(B)
    for κ=1:n, ℓ=1:m
        ii=κ+ℓ-u-1
        if ii <1 || ii > Bn
            vb[κ+n*(ℓ-1)] = 0
        end
    end
    vb
end

sparse(B::BandedMatrix) = sparse(_banded_rowval(B), _banded_colval(B), _banded_nzval(B), size(B)...)

function _bidiagonalize!(A::AbstractMatrix{T}, M::BandedColumnMajor) where T
    m, n = size(A)
    mn = min(m, n)
    d = Vector{T}(undef, mn)
    e = Vector{T}(undef, mn-1)
    Q = Matrix{T}(undef, 0, 0)
    Pt = Matrix{T}(undef, 0, 0)
    C = Matrix{T}(undef, 0, 0)
    work = Vector{T}(undef, 2*max(m, n))
    gbbrd!('N', m, n, 0, bandwidth(A, 1), bandwidth(A, 2), bandeddata(A), d, e, Q, Pt, C, work)
    Bidiagonal(d, e, :U)
end

bidiagonalize!(A::AbstractMatrix) = _bidiagonalize!(A, MemoryLayout(typeof(A)))
bidiagonalize(A::AbstractMatrix) = bidiagonalize!(copy(A))

svdvals!(A::BandedMatrix) = svdvals!(bidiagonalize!(A))
svdvals(A::BandedMatrix) = svdvals!(copy(A))

# We turn off bound checking to allow nicer syntax without branching
#setindex!(A::BandedMatrix,v,k::Integer,j::Integer)=((A.l≤j-k≤A.u)&&k≤A.n)?ussetindex!(A,v,k,j):throw(BoundsError())
#setindex!(A::BandedMatrix,v,kr::AbstractRange,j::Integer)=(A.l≤j-kr[end]≤j-kr[1]≤A.u&&kr[end]≤A.n)?ussetindex!(A,v,kr,j):throw(BoundsError())

## ALgebra and other functions

function fill!(A::BandedMatrix{T}, x) where T
    x == zero(T) || throw(BandError(A))
    fill!(A.data, x)
    A
end

function diag(A::BandedMatrix{T}) where {T}
    n=size(A,1)
    @assert n==size(A,2)

    vec(A.data[A.u+1,1:n])
end

#implements fliplr(flipud(A))
function fliplrud(A::BandedMatrix)
    n,m=size(A)
    l=A.u+n-m
    u=A.l+m-n
    ret=BandedMatrix(eltype(A),n,m,l,u)
    for j = 1:size(ret,2), k = colrange(ret,j)
        @inbounds ret[k,j] = A[n-k+1,m-j+1]
    end
    ret
end

for OP in (:real, :imag)
    @eval $OP(A::BandedMatrix) = _BandedMatrix($OP(A.data),A.raxis,A.l,A.u)
end

## BandedSubBandedMatrix routines

# gives the band which is diagonal for the parent
bandshift(a::AbstractRange, b::AbstractRange) = first(a)-first(b)
bandshift(::Slice{OneTo{Int}}, b::AbstractRange) = 1-first(b)
bandshift(a::AbstractRange, ::Slice{OneTo{Int}}) = first(a)-1
bandshift(::Slice{OneTo{Int}}, b::Slice{OneTo{Int}}) = 0
bandshift(S) = bandshift(parentindices(S)[1],parentindices(S)[2])

# BandedMatrix with unit range indexes is also banded
const BandedSubBandedMatrix{T, C, R, I1<:AbstractUnitRange, I2<:AbstractUnitRange, t} =
    SubArray{T,2,BandedMatrix{T, C, R},Tuple{I1,I2},t}

isbanded(::BandedSubBandedMatrix) = true
sublayout(::BandedColumns{L}, ::Type{<:Tuple{AbstractUnitRange,J}}) where {L,J<:AbstractUnitRange} =
    bandedcolumns(sublayout(L(),Tuple{Slice{OneTo{Int}},J}))
BroadcastStyle(::Type{<:BandedSubBandedMatrix}) = BandedStyle()

function _shift(bm::BandedSubBandedMatrix)
    kr,jr=parentindices(bm)
    kr[1]-jr[1]
end

function similar(bm::BandedSubBandedMatrix, T::Type=eltype(bm),
                      n::Integer=size(bm, 1), m::Integer=size(bm, 2),
                      l::Integer=max(0, bandwidth(parent(bm), 1) - _shift(bm)),
                      u::Integer=max(0, bandwidth(parent(bm), 2) + _shift(bm)))
    similar(bm.parent, T, n, m, l, u)
end

bandeddata(V::SubArray) = view(bandeddata(parent(V)), :, parentindices(V)[2])

bandwidths(S::SubArray{T,2,<:AbstractMatrix,I}) where {T,I<:Tuple{Vararg{AbstractUnitRange}}} =
    bandwidths(parent(S)) .+ (-1,1) .* bandshift(S)

if VERSION < v"1.2-"
    @inline function inbands_getindex(S::BandedSubBandedMatrix{T}, k::Integer, j::Integer) where {T}
        @inbounds r = inbands_getindex(parent(S), reindex(S, parentindices(S), (k, j))...)
        r
    end

    @inline function inbands_setindex!(S::BandedSubBandedMatrix{T}, v, k::Integer, j::Integer) where {T}
        @inbounds r = inbands_setindex!(parent(S), v, reindex(S, parentindices(S), (k, j))...)
        r
    end
else
    @inline function inbands_getindex(S::BandedSubBandedMatrix{T}, k::Integer, j::Integer) where {T}
        @inbounds r = inbands_getindex(parent(S), reindex(parentindices(S), (k, j))...)
        r
    end

    @inline function inbands_setindex!(S::BandedSubBandedMatrix{T}, v, k::Integer, j::Integer) where {T}
        @inbounds r = inbands_setindex!(parent(S), v, reindex(parentindices(S), (k, j))...)
        r
    end
end

function convert(::Type{BandedMatrix}, S::BandedSubBandedMatrix)
    A=parent(S)
    kr,jr=parentindices(S)
    shft=kr[1]-jr[1]
    l,u=bandwidths(A)
    if -u ≤ shft ≤ l
        data = A.data[:, jr]
    elseif shft > l
        # need to add extra zeros at top since negative bandwidths not supported
        # new bandwidths = (0,u+shft)
        data = fill!(similar(A.data,u+shft+1,length(jr)), zero(eltype(A.data)))
        data[1:l+u+1,:] = A.data[:,jr]
    else  # shft < -u
        data = fill!(similar(A.data,l-shft+1,length(jr)), zero(eltype(A.data)))
        data[-shft-u+1:end,:] = A.data[:,jr]  # l-shft+1 - (-shft-u) == l+u+1
    end
    _BandedMatrix(data,length(kr),max(0, l-shft),max(0, u+shft))
end


text:
module BilliardModels

VERSION < v"0.4-" && using Docile

include("Vector2d.jl")
include("BilliardModel.jl")

end

ext: jl

max_stars_repo_name: dpsanders/BilliardModel

max_issues_repo_name: dpsanders/BilliardModel

max_forks_repo_name: dpsanders/BilliardModel


text:
using Documenter
using Altro

makedocs(
    sitename = "Altro",
    format = Documenter.HTML(prettyurls = false),
    pages = [
        "Introduction" => "index.md",
        "Getting Started" => "quickstart.md",
        "Solver Options" => "options.md",
        "Advanced Options" => "advanced.md",
        "API" => "api.md"
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/RoboticExplorationLab/Altro.jl.git",
)


text:
using DiffEqBase: remake, solve, ODEProblem, EnsembleProblem, EnsembleThreads

####################################################################################################
# this function takes into accound a parameter passed to the vector field
# Putting the options `save_start = false` seems to give bugs with Sundials
function flowTimeSol(x, p, tm, pb::ODEProblem; alg = Euler(), kwargs...)
	_prob = remake(pb; u0 = x, tspan = (zero(eltype(tm)), tm), p = p)
	# the use of concrete_solve makes it compatible with Zygote
	sol = solve(_prob, alg; save_everystep = false, kwargs...)
	return (t = sol.t[end], u = sol[end])
end

# this function is a bit different from the previous one as it is geared towards parallel computing of the flow.
function flowTimeSol(x::AbstractMatrix, p, tm, epb::EnsembleProblem; alg = Euler(), kwargs...)
	# modify the function which asigns new initial conditions
	# see docs at https://docs.sciml.ai/dev/features/ensemble/#Performing-an-Ensemble-Simulation-1
	_prob_func = (prob, ii, repeat) -> prob = remake(prob, u0 = x[:, ii], tspan = (zero(eltype(tm[ii])), tm[ii]), p = p)
	_epb = setproperties(epb, output_func = (sol, i) -> ((t = sol.t[end], u = sol[end]), false), prob_func = _prob_func)
	sol = solve(_epb, alg, EnsembleThreads(); trajectories = size(x, 2), save_everystep = false, kwargs...)
	# sol.u contains a vector of tuples (sol_i.t[end], sol_i[end])
	return sol.u
end

flow(x, p, tm, pb::ODEProblem; alg = Euler(), kwargs...) = flowTimeSol(x, p, tm, pb; alg = alg, kwargs...).u
flow(x, p, tm, pb::EnsembleProblem; alg = Euler(), kwargs...) = flowTimeSol(x, p, tm, pb; alg = alg, kwargs...)
flow(x, tm, pb::Union{ODEProblem, EnsembleProblem}; alg = Euler(), kwargs...) =  flow(x, nothing, tm, pb; alg = alg, kwargs...)
####################################################################################################
# function used to compute the derivative of the flow, so pb encodes the variational equation
function dflow(x::AbstractVector, p, dx, tm, pb::ODEProblem; alg = Euler(), kwargs...)
	n = length(x)
	_prob = remake(pb; u0 = vcat(x, dx), tspan = (zero(eltype(tm)), tm), p = p)
	# the use of concrete_solve makes it compatible with Zygote
	sol = solve(_prob, alg, save_everystep = false; kwargs...)[end]
	return (t = tm, u = sol[1:n], du = sol[n+1:end])
end

# same for Parallel computing
function dflow(x::AbstractMatrix, p, dx, tm, epb::EnsembleProblem; alg = Euler(), kwargs...)
	N = size(x,1)
	_prob_func = (prob, ii, repeat) -> prob = remake(prob, u0 = vcat(x[:, ii], dx[:, ii]), tspan = (zero(eltype(tm[ii])), tm[ii]), p = p)
	_epb = setproperties(epb, output_func = (sol,i) -> ((t = sol.t[end], u = sol[end][1:N], du = sol[end][N+1:end]), false), prob_func = _prob_func)
	sol = solve(_epb, alg, EnsembleThreads(); trajectories = size(x, 2), save_everystep = false, kwargs...)
	return sol.u
end

dflow(x, dx, tspan, pb::Union{ODEProblem, EnsembleProblem}; alg = Euler(), kwargs...) = dflow(x, nothing, dx, tspan, pb; alg = alg, kwargs...)
####################################################################################################
# this function takes into accound a parameter passed to the vector field
function dflow_fd(x, p, dx, tm, pb::ODEProblem; alg = Euler(), δ = 1e-9, kwargs...)
	sol1 = flow(x .+ δ .* dx, p, tm, pb; alg = alg, kwargs...)
	sol2 = flow(x 			, p, tm, pb; alg = alg, kwargs...)
	return (t = tm, u = sol2, du = (sol1 .- sol2) ./ δ)
end

function dflow_fd(x, p, dx, tm, pb::EnsembleProblem; alg = Euler(), δ = 1e-9, kwargs...)
	sol1 = flow(x .+ δ .* dx, p, tm, pb; alg = alg, kwargs...)
	sol2 = flow(x 			, p, tm, pb; alg = alg, kwargs...)
	return [(t = sol1[ii][1], u = sol2[ii][2], du = (sol1[ii][2] .- sol2[ii][2]) ./ δ) for ii = 1:size(x,2) ]
end
dflow_fd(x, dx, tm, pb::Union{ODEProblem, EnsembleProblem}; alg = Euler(), δ = 1e-9, kwargs...) = dflow_fd(x, nothing, dx, tm, pb; alg = alg, δ = δ, kwargs...)
####################################################################################################
# this gives access to the full solution, convenient for Poincaré shooting
# this function takes into accound a parameter passed to the vector field and returns the full solution from the ODE solver. This is useful in Poincare Shooting to extract the period.
function flowFull(x, p, tm, pb::ODEProblem; alg = Euler(), kwargs...)
	_prob = remake(pb; u0 = x, tspan = (zero(tm), tm), p = p)
	sol = solve(_prob, alg; kwargs...)
end

function flowFull(x, p, tm, epb::EnsembleProblem; alg = Euler(), kwargs...)
	_prob_func = (prob, ii, repeat) -> prob = remake(prob, u0 = x[:, ii], tspan = (zero(eltype(tm[ii])), tm[ii]), p = p)
	_epb = setproperties(epb, prob_func = _prob_func)
	sol = solve(_epb, alg, EnsembleThreads(); trajectories = size(x, 2), kwargs...)
end
flowFull(x, tm, pb::Union{ODEProblem, EnsembleProblem}; alg = Euler(), kwargs...) = flowFull(x, nothing, tm, pb; alg = alg, kwargs...)
####################################################################################################
# Structures related to computing ODE/PDE Flows
"""
$(TYPEDEF)
$(TYPEDFIELDS)

# Simplified constructor(s)
We provide a simple constructor where you only pass the vector field `F`, the flow `ϕ` and its differential `dϕ`:

	fl = Flow(F, ϕ, dϕ)

# Simplified constructors for DifferentialEquations.jl

There are some simple constructors for which you only have to pass a `prob::ODEProblem` or `prob::EnsembleProblem` (for parallel computation) from `DifferentialEquations.jl` and an ODE time stepper like `Tsit5()`. Hence, you can do for example

	fl = Flow(F, prob, Tsit5(); kwargs...)

where `kwargs` is passed to `DiffEqBase::solve`. If your vector field depends on parameters `p`, you can define a `Flow` using

	fl = Flow(F, p, prob, Tsit5(); kwargs...)

Finally, you can pass two `ODEProblem` where the second one is used to compute the variational equation:

	fl = Flow(F, p, prob1::ODEProblem, alg1, prob2::ODEProblem, alg2; kwargs...)

"""
@with_kw struct Flow{TF, Tf, Tts, Tff, Td, Tse, Tprob, TprobMono, Tfs}
	"The vector field `(x, p) -> F(x, p)` associated to a Cauchy problem,"
	F::TF = nothing

	"The flow (or semigroup) associated to the Cauchy problem `(x, p, t) -> flow(x, p, t)`. Only the last time point must be returned."
	flow::Tf = nothing

	"Flow which returns the tuple (t, u(t)). Optional, mainly used for plotting on the user side. Please use `nothing` as default."
	flowTimeSol::Tts = nothing

	"The flow (or semigroup) associated to the Cauchy problem `(x, p, t) -> flow(x, p, t)`. The whole solution on the time interval [0,t] must be returned. It is not strictly necessary to provide this, mainly used for plotting on the user side. Please use `nothing` as default."
	flowFull::Tff = nothing

	"The differential `dflow` of the flow w.r.t. `x`, `(x, p, dx, t) -> dflow(x, p, dx, t)`. One important thing is that we require `dflow(x, dx, t)` to return a Named Tuple: `(t = t, u = flow(x, p, t), du = dflow(x, p, dx, t))`, the last component being the value of the derivative of the flow."
	dflow::Td = nothing

	"Serial version of dflow. Used internally when using parallel multiple shooting. Please use `nothing` as default."
	dfSerial::Tse = nothing

	"[Internal] store the ODEProblem associated to the flow of the Cauchy problem"
	prob::Tprob = nothing

	"[Internal] store the ODEProblem associated to the flow of the variational problem"
	probMono::TprobMono = nothing

	"[Internal] Serial version of the flow"
	flowSerial::Tfs = nothing
end

# constructors
Flow(F, fl, df = nothing) = Flow(F = F, flow = fl, dflow = df, dfSerial = df)

# callable struct
(fl::Flow)(x, p, t; k...)     			  			= fl.flow(x, p, t; k...)
(fl::Flow)(x, p, dx, t; k...) 	  					= fl.dflow(x, p, dx, t; k...)
(fl::Flow)(::Val{:Full}, x, p, t; k...) 	  		= fl.flowFull(x, p, t; k...)
(fl::Flow)(::Val{:TimeSol}, x, p, t; k...)  		= fl.flowTimeSol(x, p, t; k...)
(fl::Flow)(::Val{:SerialTimeSol}, x, p, t; k...)   	= fl.flowSerial(x, p, t; k...)
(fl::Flow)(::Val{:SerialdFlow}, x, p, dx, t; k...)  = fl.dfSerial(x, p, dx, t; k...)

"""
Creates a Flow variable based on a `prob::ODEProblem` and ODE solver `alg`. The vector field `F` has to be passed, this will be resolved in the future as it can be recovered from `prob`. Also, the derivative of the flow is estimated with finite differences.
"""
# this constructor takes into accound a parameter passed to the vector field
function Flow(F, p, prob::Union{ODEProblem, EnsembleProblem}, alg; kwargs...)
	probserial = prob isa EnsembleProblem ? prob.prob : prob
	return Flow(F = F,
		flow = (x, p, t; kw2...) -> flow(x, p, t, prob; alg = alg, kwargs..., kw2...),

		flowTimeSol = (x, p, t; kw2...) -> flowTimeSol(x, p, t, prob; alg = alg, kwargs..., kw2...),

		flowFull = (x, p, t; kw2...) -> flowFull(x, p, t, prob; alg = alg, kwargs..., kw2...),

		dflow = (x, p, dx, t; kw2...) -> dflow_fd(x, p, dx, t, prob; alg = alg, kwargs..., kw2...),

		# serial version of dflow. Used for the computation of Floquet coefficients
		dfSerial = (x, p, dx, t; kw2...) -> dflow_fd(x, p, dx, t, probserial; alg = alg, kwargs..., kw2...),

		flowSerial = (x, p, t; kw2...) -> flowTimeSol(x, p, t, probserial; alg = alg, kwargs..., kw2...),

		prob = prob, probMono = nothing,
		)
end

function Flow(F, p, prob1::Union{ODEProblem, EnsembleProblem}, alg1, prob2::Union{ODEProblem, EnsembleProblem}, alg2; kwargs...)
	probserial1 = prob1 isa EnsembleProblem ? prob1.prob : prob1
	probserial2 = prob2 isa EnsembleProblem ? prob2.prob : prob2
	return Flow(F = F,
		flow = (x, p, t; kw2...) -> flow(x, p, t, prob1, alg = alg1; kwargs..., kw2...),

		flowTimeSol = (x, p, t; kw2...) -> flowTimeSol(x, p, t, prob1; alg = alg1, kwargs..., kw2...),

		flowFull = (x, p, t; kw2...) -> flowFull(x, p, t, prob1, alg = alg1; kwargs..., kw2...),

		dflow = (x, p, dx, t; kw2...) -> dflow(x, p, dx, t, prob2; alg = alg2, kwargs..., kw2...),

		# serial version of dflow. Used for the computation of Floquet coefficients
		dfSerial = (x, p, dx, t; kw2...) -> dflow(x, p, dx, t, probserial2; alg = alg2, kwargs..., kw2...),

		flowSerial = (x, p, t; kw2...) -> flowTimeSol(x, p, t, probserial1; alg = alg1, kwargs..., kw2...),

		prob = prob1, probMono = prob2,
		)
end


text:
# Autogenerated wrapper script for S2Geometry_jll for x86_64-apple-darwin
export libs2

using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("S2Geometry")
JLLWrappers.@declare_library_product(libs2, "@rpath/libs2.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libs2,
        "lib/libs2.0.10.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Balloons

export deflate, inflate!, inflate

abstract type CompressionAlgorithm end

"""
    Packbits Compression

A simple byte-oriented run-length encoding as defined by the TIFF Spec 6.0
"""
struct Packbits <: CompressionAlgorithm end

inflate!(pb::Packbits, input, output) = inflate!(pb, reinterpret(UInt8, input), reinterpret(UInt8, output))

function inflate!(::Packbits, input::AbstractVector{UInt8}, output::AbstractVector{UInt8})
    ipos, opos = 1, 1

    while opos <= length(output)
        n = reinterpret(Int8, input[ipos])
        ipos += 1
        if 0 <= n <= 127
            output[opos:opos+n-1] = view(input, ipos:ipos+n-1)
            ipos += n
            opos += n
        elseif -127 <= n <= -1
            nxt = input[ipos]
            output[opos:(opos-n-1)] .= nxt
            opos += -n
            ipos += 1
        end
    end
end

"""
    Lempel-Ziv-Welch Compression

A robust universal lossless compression algorithm. This implementation supports
variable-width compression codes.
"""
struct LZW <: CompressionAlgorithm
    initial_codewidth::Int
    initial_dictsize::Int

    clearcode::UInt32
    endcode::UInt32
end

LZW() = LZW(9, 256, 256, 257)

mask(t::Type{T}, p::Integer, l::Integer) where {T} = (typemax(T) << (sizeof(T)*8-l)) >> p

extract(packed::T, p::Integer, l::Integer) where {T} = (packed & mask(T, p, l)) >> (sizeof(T)*8 - p - l)

"""
    stuff(v, a, loc)

Stuff bits of `a` into 32-bit unsigned integer `v` at a given byte location `loc`

```jldoctest
julia> Balloons.stuff(zero(UInt32), 0xff, 1)
0xff000000

julia> Balloons.stuff(zero(UInt32), 0xff, 4)
0x000000ff
```
"""
function stuff(v::UInt32, a::UInt8, loc::T) where {T <: Integer}
    @assert 0 < loc <= 4
    v += (zero(UInt32) + a) << ((4-loc)*8)
end

function initdict!(lzw::LZW, dict)
    empty!(dict)

    for i in 0:(lzw.initial_dictsize - 1)
        dict[i] = UInt8[i]
    end
    dict[lzw.clearcode] = [0x00]
    dict[lzw.endcode] = [0x00]
end

function nextchar(input, ρ, bitwidth)
    i = ρ÷8+1
    v = stuff(zero(UInt32), input[i], 1)

    for j in 1:min(3, length(input)-i)
        v = stuff(v, input[i+j], j+1)
    end

    extract(v, mod(ρ, 8), bitwidth)
end

function inflate(lzw::LZW, input::Vector{UInt8})
    bitwidth = lzw.initial_codewidth

    ρ = 0 # position in bits

    output = IOBuffer()

    dict = Dict{UInt32, Vector{UInt8}}()
    k = lzw.clearcode

    entry = UInt8[]
    prev_entry = UInt8[]

    while true
        # increase bitwidth if we're out of space
        if length(dict) >= 2^bitwidth - 1
            bitwidth += 1
        end

        k = nextchar(input, ρ, bitwidth)
        if k == lzw.endcode
            break
        end

        # reset dictionary and variables on detecting a clear code
        if k == lzw.clearcode
            ρ += bitwidth

            # reset everything back to baseline
            initdict!(lzw, dict)
            sz = length(dict)
            bitwidth = lzw.initial_codewidth

            # get next character
            k = nextchar(input, ρ, bitwidth)
            (k == lzw.endcode) && break
            entry = dict[k]
            write(output, entry)
            prev_entry = entry
        else
            if haskey(dict, k)
                entry = dict[k]
            else
                entry = vcat(prev_entry, first(prev_entry))
            end

            write(output, entry)
            if length(prev_entry) != 0
                dict[length(dict)] = vcat(prev_entry, first(entry))
            end
            prev_entry = entry
        end

        ρ += bitwidth
    end

    seekstart(output)
    readavailable(output)
end

end


text:
struct DAESolution{T,N,uType,duType,uType2,DType,tType,P,A,ID} <: AbstractDAESolution{T,N}
  u::uType
  du::duType
  u_analytic::uType2
  errors::DType
  t::tType
  prob::P
  alg::A
  interp::ID
  dense::Bool
  tslocation::Int
  destats::DEStats
  retcode::Symbol
end
(sol::DAESolution)(t,deriv::Type=Val{0};idxs=nothing,continuity=:left) = sol.interp(t,idxs,deriv,sol.prob.p,continuity)
(sol::DAESolution)(v,t,deriv::Type=Val{0};idxs=nothing,continuity=:left) = sol.interp(v,t,idxs,deriv,sol.prob.p,continuity)

function build_solution(
        prob::AbstractDAEProblem{uType,duType,tType,isinplace},
alg,t,u;dense=false,du=[],
interp = !isempty(du) ? HermiteInterpolation(t,u,du) : LinearInterpolation(t,u),
timeseries_errors=true,dense_errors=true, retcode = :Default, destats=DEStats(), kwargs...) where {uType,duType,tType,isinplace}

  T = eltype(eltype(u))
  N = length((size(u[1])..., length(u)))

  if has_analytic(prob.f)
    u_analytic = Vector{typeof(prob.u0)}()
    for i in 1:size(u,1)
      push!(u_analytic,prob.analytic(prob.u0,t[i]))
    end

    save_everystep = length(u) > 2

    errors = Dict{Symbol,real(eltype(prob.u0))}()
    if !isempty(u_analytic)
      errors[:final] = norm(recursive_mean(abs.(u[end]-u_analytic[end])))

      if save_everystep && timeseries_errors
        errors[:l∞] = norm(maximum(vecvecapply((x)->abs.(x),u-u_analytic)))
        errors[:l2] = norm(sqrt(recursive_mean(vecvecapply((x)->float.(x).^2,u-u_analytic))))
        if dense && dense_errors
          densetimes = collect(range(t[1], stop=t[end], length=100))
          interp_u = interp(densetimes)
          interp_analytic = [prob.analytic(t,u[1]) for t in densetimes]
          errors[:L∞] = norm(maximum(vecvecapply((x)->abs.(x),interp_u-interp_analytic)))
          errors[:L2] = norm(sqrt(recursive_mean(vecvecapply((x)->float.(x).^2,interp_u-interp_analytic))))
        end
      end
    end
    DAESolution{T,N,typeof(u),typeof(du),typeof(u_analytic),typeof(errors),typeof(t),
                       typeof(prob),typeof(alg),typeof(interp)}(u,du,u_analytic,errors,t,prob,alg,interp,dense,0,destats,retcode)
  else
    DAESolution{T,N,typeof(u),typeof(du),Nothing,Nothing,typeof(t),
                       typeof(prob),typeof(alg),typeof(interp)}(u,du,nothing,nothing,t,prob,alg,interp,dense,0,destats,retcode)
  end
end

function build_solution(sol::AbstractDAESolution{T,N},u_analytic,errors) where {T,N}
  DAESolution{T,N,typeof(sol.u),typeof(sol.du),typeof(u_analytic),typeof(errors),typeof(sol.t),
                     typeof(sol.prob),typeof(sol.alg),typeof(sol.interp)}(
                     sol.u,sol.du,u_analytic,errors,sol.t,
              sol.prob,sol.alg,sol.interp,sol.dense,sol.tslocation,sol.destats,sol.retcode)
end

function solution_new_retcode(sol::AbstractDAESolution{T,N},retcode) where {T,N}
  DAESolution{T,N,typeof(sol.u),typeof(sol.du),typeof(sol.u_analytic),
              typeof(sol.errors),typeof(sol.t),
              typeof(sol.prob),typeof(sol.alg),typeof(sol.interp)}(
              sol.u,sol.du,sol.u_analytic,sol.errors,sol.t,
              sol.prob,sol.alg,sol.interp,sol.dense,sol.tslocation,sol.destats,retcode)
end

function solution_new_tslocation(sol::AbstractDAESolution{T,N},tslocation) where {T,N}
  DAESolution{T,N,typeof(sol.u),typeof(sol.du),typeof(sol.u_analytic),
              typeof(sol.errors),typeof(sol.t),
              typeof(sol.prob),typeof(sol.alg),typeof(sol.interp)}(
              sol.u,sol.du,sol.u_analytic,sol.errors,sol.t,
              sol.prob,sol.alg,sol.interp,sol.dense,tslocation,sol.destats,sol.retcode)
end

function solution_slice(sol::AbstractDAESolution{T,N},I) where {T,N}
  DAESolution{T,N,typeof(sol.u),typeof(sol.du),typeof(sol.u_analytic),
              typeof(sol.errors),typeof(sol.t),
              typeof(sol.prob),typeof(sol.alg),typeof(sol.interp)}(
              sol.u[I],sol.du[I],
              sol.u_analytic === nothing ? nothing : sol.u_analytic[I],
              sol.errors,sol.t[I],
              sol.prob,sol.alg,
              sol.interp,false,
              sol.tslocation,sol.destats,sol.retcode)
end


text:
using RPG

RPG.julia_main([""])


text:
function phi(z,a,interval,x)
    z = [z; -z];
    a = [a; -a];

    r = (interval[2] - interval[1]) / 2;
    Emid = interval[1] + r;
    z = z*r .+ Emid;
    a = a*r;

    y = zeros(size(x)...);
    for i = 1:length(z)
        y = y + a[i] ./ (x .- z[i]) + conj(a[i]) ./ (x .- conj(z[i]));
    end
    y
end


text:
__precompile__()

module Wavelets

include("mod/Util.jl")
include("mod/WT.jl")
include("mod/Transforms.jl")
include("mod/Threshold.jl")
include("mod/Plot.jl")

using Reexport
@reexport using .Util, .WT, .Transforms, .Threshold, .Plot

end


text:
using Test

using Interact

@test gettheme() == Interact.InteractBulma.Bulma()
settheme!(:nativehtml)
@test gettheme() == Interact.InteractBase.NativeHTML()
settheme!(:bulma)
@test gettheme() == Interact.InteractBulma.Bulma()


text:
using Parameters
using DifferentialEquations
using PyPlot

include("cr_model.jl")

## CR exp resource, type I consumer
let
    u0 = [0.5, 0.5]
    t_span = (0.0, 50.0)
    p = CRPar()

    prob = ODEProblem(cr_expI_model, u0, t_span, p)
    sol = DifferentialEquations.solve(prob, reltol = 1e-8)

    cr_expI_ts = figure()
    plot(sol.t, sol.u)
    xlabel("time")
    ylabel("Density")
    return cr_expI_ts
end
## CR log resource, type I consumer
let
    u0 = [0.5, 0.5]
    t_span = (0.0, 50.0)
    p = CRPar()

    prob = ODEProblem(cr_logI_model, u0, t_span, p)
    sol = DifferentialEquations.solve(prob, reltol = 1e-8)

    cr_logI_ts = figure()
    plot(sol.t, sol.u)
    xlabel("time")
    ylabel("Density")
    return cr_logI_ts
end

## CR log resource, type II consumer
let
    u0 = [0.5, 0.5]
    t_span = (0.0, 75.0)
    p = CRPar(e=0.5)

    prob = ODEProblem(cr_logII_model, u0, t_span, p)
    sol = DifferentialEquations.solve(prob, reltol = 1e-8)
    cr_logII_ts = figure()
    plot(sol.t, sol.u)
    xlabel("time")
    ylabel("Density")
    return cr_logII_ts
end

## CR exp resource, type II consumer

let
    u0 = [0.5,0.5]
    t_span = (0.0, 50.0)
    p = CRPar(r=0.5)

    prob = ODEProblem(cr_expII_model, u0, t_span, p)
    sol = DifferentialEquations.solve(prob, reltol = 1e-8)

    cr_expII_ts = figure()
    plot(sol.t, sol.u)
    xlabel("time")
    ylabel("Density")
    return cr_expII_ts
end

let
    Rrange = 0.0:0.1:30.0
    u0 = [0.5,0.5]
    t_span = (0.0, 50.0)
    p = CRPar(r=0.549)

    prob = ODEProblem(cr_expII_model, u0, t_span, p)
    sol = DifferentialEquations.solve(prob, reltol = 1e-8)

    cr_expII_phasets = figure()
    plot(collect(Rrange), [expII_resiso(R, p) for R in Rrange], color="blue")
    vlines(expII_coniso(p), 0.0, 30.0, colors="orange")
    plot(sol[1,:], sol[2,:])
    xlim(0.0,30.0)
    ylim(0.0,30.0)
    xlabel("Resource")
    ylabel("Consumer")
    return cr_expII_phasets
end

max_stars_repo_name: McCannLab/TheoryClub

max_issues_repo_name: McCannLab/TheoryClub

max_forks_repo_name: McCannLab/TheoryClub


text:
#! /usr/local/julia-1.1.1/bin/julia

using Pkg

function recompile()
    for pkg in Pkg.installed()
        try
            pkgsym = Symbol(pkg[1])
            eval(:(using $pkgsym))
        catch
            println("ERROR!")
        end
    end
end

Pkg.add(map(x->string(x), ARGS))
recompile()

max_stars_repo_name: Losses/powerData

max_issues_repo_name: Losses/powerData

max_forks_repo_name: Losses/powerData


text:
import OpenSpiel
using ..AlphaZero

spiel_game = OpenSpiel.load_game("tic_tac_toe")

function vectorize_state(spec, state)
  if !OpenSpiel.is_terminal(state)
    obsten = OpenSpiel.observation_tensor(state)
    empt =  obsten[1:9]
    black = obsten[10:18]
    white = obsten[19:27]
    # in AlphaZero, oracle always see game from current player's point of view
    if OpenSpiel.current_player(state) == 1
      white, black = black, white
    end
    ans = convert(Array{Float32}, vcat(empt, white, black))

    return reshape(ans, 3,3,3)
  else
    return cat(
      zeros(Float32, 3,3),
      ones(Float32,3,3),
      ones(Float32,3,3);
      dims=3
    )
  end
end

#####
##### Interaction API (from AlphaZero's tictactoe)
#####

function action_string(spec, a)
  string(Char(Int('A') + a))
end

function parse_action(spec, str)
  length(str) == 1 || (return nothing)
  x = Int(uppercase(str[1])) - Int('A')
  (0 <= x < 9) ? x : nothing
end

GameSpec() = AlphaZero.OpenSpielWrapper.Spec(
  spiel_game;
  vectorize_state,
  action_string,
  parse_action,
  suppress_warnings=true
)


text:
using Test
using NSGAII
using Statistics
import Base.copy

include("./cornercases.jl")

fitness(x::FloatIndividual)=[mean(x.genes),std(x.genes),median(x.genes)]

function copy(ind::FloatIndividual)
    genes=deepcopy(ind.genes)
    fitness=zeros(length(ind.fitness))
    FloatIndividual(genes,fitness)
end

@testset "Populate" begin
    cfg=Cambrian.get_config("./test.yaml";n_population= 7,n_offsprings=5,d_fitness=3)
    e=NSGA2Evolution{FloatIndividual}(cfg,fitness)
    oldPop=copy(e.population)
    NSGA2Populate(e)
    for x in oldPop
        @test x in e.population
    end
    for x in e.population
        if x in oldPop == false
            @test e.offsprings[objectid(x)]==true
        end
    end
    @test length(e.population)==length(oldPop)+5
end

@testset "Evaluate" begin
    cfg=Cambrian.get_config("./test.yaml";n_population= 7,n_offsprings=5,d_fitness=3)
    e=NSGA2Evolution{FloatIndividual}(cfg,fitness)
    L=[]
    for i in 1:3
        e.offsprings[objectid(e.population[i])]=false
        push!(L,copy(e.population[i].fitness))
    end
    for i in 4:7
        push!(L,copy(e.population[i].fitness))
    end
    evaluate(e)
    for i in 1:3
        @test all(e.population[i].fitness.==L[i])
    end
    for i in 4:7
        @test all(e.population[i].fitness.!=L[i])
    end

end

@testset "Generation" begin
    cfg=Cambrian.get_config("./test.yaml";n_population= 7,n_offsprings=5,d_fitness=3)
    e=NSGA2Evolution{FloatIndividual}(cfg,fitness)
    evaluate(e)
    max=maxPop(e.population)
    NSGA2Generation(e)
    for x in max
        @test x in e.population
    end
end

@testset "Step" begin
    cfg=Cambrian.get_config("./test.yaml";n_population= 7,n_offsprings=5,d_fitness=3)
    e=NSGA2Evolution{FloatIndividual}(cfg,fitness)
    step!(e)
    @test length(e.population)==e.config.n_population
    step!(e)
    @test length(e.population)==e.config.n_population
    for x in e.population
        @test e.offsprings[objectid(x)]==false
    end
end


text:
# ---------------------------------
# PIECEWISE LINEAR UPPER ENVELOPE |
# ____________________ __________ |
#    Author: Gustavo Pereira      |
#   Columbia University, 2020     |
# ---------------------------------

abstract type AbstractPiecewiseLinear{T} end
@inline Base.length(fun::AbstractPiecewiseLinear)=length(get_x(fun))
@inline Base.eltype(fun::AbstractPiecewiseLinear{T}) where T=T
function transfer_point!(pl_target, i_target,
                      pl_source, i_source)
    pl_target[i_target]=pl_source[i_source]
    transfer_subordinate!(pl_target, i_target, pl_source, i_source)
end
@inline Base.setindex!(fun::AbstractPiecewiseLinear, pt, i)=begin
    set_x!(fun, get_x(pt), i)
    set_y!(fun, get_y(pt), i)
end

"""
# Piecewise linear function

    mutable struct PiecewiseLinear{T}
      xcoords::Vector{T}
      ycoords::Vector{T}
    end

The purpose is to represent "piecewise linear" functions by storing the x and y
coordinates of the points that define them. The i-th point has x coordinate
xcoords[i], and y coordinate `ycoords[i]`

Convenience methods `Base.getindex` and `Base.setindex!` are defined for this
type. Accessing index `i` returns x and y coordinates grouped in a
UpperEnvelope.Point2 type. That type should also be used for insertion.

## Assumptions
**The field `xcoords` should be ordered.** Moreover, `xcoords` and `ycoords`
must be of the same length.

"""
mutable struct PiecewiseLinear{T} <: AbstractPiecewiseLinear{T}
    xcoords::Vector{T}
    ycoords::Vector{T}
    function PiecewiseLinear(xcoords, ycoords)
        length(xcoords) != length(ycoords) && error("x and y coordinates must have the same length")
        new{eltype(xcoords)}(xcoords, ycoords)
    end
end

# These methods contribute very little, except for readability later on.
@inline Base.Tuple(fun::PiecewiseLinear)=(get_x(fun), get_y(fun))
@inline get_x(fun::PiecewiseLinear, i)=fun.xcoords[i]
@inline get_x(fun::PiecewiseLinear)=fun.xcoords
function set_x!(fun::PiecewiseLinear, x, i)
    fun.xcoords[i]=x
    nothing
end
@inline get_y(fun::PiecewiseLinear, i)=fun.ycoords[i]
@inline get_y(fun::PiecewiseLinear)=fun.ycoords
function set_y!(fun::PiecewiseLinear, y, i)
    fun.ycoords[i]=y
    nothing
end
Base.@propagate_inbounds @inline Base.getindex(fun::PiecewiseLinear, i)=Point2(get_x(fun, i), get_y(fun, i))
Base.resize!(fun::PiecewiseLinear, n)=begin
    resize!(get_x(fun), n);resize!(get_y(fun), n);
    nothing
end
function init_zeros(::Type{PiecewiseLinear{T}}, n) where T
    PiecewiseLinear(fill(zero(T), n), fill(zero(T), n))
end
@inline function transfer_subordinate!(pl_target::PiecewiseLinear, i_target,
                               pl_source::PiecewiseLinear, i_source)
    nothing
end
function transfer_subordinate_combo!(fun_target::PiecewiseLinear, i_target,
                                     fun_source::PiecewiseLinear, i_source,
                                     weight) where {N,T}
    nothing
end

"""
# Extended piecewise linear function

This type extends PiecewiseLinear by including an additional set of `y values`
(elements of `subordinate_vectors`) to be discarded/extended based on the
envelope procedure applied to `base_fun`.

Indexing returns the _main_ (x,y) pair. If `F::ExtendedPiecewiseEnvelope`,
then `F[i]` returns `F.base_fun[i]`.
"""
mutable struct ExtendedPiecewiseLinear{T, N} <: AbstractPiecewiseLinear{T}
    base_fun::PiecewiseLinear{T}
    subordinate_vectors::NTuple{N, Vector{T}}
end
@inline get_x(fun::ExtendedPiecewiseLinear)=get_x(fun.base_fun)
@inline get_x(fun::ExtendedPiecewiseLinear, i)=get_x(fun.base_fun, i)
@inline set_x!(fun::ExtendedPiecewiseLinear, x, i)=set_x!(fun.base_fun, x, i)

@inline get_y(fun::ExtendedPiecewiseLinear, i)=get_y(fun.base_fun, i)
@inline set_y!(fun::ExtendedPiecewiseLinear, y, i)=set_y!(fun.base_fun, y, i)
Base.@propagate_inbounds @inline Base.getindex(fun::ExtendedPiecewiseLinear, i)=getindex(fun.base_fun, i)
@inline Base.resize!(fun::ExtendedPiecewiseLinear, n)=begin
    resize!(fun.base_fun, n);
    map(fun.subordinate_vectors) do x
        resize!(x, n)
    end
    nothing
end
@inline Base.Tuple(fun::ExtendedPiecewiseLinear)=(Tuple(fun.base_fun)..., fun.subordinate_vectors...)

function init_zeros(::Type{ExtendedPiecewiseLinear{T, N}}, n) where {T, N}
    pl = init_zeros(PiecewiseLinear{T}, n)
    subord = ntuple(Val(N)) do _
        fill(zero(T), n)
    end
    ExtendedPiecewiseLinear{T,N}(pl, subord)
end
@inline function transfer_subordinate!(pl_target::ExtendedPiecewiseLinear{T, N}, i_target,
                               pl_source::ExtendedPiecewiseLinear{T, N}, i_source) where {N, T}
    @inbounds for j in 1:N
        pl_target.subordinate_vectors[j][i_target] = pl_source.subordinate_vectors[j][i_source]
    end
    nothing
end

@inline function transfer_subordinate_combo!(fun_target::ExtendedPiecewiseLinear{T,N}, i_target,
                                     fun_source::ExtendedPiecewiseLinear{T,N}, i_source,
                                     weight) where {N,T}
    @assert i_source > 1

    @inbounds for j in 1:N
        z = fun_source.subordinate_vectors[j][i_source-1]*weight
        z += fun_source.subordinate_vectors[j][i_source]*(1-weight)
        fun_target.subordinate_vectors[j][i_target] = z
    end
end

"""

Created with sole purpose of overloading base arithmetics for 2-Tuples without
messing with Tuples.

The benefit is simply notational convenience.

"""
struct Point2{T}
    x::T
    y::T
end
@inline get_x(pt::Point2)=pt.x
@inline get_y(pt::Point2)=pt.y

Base.Tuple(u::Point2)=(u.x, u.y)
Base.:+(u::Point2, v::Point2)=Point2(u.x+v.x, u.y+v.y)
Base.:*(u::Point2, k::Number)=Point2(k*u.x, k*u.y)
Base.:*(k::Number, u::Point2)=Base.:*(u, k)
Base.:-(u::Point2, v::Point2)=Point2(u.x - v.x, u.y - v.y)

" 90deg counter clockwise rotation "
rot90cc(u::Point2)=Point2(-u.y, u.x) # rotate 90 deg counter clockwise

" Dot product of `Point2` types."
pdot(u::Point2, v::Point2)=u.x*v.x+u.y*v.y

" Rotate the second argument and take dot product with the first one. "
inner_rotate(v::Point2, w::Point2)=pdot(v, rot90cc(w))

# * Main function

"""
# Compute upper envelope of piecewise linear functions

    compute_envelope(pl_funcs::NTuple{2, PiecewiseLinear})

Computes the upper envelope of two piecewise linear functions, as represented by
type PiecewiseLinear.

## Algorithm

Modulo dealing with the non-intersecting part of domains, the idea is to fix a
pair of points (i.e., a segment) in one of the functions -- call it "benchmark"
-- and analyse all segments of the other function -- call it "comparison" --
within that segment. It is easy to check whether the subsegments are "dominated"
or not, and to compute intersections.

When all comparison segments are analysed, switch the states between functions:
the benchmark function becomes the comparison one and vice-versa. If this step
is ignored, intersections **will be** overlooked.

## Notation

- `k` denotes the current state. It is either 1 or 2.
- `other(k)` returns 2 if `k == 1`, or 1 if `k==2`.

"""
function compute_envelope(pl_funcs::NTuple{2, T}, trim_empty=true) where T <: AbstractPiecewiseLinear
    envsize_max = length(pl_funcs[1]) + length(pl_funcs[2]) +
        max(length(pl_funcs[1]), length(pl_funcs[2]))

    envelope = init_zeros(T, envsize_max)
    #  ^ Initialize an instance of piecewise linear with size envsize_max, filled with NaN

    @inbounds k = determine_first(pl_funcs)
    v_idx_visit = [1, 1]
    i_env = 1

    # Deal the left part of (possibly) non-intersecting domains
    while is_lower_x(pl_funcs, k, v_idx_visit)
        if v_idx_visit[k] < length(pl_funcs[k])
            this_ind = v_idx_visit[k]

            transfer_point!(envelope, i_env, pl_funcs[k], this_ind); i_env += 1
            # pt = pl_funcs[k][this_ind]
            # envelope[i_env]=pt;

            inc_counters!(v_idx_visit, k)
        else
            break
        end
    end

    # Equal values of x in the beginning ('regular' points)
    n_iter=1
    reached_maxlen = false
    while  !reached_maxlen && (get_x(pl_funcs[1], v_idx_visit[1]) == get_x(pl_funcs[2], v_idx_visit[2]))

        fun_better = get_y(pl_funcs[1], v_idx_visit[1]) >
            get_y(pl_funcs[2], v_idx_visit[2]) ? 1 : 2

        # envelope[i_env] = ifelse(first_is_better, pl_funcs[1][v_idx_visit[1]],
        #                          pl_funcs[2][v_idx_visit[2]])
        transfer_point!(envelope, i_env, pl_funcs[fun_better][v_idx_visit[fun_better]])
        i_env += 1

        reached_maxlen = v_idx_visit[1] == length(pl_funcs[1]) || v_idx_visit[2] == length(pl_funcs[2])

        if !reached_maxlen
            # i_env += segment_intersect!(envelope, i_env,
            #                             pl_funcs[1][n_iter], pl_funcs[1][n_iter+1],
            #                             pl_funcs[2][n_iter], pl_funcs[2][n_iter+1])
            i_env += segment_intersect!(envelope, i_env, pl_funcs, 1, n_iter+1, n_iter+1)
            inc_counters!(v_idx_visit)
        end
        n_iter += 1
    end

    if reached_maxlen
        return envelope
    end

    # Irregular x values
    @inbounds k = determine_first(pl_funcs, v_idx_visit)
    @inbounds while (v_idx_visit[1] <= length(pl_funcs[1])) && (v_idx_visit[2] <= length(pl_funcs[2]))
        i_curr = v_idx_visit[k]
        i_other = v_idx_visit[other(k)]

        # First step: compute the "score" of the current point under analysis
        # (will: (I) determine whether point is included, (II) be used in
        # computing a potential intersection)
        pt_curr = pl_funcs[k][i_curr];
        pt_other = pl_funcs[other(k)][i_other]
        pt_other_prev = pl_funcs[other(k)][i_other-1]
        sco = inner_rotate(pt_curr - pt_other_prev, pt_other - pt_other_prev)

        if sco > 0
            transfer_point!(envelope, i_env, pl_funcs[k], i_curr)
            i_env += 1
        end

        if i_curr + 1 > length(pl_funcs[k])
            k = other(k)
            break
        end

        i_env += segment_intersect!(envelope, i_env, sco, pl_funcs, k, i_curr+1, i_other)
        inc_counters!(v_idx_visit, k)

        if get_x(pl_funcs[k][i_curr+1]) >= get_x(pt_other)
            k = other(k)
        end
    end

    # Now finalize by adding all rightmost nodes
    while v_idx_visit[k] <= length(pl_funcs[k])
        transfer_point!(envelope, i_env, pl_funcs[k], v_idx_visit[k]); i_env += 1
        inc_counters!(v_idx_visit, k)
    end

    # Remove NaN entries if specified
    trim_empty && resize!(envelope, i_env-1)
    return envelope
end

# Convenience method
compute_envelope(tup_pts1, tup_pts2, trim_empty=true)=begin
    Tuple(compute_envelope((construct_piecewiselinear(tup_pts1...), construct_piecewiselinear(tup_pts2...)),
                           trim_empty))
end

function construct_piecewiselinear(a1, a2, rest...)
    ExtendedPiecewiseLinear(construct_piecewiselinear(a1, a2), rest)
end

function construct_piecewiselinear(a1, a2)
    PiecewiseLinear(a1, a2)
end

# * Auxiliary functions
" Returns 1 if `i==2`, and 2 if `i==1`."
other(i)=3-i

inc_counters!(v_idx, i)=v_idx[i] += 1
inc_counters!(v_idx)=begin
    v_idx[1] += 1; v_idx[2] += 1;
end

"""
Check that the x value is lower for the current than the alternative state at
indexes given by v_idx_visit.
"""
function is_lower_x(pl_funcs, k, v_idx_visit)
    fun_other = other(k)

    @inbounds i_curr = v_idx_visit[k]
    @inbounds i_other = v_idx_visit[fun_other]
    get_x(pl_funcs[k], i_curr) < get_x(pl_funcs[fun_other], i_other)
end

Base.@propagate_inbounds function determine_first(pl_funcs, init_ind=(1, 1))
    r = 1
    if get_x(pl_funcs[2], init_ind[2]) < get_x(pl_funcs[1], init_ind[1])
        r = 2
    end
    return r
end

function segment_intersect!(pl_target, i_target, sco, pl_funcs, k, i_curr, i_other)
    ptA1 = pl_funcs[other(k)][i_other-1]
    ptA2 = pl_funcs[other(k)][i_other]
    ptB1 = pl_funcs[k][i_curr-1]
    ptB2 = pl_funcs[k][i_curr]

    P = ptA2 - ptA1
    Q = ptB2 - ptB1
    R = ptB2 - ptA1

    sco_cmp = inner_rotate(R, P)

    if sign(sco_cmp) == sign(sco)
        return false
    end

    bb = sco_cmp/(sco_cmp - sco)
    intersection_candidate = bb * ptB1 + (1-bb) * ptB2
    r = get_x(intersection_candidate) < get_x(ptA2)
    # r = bb*get_x(ptB1) + (1 - bb)*get_x(ptB2) < get_x(ptA2)

    if r
        pl_target[i_target] = intersection_candidate
        transfer_subordinate_combo!(pl_target, i_target, pl_funcs[k], i_curr, bb)
    end

    return r
end

segment_intersect!(v_store, pos_store, ptA1, ptA2, ptB1, ptB2)=begin
    segment_intersect!(v_store, pos_store,
                       inner_rotate(ptB1-ptA1, ptA2-ptA1),
                       ptA1, ptA2, ptB1, ptB2)
end

"""
    segment_intersect!(v_store, pos_store, sco, ptA1, ptA2, ptB1, ptB2)

If an intersection between pl_funcs [A1, A2] and [B1, B2] exist, insert that in
v_store at position pos_store. Also, returns Boolean for whether the
intersection exists.

Note that ptA1, ptA2, ptB1 and ptB2 are NOT general points in the plane.
Because of the problem's structure, the following inequality holds:
get_x(ptA1) <= get_x(pt_B1)

Note also that `sco` already stores a value whose sign depends on the half space
of point ptB1 relative to the line formed by the segment ptA1-ptA2 Hence, the
first way of ruling out the intersection is to compute the same value for ptB2;

Formally, `sco` is the inner product of the 90 deg counter-clockwise
rotation of `P` (as defined below) with `B1 - A1`.

"""
function segment_intersect!(v_store, pos_store, sco, ptA1::Point2, ptA2::Point2, ptB1::Point2, ptB2::Point2)

    P = ptA2 - ptA1
    Q = ptB2 - ptB1
    R = ptB2 - ptA1

    sco_cmp = inner_rotate(R, P)

    if sign(sco_cmp) == sign(sco)
        return false
    end

    # If sign switches, this means there is _possibly_ an intersection.
    #
    # Use the geometry of the problem to find the coefficient:
    bb = sco_cmp/(sco_cmp - sco)
    intersection_candidate = bb * ptB1 + (1-bb) * ptB2

    r = get_x(intersection_candidate) < get_x(ptA2)

    if r
        v_store[pos_store]=intersection_candidate
    end

    return r
end


text:
# Autogenerated wrapper script for gb_jll for armv7l-linux-musleabihf
export libgb

using GMP_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("gb")
JLLWrappers.@declare_library_product(libgb, "libgb.so.0")
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libgb,
        "lib/libgb.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    gbₕ_free(Tₐ,Tₗ,d,Dₕ₀)
    gbₕ_free(Tₐ,Tₗ,d)

Leaf boundary layer conductance for heat under **free** convection (m s-1).

# Arguments

- `Tₐ` (°C): air temperature
- `Tₗ` (°C): leaf temperature
- `d` (m): characteristic dimension, *e.g.* leaf width (see eq. 10.9 from Monteith and Unsworth, 2013).
- `Dₕ₀ = 21.5e-6`: molecular diffusivity for heat at base temperature. Use value from
[`Constants`](@ref) if not provided.

# Note

`R` and `Dₕ₀` can be found using [`Constants`](@ref). To transform in ``mol\\ m^{-2}\\ s^{-1}``,
use [`ms_to_mol`](@ref).

# References

Leuning, R., F. M. Kelliher, DGG de Pury, et E.-D. SCHULZE. 1995. « Leaf nitrogen,
photosynthesis, conductance and transpiration: scaling from leaves to canopies ». Plant,
Cell & Environment 18 (10): 1183‑1200.

Monteith, John, et Mike Unsworth. 2013. Principles of environmental physics: plants,
animals, and the atmosphere. Academic Press. Paragraph 10.1.3, eq. 10.9.
"""
function gbₕ_free(Tₐ,Tₗ,d,Dₕ₀ = Constants().Dₕ₀)
    zeroT = zero(Tₐ) # make it type stable

    if abs(Tₗ-Tₐ) > zeroT
        Gr = 1.58e8 * d^3.0 * abs(Tₗ-Tₐ) # Grashof number (Monteith and Unsworth, 2013)
        # !Note: Leuning et al. (1995) use 1.6e8 (eq. E4).
        # Leuning et al. (1995) eq. E3:
        Gbₕ_free = 0.5 * get_Dₕ(Tₐ,Dₕ₀) * (Gr^0.25) / d
    else
        Gbₕ_free = zeroT
    end

    return Gbₕ_free
end

"""
    gbₕ_forced(Wind,d)

Boundary layer conductance for heat under **forced** convection (m s-1). See eq. E1 from
Leuning et al. (1995) for more details.

# Arguments

- `Wind` (m s-1): wind speed
- `d` (m): characteristic dimension, *e.g.* leaf width (see eq. 10.9 from Monteith and Unsworth, 2013).

# Notes

`d` is the minimal dimension of the surface of an object in contact with the air.

# References

Leuning, R., F. M. Kelliher, DGG de Pury, et E.-D. SCHULZE. 1995. « Leaf nitrogen,
photosynthesis, conductance and transpiration: scaling from leaves to canopies ». Plant,
Cell & Environment 18 (10): 1183‑1200.
"""
function gbₕ_forced(Wind,d)
    0.003 * sqrt(Wind / d)
end

"""
    get_Dₕ(T,Dₕ₀)
    get_Dₕ(T)

Dₕ -molecular diffusivity for heat at base temperature- from Dₕ₀ (corrected by temperature).
See Monteith and Unsworth (2013, eq. 3.10).

# Arguments

- `Tₐ` (°C): temperature
- `Dₕ₀`: molecular diffusivity for heat at base temperature. Use value from [`Constants`](@ref)
if not provided.

# References

Monteith, John, et Mike Unsworth. 2013. Principles of environmental physics: plants,
animals, and the atmosphere. Academic Press. Paragraph 10.1.3.
"""
function get_Dₕ(T,Dₕ₀ = Constants().Dₕ₀)
    Dₕ₀ * (1 + 0.007*T)
end


text:
@enum MatchEvent noevents matchset roundinit roundwin roundtie matchend tiledraw tiledrop meldcall riichicall doraflip playerdc playerin

    function GO(str::AbstractString, pst::PlayState)

        function getbits(code)
            digits(parse(Int, code); base = 2, pad = 8)
        end

        bits = map(Bool, parsekey(getbits, "type", str))

        (   !bits[8] && !bits[6]    ) && ( lobby = 一般 )
        (    bits[8] && !bits[6]    ) && ( lobby = 上級 )
        (   !bits[8] &&  bits[6]    ) && ( lobby = 特上 )
        (    bits[8] &&  bits[6]    ) && ( lobby = 鳳凰 )

        pst.rules = Rules(!bits[1], !bits[2], !bits[3],
                bits[4], bits[5], bits[7], lobby)

        return noevents
    end

    function UN(str::AbstractString, pst::PlayState)

        namaes = String[]
        for namekey in ("n0", "n1", "n2", "n3")
            name = parsekey((s)->decodeuri(s), namekey, str)
            if !isnothing(name)
                push!(namaes, name)
            end
        end

        if length(namaes) == 1
            nameindex = findfirst(isequal(namaes[1]), pst.table.namaes) - 1
            seatindex = findfirst(isequal(Seat(nameindex)), pst.dced)
            if !isnothing(seatindex) deleteat!(pst.dced, seatindex)
            else push!(pst.dced, Seat(nameindex)) end
            return playerin
        end

        playercount = length(namaes)
        vecranks = Vector{Dan}(undef, playercount)
        vecrates = Vector{Float32}(undef, playercount)
        vecsexes = Vector{Char}(undef, playercount)

        pst.table = Table(namaes,
            splitkey((s)->Dan(parse(Int,s)), "dan", str, vecranks),
            splitkey((s)->parse(Float32,s), "rate", str, vecrates),
            splitkey((s)->(s[1]), "sx", str, vecsexes))
        pst.dced = Seat[]

        pst.hands = [Tiles(undef, 14) for i in 1:playercount]
        pst.melds = [Melds(undef, 12) for i in 1:playercount]
        pst.discard = [Tiles(undef, 32) for i in 1:playercount]
        pst.tedashi = [Tiles(undef, 32) for i in 1:playercount]
        pst.rchtile = [Tiles(undef, 1) for i in 1:playercount]

        pst.scores = Vector{Int32}(undef, playercount)
        pst.status = Vector{State}(undef, playercount)

        return matchset
    end

    function INIT(str::AbstractString, pst::PlayState)

        roundseed = splitkey((s)->(s), "seed", str,
                Vector{String}(undef, 6))
        number, pst.honba, pst.rbets, dice01, dice02 =
                map((s)-> s[1] - '0', roundseed)
        pst.doraid = [Wall[roundseed[end]]]

        pst.turn = 0
        pst.cycle = Round(number)
        pst.rolls = Dice(dice01), Dice(dice02)

        pst.dealer = Seat(parsekey((s)->parse(Int,s), "oya", str))
        pst.player = pst.dealer
        splitkey((s)->parse(Int32,s), "ten", str, pst.scores)

        playercount = length(pst.table.namaes)
        haikey = ("hai0", "hai1", "hai2", "hai3")

        for i in 1:playercount
            fill!(pst.hands[i], missing)
            fill!(pst.melds[i], missing)
            fill!(pst.discard[i], missing)
            fill!(pst.tedashi[i], missing)
            fill!(pst.rchtile[i], missing)
            splitkey(
                (s)->Wall[s], haikey[i],
                    str, pst.hands[i]
            )
        end

        fill!(pst.status, CLOSED)
        pst.result = nothing

        return roundinit
    end

    function AGA(str::AbstractString, pst::PlayState)

        fu, pt, lh = splitkey((s)->parse(Int,s), "ten", str,
                Vector{Int}(undef, 3))

        yaku = Tuple{Yaku,Int}[]

        if occursin("yakuman", str)
            ykm::Vector{Int} = splitkey((s)->parse(Int,s), "yakuman", str)
            for index in range(1, length(ykm); step = 1)
                push!(yaku, (Yaku(ykm[index]), 1))
            end
        else
            yku::Vector{Int} = splitkey((s)->parse(Int,s), "yaku", str)
            for index in range(1, length(yku); step = 2)
                push!(yaku, (Yaku(yku[index]), yku[index + 1]))
            end
        end

        han = mapreduce(x->x[2], +, yaku)

        dora = Tile[]; splitkey((s)->Wall[s], "doraHai", str, dora)
        ura = Tile[]; splitkey((s)->Wall[s], "doraHaiUra", str, ura)

        pst.result = RoundWin(pt, (han, fu), Limit(lh), yaku, dora, ura,
            Seat(parsekey((s)->parse(Int,s), "who", str)),
            Seat(parsekey((s)->parse(Int,s), "fromWho", str))
        )

        sc = splitkey((s)->parse(Int,s), "sc", str, Vector{Int}(undef, 8))
        for i in range(1, Int(length(sc) / 2); step = 1)
            pst.scores[i] = pst.scores[i] + sc[2*i]
        end

        return roundwin
    end

    function RYU(str::AbstractString, pst::PlayState)

        tierule::Ryuukyoku = tsuujou
        reveal = Seat[]

        if occursin("type", str)
            gettierule(s) = Ryuukyoku(findfirst(isequal(s),
                ("yao9", "reach4", "ron3", "kan4", "kaze4", "nm")
            ))
            tierule = parsekey(gettierule, "type", str)
        end

        occursin("hai0", str) && push!(reveal, Seat(0))
        occursin("hai1", str) && push!(reveal, Seat(1))
        occursin("hai2", str) && push!(reveal, Seat(2))
        occursin("hai3", str) && push!(reveal, Seat(3))

        pst.result = RoundTie(tierule, reveal)

        sc = splitkey((s)->parse(Int,s), "sc", str, Vector{Int}(undef, 8))
        for i in range(1, Int(length(sc) / 2); step = 1)
            pst.scores[i] = pst.scores[i] + sc[2*i]
        end

        return roundtie
    end

    function DRAW(str::AbstractString, pst::PlayState, who::Int)

        pst.turn = pst.turn + 1
        pst.player = Seat(who - 1)
        pst.hands[who][end] = Wall[str]

        return tiledraw
    end

    function DROP(str::AbstractString, pst::PlayState, who::Int)

        droppedtile = Wall[str]

        dropindex = findfirst(isequal(missing), pst.discard[who])
        pst.discard[who][dropindex] = droppedtile

        if !isequal(pst.hands[who][end], droppedtile)

            dropindex = findfirst(isequal(missing), pst.tedashi[who])
            pst.tedashi[who][dropindex] = droppedtile

            tileindex = findfirst(isequal(droppedtile), pst.hands[who])
            pst.hands[who][tileindex] = pst.hands[who][end]
        end

        pst.hands[who][end] = missing

        return tiledrop
    end

    function MELD(str::AbstractString, pst::PlayState)

        who = parsekey((s)->parse(Int,s), "who", str) + 1
        code = parsekey((s)->parse(Int,s), "m", str)
        from = mod(who + code & 0x3, 1:4)

        if      (code & 0x04 != 0)  meld = chi(code, who, from)
        elseif  (code & 0x08 != 0)  meld = pon(code, who, from)
        elseif  (code & 0x10 != 0)  meld = kkn(code, who, from)
        elseif  (code & 0x20 != 0)  meld = pei(code, who, from)
        else    #= code unknown =#  meld = kan(code, who, from)
        end

        meldindex = findfirst(isequal(missing), pst.melds[who])
        pst.player = Seat(who - 1)
        pst.melds[who][meldindex] = meld

        for tile in meld.with
            tileindex = findfirst(isequal(tile), pst.hands[who])
            pst.hands[who][tileindex] = missing
        end

        if !ismissing(pst.hands[who][end])
            tileindex = findfirst(isequal(meld.tile), pst.hands[who])
            pst.hands[who][tileindex] = pst.hands[who][end]
            pst.hands[who][end] = missing
        end

        if (meld.play != 暗槓)
            pst.status[who] = OPENED
        end

        return meldcall
    end

    function RCH(str::AbstractString, pst::PlayState)

        who = parsekey((s)->parse(Int,s), "who", str) + 1

        if parsekey((s)->parse(Int,s), "step", str) == 1
            pst.status[who] = RIICHI
        else
            pst.scores[who] -= 10
            tileindex = findfirst(isequal(missing), pst.discard[who])
            pst.rchtile[who][begin] = pst.discard[who][tileindex-1]
        end

        return riichicall
    end

    function DORA(str::AbstractString, pst::PlayState)

        push!(pst.doraid, parsekey((s)->Wall[s], "hai", str))

        return doraflip
    end

    function BYE(str::AbstractString, pst::PlayState)

        seat = Seat(parsekey((s)->parse(Int,s), "who", str))
        dcedindex = findfirst(isequal(seat), pst.dced)
        if isnothing(dcedindex) push!(pst.dced, seat)
        else deleteat!(pst.dced, dcedindex) end

        return playerdc
    end

    function SKIP(str::AbstractString, pst::PlayState)

        return noevents
    end

    function END(str::AbstractString, pst::PlayState)

        return matchend
    end


text:
function spKLIEP_refit!(x::SparseIterate, Ψx, Ψy, supp::Vector{Int64})
    w = ones(Float64, length(x)) * 1e10

    for k in supp
        w[k] = 0.
    end

    f = CDKLIEPLoss(Ψx, Ψy)
    g = ProxL1(1., w)
    coordinateDescent!(x, f, g)
end

spKLIEP_refit!(x::SparseIterate, Ψx, Ψy) = spKLIEP_refit!(x::SparseIterate, Ψx, Ψy, findall(!iszero, x))

function boot_KLIEP(Ψx, Ψy; bootSamples::Int64=300)
    θ = KLIEP(Ψx, Ψy, CD_KLIEP())
    θhat = convert(Vector, θ)
    θb = Matrix{Float64}(undef, length(θhat), bootSamples)

    m, nx = size(Ψx)
    ny = size(Ψy, 2)
    bΨx = similar(Ψx)
    bΨy = similar(Ψy)
    x_ind = Vector{Int64}(undef, nx)
    y_ind = Vector{Int64}(undef, ny)

    for b=1:bootSamples
       sample!(1:nx, x_ind)
       sample!(1:ny, y_ind)

       _fill_boot_Psi!(bΨx, Ψx, x_ind)
       _fill_boot_Psi!(bΨy, Ψy, y_ind)

       KLIEP!(θ, bΨx, bΨy)
       θb[:, b] .= θ
    end

    BootstrapEstimates(θhat, θb)
end

function boot_spKLIEP(Ψx, Ψy, θfs, Hinv; bootSamples::Int64=300)

    # generate bootstrap samples first since we want to use the same
    # samples for each coordinate j=1...m
    m, nx = size(Ψx)
    ny = size(Ψy, 2)

    x_ind = Matrix{Int16}(undef, nx, bootSamples)
    y_ind = Matrix{Int16}(undef, ny, bootSamples)

    for b=1:bootSamples
       sample!(1:nx, view(x_ind, :, b))
       sample!(1:ny, view(y_ind, :, b))
    end
    θhat = Vector{Float64}(undef, m)
    θb = Matrix{Float64}(undef, length(θhat), bootSamples)

    for j=1:m
        supp3 = findSupp3(θfs, Hinv[j], j)

        # obtain θhat_j
        bΨx = Ψx[supp3, :]
        bΨy = Ψy[supp3, :]
        bbΨx = similar(bΨx)
        bbΨy = similar(bΨy)
        θ = KLIEP(bΨx, bΨy, CD_KLIEP())
        θhat[j] = θ[end]

        # bootstrap
        for b=1:bootSamples
           _fill_boot_Psi!(bbΨx, bΨx, view(x_ind, :, b))
           _fill_boot_Psi!(bbΨy, bΨy, view(y_ind, :, b))

           KLIEP!(θ, bbΨx, bbΨy)
           θb[j, b] = θ[end]
        end
    end

    BootstrapEstimates(θhat, θb)
end

function boot_spKLIEPfull(Ψx, Ψy, θfs, Hinv, λ1, λ2; bootSamples::Int64=300)

    m, nx = size(Ψx)
    ny = size(Ψy, 2)

    θhat = Vector{Float64}(undef, m)
    θb = Matrix{Float64}(undef, length(θhat), bootSamples)

    # compute third stage estimator first
    for j=1:m
        supp3 = findSupp3(θfs, Hinv[j], j)

        # obtain θhat_j
        bΨx = Ψx[supp3, :]
        bΨy = Ψy[supp3, :]
        θ = KLIEP(bΨx, bΨy, CD_KLIEP())
        θhat[j] = θ[end]
    end

    x_ind = Vector{Int64}(undef, nx)
    y_ind = Vector{Int64}(undef, ny)
    Ψx_boot = similar(Ψx)
    Ψy_boot = similar(Ψy)

    for b=1:bootSamples
        sample!(1:nx, x_ind)
        sample!(1:ny, y_ind)

        _fill_boot_Psi!(Ψx_boot, Ψx, x_ind)
        _fill_boot_Psi!(Ψy_boot, Ψy, y_ind)

        # first step
        θfs_boot = deepcopy(θfs)
        spKLIEP!(θfs_boot, Ψx_boot, Ψy_boot, λ1, CD_KLIEP())
        spKLIEP_refit!(θfs_boot, Ψx_boot, Ψy_boot)

        # second and third step
        H_boot = KLIEP_Hessian(θfs_boot, Ψy_boot)
        for j=1:m
            ω = Hinv_row(H_boot, j, λ2)
            supp3_boot = findSupp3(θfs_boot, ω, j)

            bΨx_boot = Ψx_boot[supp3_boot, :]
            bΨy_boot = Ψy_boot[supp3_boot, :]

            θ = KLIEP(bΨx_boot, bΨy_boot, CD_KLIEP())
            θb[j, b] = θ[end]
        end
    end

    BootstrapEstimates(θhat, θb)
end

# S_delta --- support of Delta
function boot_oracleKLIEP(Ψx, Ψy, S_delta; bootSamples::Int64=300)

    # generate bootstrap samples first since we want to use the same
    # samples for each coordinate j=1...m
    m, nx = size(Ψx)
    ny = size(Ψy, 2)

    x_ind = Matrix{Int16}(undef, nx, bootSamples)
    y_ind = Matrix{Int16}(undef, ny, bootSamples)

    for b=1:bootSamples
       sample!(1:nx, view(x_ind, :, b))
       sample!(1:ny, view(y_ind, :, b))
    end
    θhat = Vector{Float64}(undef, m)
    θb = Matrix{Float64}(undef, length(θhat), bootSamples)

    for j=1:m
        S = copy(S_delta)
        pos_j = findfirst(isequal(j), S)
        if pos_j == nothing
            push!(S, j)
        else
            S[pos_j], S[end] = S[end], S[pos_j]
        end

        # obtain θhat_j
        bΨx = Ψx[S, :]
        bΨy = Ψy[S, :]
        bbΨx = similar(bΨx)
        bbΨy = similar(bΨy)
        θ = KLIEP(bΨx, bΨy, CD_KLIEP())
        θhat[j] = θ[end]

        # bootstrap
        for b=1:bootSamples
           _fill_boot_Psi!(bbΨx, bΨx, view(x_ind, :, b))
           _fill_boot_Psi!(bbΨy, bΨy, view(y_ind, :, b))

           KLIEP!(θ, bbΨx, bbΨy)
           θb[j, b] = θ[end]
        end
    end

    BootstrapEstimates(θhat, θb)
end

function boot_gaussKLIEP(Ψx, Ψy, θhat, Hinv; bootSamples::Int64=300)

    # generate bootstrap samples first since we want to use the same
    # samples for each coordinate j=1...m
    m, nx = size(Ψx)
    ny = size(Ψy, 2)

    # compute weights corresponding to rhat
    w = transpose(Ψy) * θhat
    w .= exp.(w)
    w ./= mean(w)

    # means
    mux = vec( mean(Ψx, dims = 2) )
    muy = vec( mean(Ψy, weights(w), 2) )

    θb = Matrix{Float64}(undef, length(θhat), bootSamples)

    # store Gaussian multipliers
    gx = Vector{Float64}(undef, nx)
    gy = Vector{Float64}(undef, ny)

    tmp1 = Vector{Float64}(undef, m)
    tmp2 = Vector{Float64}(undef, m)
    # bootstrap
    for b=1:bootSamples
        rand!(Normal(), gx)
        rand!(Normal(), gy)

        fill!(tmp1, 0.)
        fill!(tmp2, 0.)
        for i=1:nx
            @. tmp1 += (Ψx[:, i] - mux) * gx[i]
        end

        for i=1:ny
            @. tmp2 += (Ψy[:, i] * w[i] - muy) * gy[i]
        end
        @. tmp1 = tmp1 / nx + tmp2 / ny

        for j=1:m
            θb[j, b] = θhat[j] + dot(Hinv[j], tmp1)
        end
    end

    BootstrapEstimates(θhat, θb)
end


text:
module WeyandDB

include("./CSV.jl")
include("./FileOps.jl")
using .CSV, .FileOps

export eachpair

const DATEDICT = Dict{Vector{String}, Int}(
    ["ja", "jan", "january"]   => 01,
    ["fe", "feb", "february"]  => 02,
    ["ma", "mar", "march"]     => 03,
    ["ap", "apr", "april"]     => 04,
    ["may"]                    => 05,
    ["ju", "jun", "june"]      => 06,
    ["jl", "jul", "july"]      => 07,
    ["au", "aug", "august"]    => 08,
    ["se", "sep", "september"] => 09,
    ["oc", "oct", "october"]   => 10,
    ["no", "nov", "november"]  => 11,
    ["de", "dec", "december"]  => 12
)
# ============================================================================ #
struct DB
    db::Vector{Pair{Int,String}}
end
get_ids(db::DB) = first.(db.db)
Base.length(db::DB) = length(db.db)
Base.getindex(db::DB, k::Integer) = db.db[k]
function Base.getindex(db::DB; id::Integer=-1)
    k = findfirst(isequal(id), get_ids(db))
    k == nothing && error("No entry with id $(id) exists in the database")
    return db.db[k]
end
# ============================================================================ #
struct DBIterator
    db::DB
end

eachpair(db::DB) = DBIterator(db)

Base.IteratorSize() = HasLength()
Base.IteratorEltype() = HasEltype()
Base.length(db::DBIterator) = length(db.db)
Base.eltype(db::DBIterator) = Tuple{Vector{Float64}, Vector{Float64}}
Base.iterate(db::DBIterator) = Base.iterate(db, 1)
function Base.iterate(db::DBIterator, state::Integer)
    state > length(db) && return nothing
    return get_data(db.db, state), state+1
end
# ============================================================================ #
function get_month_id(str::AbstractString)
    for (k, v) in DATEDICT
        if str in k
            return v
        end
    end
    return -1
end
# ============================================================================ #
function get_id(file::AbstractString)
    idstr = replace(basename(file), "S_SPK.csv"=>"")
    m = match(r"(\d{4})([A-Za-z]{2,3})(\d+)", idstr)
    m == nothing && error("Failed to parse file name into id")

    yr = parse(Int, m[1])
    mo = get_month_id(lowercase(m[2]))
    id = parse(Int, m[3])

    return yr * Int(1e5) + mo * Int(1e3) + id
end
# ============================================================================ #
function get_database(::String, f::Function=x->true)
    files = find_files(joinpath(@__DIR__, "..", "data"), r".*\.csv")
    filter!(files) do file
        f(get_id(file))
    end
    ids = get_id.(files)
    ks = sortperm(ids)
    return DB(map((x,y)->x=>y, ids[ks], files[ks]))
end
# ============================================================================ #
function get_data(db::DB, idx::Integer=0; id::Integer=-1)
    MS2SEC = 0.001

    ((idx < 1 && id < 0) || idx > length(db)) && error("Invalid index / id!")
    if idx > 0
        file = db[idx].second
    elseif id > -1
        file = db[id=id].second
    end

    d = CSV.parse(file, ',', [Float64, Int])

    lgn = d["timestamps_ms"][findall(isequal(0), d["ids"])] .* MS2SEC
    ret = d["timestamps_ms"][findall(>(0), d["ids"])] .* MS2SEC

    # some records have negative timestamps, so shift both accordingly (plus
    # 100ms so we don't have any timestamps at exactly 0)
    t0 = min(minimum(ret), minimum(lgn))
    t0 = t0 < 0 ? abs(t0) + 0.1 : 0.1

    # retina recordings are S-potentials, so they do NOT include conduction
    # delays, thus we shift the retinal timestamps backwards by 2.4 ms (relative
    # to the LGN) to avoid any weirdness that the short delay may introduce
    # NOTE: 2.4 ms is the shifte required to match the median
    # S-potential -> LGN delay (0.4 ms) to the median RGC->LGN delay from
    # PairsDB (2.8 ms)
    return ret .+ (t0 - 0.0024), lgn .+ t0
end
# ============================================================================ #
end


text:
function addseedingdateoftransplant!(df)
    df[!,:seeding_date_of_transplant] = Array{Union{Missing, Date},1}(missing, nrow(df))
    for row in eachrow(df)
        if !ismissing(row[:planting])
            if (row[:planting] == "T") & !ismissing(row[:number_of_days_in_cell] & !ismissing(row[:date_in]))
                row[:seeding_date_of_transplant] = row[:date_in] - Day(row[:number_of_days_in_cell])
            end
        end
    end
end


text:
#! /usr/bin/julia

# Rosetta Code, Sorting algorithms/Selection sort

function properdivisors{T<:Integer}(n::T)
    0 < n || throw(ArgumentError("number to be factored must be ≥ 0, got $n"))
    1 < n || return T[]
    !isprime(n) || return T[one(T), n]
    f = factor(n)
    d = T[one(T)]
    for (k, v) in f
        c = T[k^i for i in 0:v]
        d = d*c'
        d = reshape(d, length(d))
    end
    sort!(d)
    return d[1:end-1]
end

lo = 1
hi = 10
println("List the proper divisors for ", lo, " through ", hi, ".")
for i in lo:hi
    println(@sprintf("%4d", i), " ", properdivisors(i))
end

hi = 2*10^4
println("\nFind the numbers within [", lo, ",", hi, "] having the most divisors.")

maxdiv = 0
nlst = Int[]

for i in lo:hi
    ndiv = length(properdivisors(i))
    if ndiv > maxdiv
        maxdiv = ndiv
        nlst = [i]
    elseif ndiv == maxdiv
        push!(nlst, i)
    end
end

println(nlst, " have the maximum proper divisor count of ", maxdiv, ".")

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using BinDeps
if !ispath(abspath(first(DEPOT_PATH), "config","DFControl"))
    mkpath(abspath(first(DEPOT_PATH), "config","DFControl"))
end
if !ispath(abspath(first(DEPOT_PATH), "config","DFControl", "user_defaults.jl"))
    open(abspath(first(DEPOT_PATH), "config","DFControl", "user_defaults.jl"), "w") do f
        write(f,"#Generated by DFControl.")
    end
end

relpath = x -> joinpath(@__DIR__, x)

pythonpath = relpath("python2")
prevpythfound = false
package_basepath = joinpath(@__DIR__, "../../")
for d in readdir(package_basepath)
    test_pythonpath = joinpath(package_basepath, d, "deps/python2/")
    if ispath(test_pythonpath) && !ispath(pythonpath) && !isempty(readdir(test_pythonpath))
        @info "Previous python2 found at $test_pythonpath."
        cp(test_pythonpath, pythonpath)
        prevpythonfound = true
        break
    end
end
if !ispath(pythonpath)
    mkdir(pythonpath)
    dlpath = relpath("downloads")
    if !ispath(dlpath)
        mkdir(dlpath)
    end
    #From Conda installation
    url = "https://repo.continuum.io/miniconda/Miniconda2-latest-"
    if Sys.isapple()
        url *= "MacOSX"
    elseif Sys.islinux()
        url *= "Linux"
    elseif Sys.iswindows()
        url = "https://repo.continuum.io/miniconda/Miniconda2-4.5.4-"
        url *= "Windows"
    else
        error("Unsuported OS.")
    end
    url *= Sys.WORD_SIZE == 64 ? "-x86_64" : "-x86"
    url *= Sys.iswindows() ? ".exe" : ".sh"

    if Sys.isunix()
        installer = joinpath(dlpath, "installer.sh")
    end
    if Sys.iswindows()
        installer = joinpath(dlpath, "installer.exe")
    end
    download(url, installer)
    if Sys.isunix()
        chmod(installer, 33261)  # 33261 corresponds to 755 mode of the 'chmod' program
        run(`$installer -b -f -p $pythonpath`)
    end
    if Sys.iswindows()
        run(Cmd(`$installer /S /AddToPath=0 /RegisterPython=0 /D=$pythonpath`, windows_verbatim=true))
    end

    tarpath = relpath("cif2cell.tar.gz")
    download("https://sourceforge.net/projects/cif2cell/files/latest/download", tarpath)
    run(unpack_cmd("cif2cell.tar.gz", @__DIR__, ".gz",".tar"))
    cif2celldir = relpath("cif2cell-1.2.10")
    cd(cif2celldir)
    pyex = Sys.iswindows() ? joinpath(pythonpath, "python") : joinpath(pythonpath, "bin", "python2")
    run(`$pyex setup.py install --prefix=$pythonpath`)
    cd("..")
#stupid urlopen
    starfile = Sys.iswindows() ? joinpath(pythonpath, "lib", "site-packages", "StarFile.py") : joinpath(pythonpath, "lib", "python2.7", "site-packages", "StarFile.py")
    starsource = read(starfile, String)
    starsource = replace(starsource, "filestream = urlopen(filename)"=>"filestream = urlopen('file:' + filename)")
    write(starfile, starsource)
#
    rm(tarpath)
    rm(relpath("cif2cell-1.2.10"), recursive=true)
    rm(relpath("downloads"), recursive=true)
end

# spglibpath = relpath("spglib/lib/libsymspg.so")
# if !ispath(spglibpath)
#     spglibdir = relpath("spglib")
#     if !ispath(spglibdir)
#         mkdir(spglibdir)
#     end

#     spgdownloads = joinpath(spglibdir, "downloads")
#     if !ispath(spgdownloads)
#         mkdir(spgdownloads)
#     end
#     spgzip = joinpath(spgdownloads, "spglib-master.zip")
#     download("https://github.com/atztogo/spglib/archive/master.zip", spgzip)
#     run(unpack_cmd(spgzip, spglibdir, ".zip",""))
#     cd(joinpath(spglibdir, "spglib-master"))
#     run(`mkdir build`)
#     cd("build")
#     run(`cmake -DCMAKE_INSTALL_PREFIX="" ..`)
#     run(`make`)
#     run(`make DESTDIR=$spglibdir install`)
# end

# cd(@__DIR__)
# open("deps.jl", "w") do f
#     write(f, """const SPGLIB = joinpath(@__DIR__, "spglib","lib","libsymspg.so")\n""")
# end

include("asset_init.jl")


text:
using Plots
using ImplicitEquations

f(x, y) = y^2 - x^2 - x^3
g(x, y) = y^2 - x^3
h(x, y) = y^3 - x^2 * y + x^4

plot()
plot!(Eq(f, 0), color=:red)
plot!(Eq(g, 0))
plot!(Eq(h, 0))

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
##
using Roots
function GE(π, γ, ω1A, ω1B, ω2A, ω2B)
    f(p) = π * (ω1A + p*ω1B) + ((ω2A + p*ω2B) / (1 + p * ((1 - π)/(p * π))^(1/γ))) - ω1A - ω2A
    p = round(find_zero(f, 1.5), digits= 3)
    C1A = π * (ω1A + p*ω1B)
    C1B = ((1 - π) / p) * (ω1A + p*ω1B)
    C2A = (ω2A + p*ω2B) / (1 + p * ((1 - π)/(p * π))^(1/γ))
    C2B = ((ω2A + p*ω2B) / (1 + p * ((1 - π)/(p * π))^(1/γ))) * ((1 - π)/(p * π))^(1/γ)
    print("p = $p \n", "C1A = $C1A \n", "C1B = $C1B \n", "C2A = $C2A \n", "C2B = $C2B \n")
end
##
#4(e) No aggregate uncertanity
GE(0.5, 2, 2, 5, 5, 2)
#4(f) There is no aggregate uncertainty but for person A the good world is more likely to happen
GE(0.25, 2, 2, 5, 5, 2)
#4(f) There is no aggregate uncertainty but for person B the good world is more likely to happen
GE(0.75, 2, 2, 5, 5, 2)
#4(g)(e) γ=1.5
GE(0.5, 1.5, 2, 5, 5, 2)
#4(g)(e) γ=3
GE(0.5, 3, 2, 5, 5, 2)
#4(g)(f) γ=1.5 π=0.25
GE(0.25, 1.5, 2, 5, 5, 2)
#4(g)(f) γ=1.5 π=0.75
GE(0.75, 1.5, 2, 5, 5, 2)
#4(g)(f) γ=3 π=0.25
GE(0.25, 3, 2, 5, 5, 2)
#4(g)(f) γ=3 π=0.75
GE(0.75, 3, 2, 5, 5, 2)
#4(g)(e) γ=2 π=0.5
GE(0.5, 2, 2, 3, 5, 2)
#4(g)(e) γ=1.5 π=0.5
GE(0.5, 1.5, 2, 3, 5, 2)
#4(g)(e) γ=3 π=0.5
GE(0.5, 3, 2, 3, 5, 2)

max_stars_repo_path: GE solver.jl

max_issues_repo_path: GE solver.jl

max_forks_repo_path: GE solver.jl