text:
function check_data_sizes(data::AbstractMatrix, raxis, (l,u), (λ,μ))
    if blocksize(data,1) ≠ l + u + 1 && !(blocksize(data,1) == 0 && (-l > u || -λ > μ))
        throw(ArgumentError("Data matrix must have number of row blocks equal to number of block bands"))
    end
    for K = blockaxes(data,1)
        if length(axes(data,1)[K]) ≠ λ + μ + 1 && !(-λ > μ)
            throw(ArgumentError("Data matrix must have row block sizes equal to number of subblock bands"))
        end
    end
end

function _BandedBlockBandedMatrix end

# Represents a block banded matrix with banded blocks
#   similar to BandedMatrix{BandedMatrix{T,Matrix{T}}}
# Here the data is stored by blocks, in a way that is consistent with
# BandedMatrix
#

struct BandedBlockBandedMatrix{T, BLOCKS, RAXIS<:AbstractUnitRange{Int}} <: AbstractBlockBandedMatrix{T}
    data::BLOCKS
    raxis::RAXIS

    l::Int  # block lower bandwidth
    u::Int  # block upper bandwidth
    λ::Int  # sub lower bandwidth
    μ::Int  # sub upper bandwidth

    global function _BandedBlockBandedMatrix(data::AbstractMatrix,
                                             raxis::AbstractUnitRange{Int},
                                             lu::NTuple{2,Int}, λμ::NTuple{2,Int})
      @boundscheck check_data_sizes(data, raxis, lu, λμ)
      new{eltype(data), typeof(data), typeof(raxis)}(data, raxis, lu..., λμ...)
    end
end

const DefaultBandedBlockBandedMatrix{T} = BandedBlockBandedMatrix{T, PseudoBlockMatrix{T, Matrix{T}, NTuple{2,DefaultBlockAxis}}, DefaultBlockAxis}

@inline _BandedBlockBandedMatrix(data::AbstractMatrix, axes::NTuple{2,AbstractUnitRange{Int}}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) =
    _BandedBlockBandedMatrix(PseudoBlockArray(data,(blockedrange(Fill(sum(λμ)+1,sum(lu)+1)),axes[2])), axes[1], lu, λμ)

_BandedBlockBandedMatrix(data::AbstractMatrix,rblocksizes::AbstractVector{Int}, cblocksizes::AbstractVector{Int}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) =
    _BandedBlockBandedMatrix(data, (blockedrange(rblocksizes),blockedrange(cblocksizes)), lu, λμ)

_blocklengths2blocklasts(b::Fill) = iszero(FillArrays.getindex_value(b)) ? (1:1:0) : cumsum(b)
_bbb_data_axes(caxes, lu, λμ) = (blockedrange(Fill(max(0,sum(λμ)+1),max(0,sum(lu)+1))),caxes)

BandedBlockBandedMatrix{T,B,R}(::UndefInitializer, axes::NTuple{2,AbstractUnitRange{Int}}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B,R<:AbstractUnitRange{Int}} =
    _BandedBlockBandedMatrix(B(undef, _bbb_data_axes(axes[2],lu,λμ)), axes[1], lu, λμ)

BandedBlockBandedMatrix{T,B,R}(::UndefInitializer, rblocksizes::AbstractVector{Int}, cblocksizes::AbstractVector{Int}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B,R<:AbstractUnitRange{Int}} =
    BandedBlockBandedMatrix{T,B,R}(undef, (blockedrange(rblocksizes),blockedrange(cblocksizes)), lu, λμ)

BandedBlockBandedMatrix{T,B}(::UndefInitializer, axes::NTuple{2,AbstractUnitRange{Int}}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B} =
    BandedBlockBandedMatrix{T,B,typeof(axes[1])}(undef,axes,lu,λμ)

BandedBlockBandedMatrix{T,B}(::UndefInitializer, rblocksizes::AbstractVector{Int}, cblocksizes::AbstractVector{Int}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B} =
    BandedBlockBandedMatrix{T,B}(undef, (blockedrange(rblocksizes),blockedrange(cblocksizes)), lu, λμ)

"""
    BandedBlockBandedMatrix{T}(M::Union{UndefInitializer,UniformScaling,AbstractMatrix},
                               rows, cols, (l, u), (λ, μ))

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

# Examples

```jldoctest
julia> BandedBlockBandedMatrix(I, [3,4,3], [3,4,3], (1,1), (1,1))
3×3-blocked 10×10 BandedBlockBandedMatrix{Bool,BlockArrays.PseudoBlockArray{Bool,2,Array{Bool,2},Tuple{BlockArrays.BlockedUnitRange{Array{Int64,1}},BlockArrays.BlockedUnitRange{Array{Int64,1}}}},BlockArrays.BlockedUnitRange{Array{Int64,1}}}:
 1  0  ⋅  │  0  0  ⋅  ⋅  │  ⋅  ⋅  ⋅
 0  1  0  │  0  0  0  ⋅  │  ⋅  ⋅  ⋅
 ⋅  0  1  │  ⋅  0  0  0  │  ⋅  ⋅  ⋅
 ─────────┼──────────────┼─────────
 0  0  ⋅  │  1  0  ⋅  ⋅  │  0  0  ⋅
 0  0  0  │  0  1  0  ⋅  │  0  0  0
 ⋅  0  0  │  ⋅  0  1  0  │  ⋅  0  0
 ⋅  ⋅  0  │  ⋅  ⋅  0  1  │  ⋅  ⋅  0
 ─────────┼──────────────┼─────────
 ⋅  ⋅  ⋅  │  0  0  ⋅  ⋅  │  1  0  ⋅
 ⋅  ⋅  ⋅  │  0  0  0  ⋅  │  0  1  0
 ⋅  ⋅  ⋅  │  ⋅  0  0  0  │  ⋅  0  1

julia> BandedBlockBandedMatrix(Ones{Int}(10,13), [3,4,3], [4,5,4], (1,1), (1,1))
3×3-blocked 10×13 BandedBlockBandedMatrix{Int64,BlockArrays.PseudoBlockArray{Int64,2,Array{Int64,2},Tuple{BlockArrays.BlockedUnitRange{Array{Int64,1}},BlockArrays.BlockedUnitRange{Array{Int64,1}}}},BlockArrays.BlockedUnitRange{Array{Int64,1}}}:
 1  1  ⋅  ⋅  │  1  1  ⋅  ⋅  ⋅  │  ⋅  ⋅  ⋅  ⋅
 1  1  1  ⋅  │  1  1  1  ⋅  ⋅  │  ⋅  ⋅  ⋅  ⋅
 ⋅  1  1  1  │  ⋅  1  1  1  ⋅  │  ⋅  ⋅  ⋅  ⋅
 ────────────┼─────────────────┼────────────
 1  1  ⋅  ⋅  │  1  1  ⋅  ⋅  ⋅  │  1  1  ⋅  ⋅
 1  1  1  ⋅  │  1  1  1  ⋅  ⋅  │  1  1  1  ⋅
 ⋅  1  1  1  │  ⋅  1  1  1  ⋅  │  ⋅  1  1  1
 ⋅  ⋅  1  1  │  ⋅  ⋅  1  1  1  │  ⋅  ⋅  1  1
 ────────────┼─────────────────┼────────────
 ⋅  ⋅  ⋅  ⋅  │  1  1  ⋅  ⋅  ⋅  │  1  1  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  │  1  1  1  ⋅  ⋅  │  1  1  1  ⋅
 ⋅  ⋅  ⋅  ⋅  │  ⋅  1  1  1  ⋅  │  ⋅  1  1  1
```
"""
BandedBlockBandedMatrix

BandedBlockBandedMatrix{T}(::UndefInitializer, args...) where T =
    DefaultBandedBlockBandedMatrix{T}(undef, args...)

# Auxiliary outer constructors

function convert(::Type{<:BandedBlockBandedMatrix}, B::BandedMatrix)
    if isdiag(B) # TODO really should be one block
        _BandedBlockBandedMatrix(copy(B.data), Fill(1,size(B,1)), Fill(1,size(B,2)), (0,0), (0,0))
    else
        _BandedBlockBandedMatrix(copy(B.data), [size(B,1)], [size(B,2)], (0,0), (B.l,B.u))
    end
end

convert(::Type{BandedBlockBandedMatrix{T,BLOCKS,RAXIS}}, A::BandedBlockBandedMatrix) where {T,BLOCKS,RAXIS} =
    _BandedBlockBandedMatrix(convert(BLOCKS, A.data), convert(RAXIS, A.raxis), (A.l, A.u), (A.λ, A.μ))

function BandedBlockBandedMatrix{T,B,R}(Z::Zeros, axes::NTuple{2,AbstractUnitRange{Int}},
                                       lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B,R<:AbstractUnitRange{Int}}
   if size(Z) ≠ map(length,axes)
       throw(DimensionMismatch())
   end
   blocks = fill!(B(undef, _bbb_data_axes(axes[2],lu,λμ)), zero(T))
   _BandedBlockBandedMatrix(blocks, convert(R, axes[1]), lu, λμ)
end

function BandedBlockBandedMatrix{T,B,R}(E::Eye, axes::NTuple{2,AbstractUnitRange{Int}},
                                       lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B,R<:AbstractUnitRange{Int}}
    if size(E) ≠ map(length,axes)
        throw(DimensionMismatch())
    end
    ret = BandedBlockBandedMatrix{T,B,R}(Zeros{T}(E), axes, lu, λμ)
    ret[diagind(ret)] .= one(T)
    ret
end

function BandedBlockBandedMatrix{T,B,R}(A::UniformScaling, axes::NTuple{2,AbstractUnitRange{Int}},
                                       lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B,R<:AbstractUnitRange{Int}}
    ret = BandedBlockBandedMatrix{T,B,R}(Zeros{T}(map(length,axes)), axes, lu, λμ)
    ret[diagind(ret)] .= convert(T, A.λ)
    ret
end

BandedBlockBandedMatrix{T,B}(m::Union{AbstractMatrix, UniformScaling},
                            axes::NTuple{2,AbstractUnitRange{Int}},
                           lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B} =
  BandedBlockBandedMatrix{T,B,typeof(axes[1])}(m, axes, lu, λμ)

BandedBlockBandedMatrix{T}(m::Union{AbstractMatrix, UniformScaling},
                            axes::NTuple{2,AbstractUnitRange{Int}},
                           lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where T =
  DefaultBandedBlockBandedMatrix{T}(m, axes, lu, λμ)

BandedBlockBandedMatrix{T,B,R}(m::Union{AbstractMatrix, UniformScaling},
                            rdims::AbstractVector{Int}, cdims::AbstractVector{Int},
                           lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B,R<:AbstractUnitRange{Int}} =
  BandedBlockBandedMatrix{T,B,R}(m, (blockedrange(rdims),blockedrange(cdims)), lu, λμ)

BandedBlockBandedMatrix{T,B}(m::Union{AbstractMatrix, UniformScaling},
                            rdims::AbstractVector{Int}, cdims::AbstractVector{Int},
                           lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B} =
  BandedBlockBandedMatrix{T,B}(m, (blockedrange(rdims),blockedrange(cdims)), lu, λμ)

BandedBlockBandedMatrix{T}(m::Union{AbstractMatrix, UniformScaling},
                            rdims::AbstractVector{Int}, cdims::AbstractVector{Int},
                           lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where T =
  BandedBlockBandedMatrix{T}(m, (blockedrange(rdims),blockedrange(cdims)), lu, λμ)

BandedBlockBandedMatrix(A::Union{AbstractMatrix,UniformScaling},
                        axes::NTuple{2, AbstractUnitRange{Int}}, lu::NTuple{2,Int},
                        λμ::NTuple{2,Int}) =
    BandedBlockBandedMatrix{eltype(A)}(A, axes, lu, λμ)
BandedBlockBandedMatrix(A::Union{AbstractMatrix,UniformScaling},
                        rdims::AbstractVector{Int}, cdims::AbstractVector{Int},
                        lu::NTuple{2,Int}, λμ::NTuple{2,Int}) =
    BandedBlockBandedMatrix(A, (blockedrange(rdims),blockedrange(cdims)), lu, λμ)

function BandedBlockBandedMatrix{T,Blocks,RR}(A::AbstractMatrix, axes::NTuple{2,AbstractUnitRange{Int}}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,Blocks,RR<:AbstractUnitRange{Int}}
    ret = BandedBlockBandedMatrix{T,Blocks,RR}(Zeros{T}(size(A)), axes, lu, λμ)
    L,M = λμ
    for J = blockaxes(ret,2), K = blockcolsupport(ret, J)
        kr, jr = axes[1][K], axes[2][J]

        # We have the correct block - now we need to only add the entries from
        # the correct bands
        B = view(A, kr, jr)
        R = view(ret, K, J)
        rows, cols = size(B)
        for λ = 1:L, j = 1:min(rows-λ, cols)
            view(R, j+λ, j) .= view(B, j+λ, j)
        end
        for i = 1:min(rows, cols)
            view(R, i, i) .= view(B, i, i)
        end
        for μ = 1:M, k = 1:min(rows, cols-μ)
            view(R, k, k+μ) .= view(B, k, k+μ)
        end
    end
    ret
end

BandedBlockBandedMatrix{T}(A::AbstractMatrix) where T =
    copyto!(BandedBlockBandedMatrix{T}(undef, axes(A), blockbandwidths(A), subblockbandwidths(A)), A)

BandedBlockBandedMatrix(A::AbstractMatrix{T}) where T = BandedBlockBandedMatrix{T}(A)

BandedBlockBandedMatrix(A::AbstractMatrix, rdims::AbstractVector{Int}, cdims::AbstractVector{Int}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) =
    BandedBlockBandedMatrix(A, (blockedrange(rdims), blockedrange(cdims)), lu, λμ)

similar(A::BandedBlockBandedMatrix, ::Type{T}, axes::NTuple{2,AbstractUnitRange{Int}}) where T =
    BandedBlockBandedMatrix{T}(undef, axes, blockbandwidths(A), subblockbandwidths(A))

@inline similar(A::BandedBlockBandedMatrix, ::Type{T}, axes::Tuple{BlockedUnitRange,AbstractUnitRange{Int}}) where T =
    BandedBlockBandedMatrix{T}(undef, axes, blockbandwidths(A), subblockbandwidths(A))
@inline similar(A::BandedBlockBandedMatrix, ::Type{T}, axes::Tuple{AbstractUnitRange{Int},BlockedUnitRange}) where T =
    BandedBlockBandedMatrix{T}(undef, axes, blockbandwidths(A), subblockbandwidths(A))
@inline similar(A::BandedBlockBandedMatrix, ::Type{T}, axes::Tuple{BlockedUnitRange,BlockedUnitRange}) where T =
    BandedBlockBandedMatrix{T}(undef, axes, blockbandwidths(A), subblockbandwidths(A))

similar(A::BandedBlockBandedMatrix{T}, axes::NTuple{2,AbstractUnitRange{Int}}) where T =
    similar(Matrix{T}, map(length,axes)...)

axes(A::BandedBlockBandedMatrix) = (A.raxis, axes(A.data,2))

function sparse(A::BandedBlockBandedMatrix)
    i = Vector{Int}()
    j = Vector{Int}()
    z = Vector{eltype(A)}()
    for J = blockaxes(A,2), K = blockcolsupport(A, J)
        B = view(A, K, J)
        ĩ = _banded_rowval(B)
        j̃ = _banded_colval(B)
        z̃ = _banded_nzval(B)
        ĩ .+= first(axes(A,1)[K])-1
        j̃ .+= first(axes(A,2)[J])-1
        append!(i, ĩ)
        append!(j, j̃)
        append!(z, z̃)
    end
    sparse(i, j, z, size(A)...)
end

################################
# BandedBlockBandedMatrix Interface #
################################

bandedblockbandedcolumns(L::AbstractColumnMajor) = BandedBlockBandedColumnMajor()
bandedblockbandedcolumns(_) = BandedBlockBandedColumns{UnknownLayout}()

MemoryLayout(::Type{<:BandedBlockBandedMatrix{<:Any,BLOCKS}}) where BLOCKS =
    bandedblockbandedcolumns(MemoryLayout(BLOCKS))
bandedblockbandedbroadcaststyle(_) = BandedBlockBandedStyle()
BroadcastStyle(::Type{<:BandedBlockBandedMatrix{<:Any,BLOCKS}}) where BLOCKS =
    bandedblockbandedbroadcaststyle(BroadcastStyle(BLOCKS))

isbandedblockbanded(_) = false
isbandedblockbanded(::BandedBlockBandedMatrix) = true

blockbandwidths(A::BandedBlockBandedMatrix) = (A.l, A.u)

"""
    subblockbandwidths(A)

returns the sub-block bandwidths of `A`, where `A` is a banded-block-banded
matrix. In other words, `A[Block(K,J)]` will return a `BandedMatrix` with
bandwidths given by `subblockbandwidths(A)`.
"""
subblockbandwidths(A::BandedBlockBandedMatrix) = (A.λ, A.μ)

# default is to use bandwidth
subblockbandwidths(A::AbstractMatrix) = bandwidths(A)

"""
    subblockbandwidth(A, i)

returns the sub-block lower (`i == 1`) or upper (`i == 2`) bandwidth of `A`,
where `A` is a banded-block-banded matrix. In other words, `A[Block(K,J)]` will
return a `BandedMatrix` with the returned lower/upper bandwidth.
"""
subblockbandwidth(A::AbstractArray, k::Integer) = subblockbandwidths(A)[k]

isdiag(A::BandedBlockBandedMatrix) = A.λ == A.μ == A.l == A.u

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

zeroblock(A::BandedBlockBandedMatrix, K::Int, J::Int) =
    BandedMatrix(Zeros{eltype(A)}(length.(getindex.(axes(A),(Block(K),Block(J))))), (A.λ, A.μ))

@inline function getblock(A::BandedBlockBandedMatrix, K::Int, J::Int)
    @boundscheck blockcheckbounds(A, K, J)
    if -A.l ≤ J - K ≤ A.u
        convert(BandedMatrix, view(A, Block(K, J)))
    else
        zeroblock(A, K, J)
    end
end

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

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

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

Base.size(arr::BandedBlockBandedMatrix) =
    @inbounds return map(length,axes(arr))

@inline function getindex(A::BandedBlockBandedMatrix{T}, i::Int, j::Int) where T
    @boundscheck checkbounds(A, i, j)
    bi = findblockindex.(axes(A), (i,j))
    @inbounds v = view(A, block.(bi)...)[blockindex.(bi)...]
    return v
end

@inline function setindex!(A::BandedBlockBandedMatrix{T}, v, i::Int, j::Int) where T
    @boundscheck checkbounds(A, i, j)
    BI,BJ = findblockindex.(axes(A), (i,j))
    if -A.l ≤ Int(block(BJ)-block(BI)) ≤ A.u
        V = view(A, block(BI),block(BJ))
         @inbounds V[blockindex(BI),blockindex(BJ)] = convert(T, v)::T
    elseif !iszero(v)
        throw(BandError(A))
    end
    return v
end

## structured matrix methods ##
function layout_replace_in_print_matrix(::AbstractBandedBlockBandedLayout, A, i, j, s)
    bi = findblockindex.(axes(A), (i,j))
    I,J = block.(bi)
    i,j = blockindex.(bi)
    l,u = blockbandwidths(A)
    λ,μ = subblockbandwidths(A)
    -l ≤ Int(J-I) ≤ u && -λ ≤ j-i ≤ μ ? s : Base.replace_with_centered_mark(s)
end

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

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

    @boundscheck (bandwidth(v, 1) > A.λ || bandwidth(v, 2) > A.μ) && throw(BandError())
    V = view(A, Block(K), Block(J))
    V .= v
    return A
end
#
# @propagate_inbounds function Base.setindex!(block_array::BlockArray{T, N}, v, block_index::BlockIndex{N}) where {T,N}
#     getblock(block_array, block_index.I...)[block_index.α...] = v
# end

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

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

##################
# BandedBlockBandedBlock
#
#   views of the blocks satisfy the BandedMatrix interface, and are memory-compatible
#   with BLASBandedMatrix.
##################

const SubBandedBlockBandedMatrix{T,R1,R2} =
    SubArray{T,2,<:BandedBlockBandedMatrix{T},<:Tuple{<:BlockSlice{R1},<:BlockSlice{R2}}}

sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{Block1}}}) = BandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) = BandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) = BandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) = BandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockRange1}}}) = BandedBlockBandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{BlockRange1}}}) = BandedBlockBandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{Block1}}}) = BandedBlockBandedColumnMajor()
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockIndexRange1}}}) = BandedBlockBandedColumnMajor()

sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{Block1}}}) = BandedLayout()
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) = BandedLayout()
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) = BandedLayout()
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) = BandedLayout()
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockRange1}}}) = BandedBlockBandedLayout()
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{BlockRange1}}}) = BandedBlockBandedLayout()
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{Block1}}}) = BandedBlockBandedLayout()
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockIndexRange1}}}) = BandedBlockBandedLayout()

sub_materialize(::AbstractBandedBlockBandedLayout, V, _) = BandedBlockBandedMatrix(V)

isbanded(A::SubArray{<:Any,2,<:BandedBlockBandedMatrix}) = MemoryLayout(typeof(A)) == BandedColumnMajor()
isbandedblockbanded(A::SubArray{<:Any,2,<:BandedBlockBandedMatrix}) = MemoryLayout(typeof(A)) == BandedBlockBandedColumnMajor()

subblockbandwidths(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{BlockRange1},<:BlockSlice{BlockRange1}}}) =
    subblockbandwidths(parent(V))

function blockbandwidths(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{BlockRange1},<:BlockSlice{Block1}}})
    A = parent(V)
    KR = parentindices(V)[1].block.indices[1]
    J = parentindices(V)[2].block
    shift = Int(KR[1])-Int(J)
    blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
end

function blockbandwidths(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{<:Block1},<:BlockSlice{<:BlockRange1}}})
    A = parent(V)
    K = parentindices(V)[1].block
    JR = parentindices(V)[2].block.indices[1]
    shift = Int(K)-Int(JR[1])
    l,u = blockbandwidths(A)
    l - shift, u + shift
end

function blockbandwidths(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{<:BlockRange1},<:BlockSlice{<:BlockRange1}}})
    A = parent(V)

    KR = parentindices(V)[1].block.indices[1]
    JR = parentindices(V)[2].block.indices[1]
    shift = Int(first(KR))-Int(first(JR))

    blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
end

const BandedBlockBandedBlock{T, BLOCKS, RAXIS} = SubArray{T,2,BandedBlockBandedMatrix{T, BLOCKS, RAXIS},<:Tuple{<:BlockSlice{Block1},<:BlockSlice{Block1}},false}

BroadcastStyle(::Type{<: BandedBlockBandedBlock}) = BandedStyle()

function inblockbands(V::SubArray{<:Any,2,<:AbstractMatrix,<:Tuple{<:BlockSlice{Block1},<:BlockSlice{Block1}},false})
    A = parent(V)
    K_sl, J_sl = parentindices(V)
    K, J = K_sl.block, J_sl.block
    l,u = blockbandwidths(A)
    -l ≤ Int(J-K) ≤ u
end

function parentblock(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) where T
    A = parent(V)
    K_sl, J_sl = parentindices(V)
    view(A, K_sl.block.block, J_sl.block.block)
end

function parentblock(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) where T
    A = parent(V)
    K_sl, J_sl = parentindices(V)
    view(A, K_sl.block, J_sl.block.block)
end

function parentblock(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) where T
    A = parent(V)
    K_sl, J_sl = parentindices(V)
    view(A, K_sl.block.block, J_sl.block)
end

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

######################################
# BandedMatrix interface  for Blocks #
######################################
@inline function bandwidths(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{Block1},BlockSlice{Block1}}}) where T
    inblockbands(V) && return subblockbandwidths(parent(V))
    (-720,-720)
end

function bandwidths(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) where T
   B = parentblock(V)
   K_sl, J_sl = parentindices(V)
   bandwidths(B) .+ (-1,1) .* bandshift(K_sl.block.indices[1],J_sl.block.indices[1])
end

function bandwidths(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) where T
    B = parentblock(V)
    K_sl, J_sl = parentindices(V)
    bandwidths(B) .+ (-1,1) .* bandshift(Base.OneTo(1),J_sl.block.indices[1])
 end

function bandwidths(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) where T
    B = parentblock(V)
    K_sl, J_sl = parentindices(V)
    bandwidths(B) .+ (-1,1) .* bandshift(K_sl.block.indices[1],Base.OneTo(1))
 end

function bandeddata(V::BandedBlockBandedBlock{T}) where T
    inblockbands(V) || return Array{T}(undef, 0, size(V,2))
    A = parent(V)
    u = A.u
    K_sl, J_sl = parentindices(V)
    K, J = K_sl.block, J_sl.block
    view(A.data, u + K - J + 1, J)
end

function bandeddata(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) where T
    A = parent(V)
    K_sl, J_sl = parentindices(V)
    view(bandeddata(parentblock(V)), :, J_sl.block.indices[1])
end

function bandeddata(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) where T
    A = parent(V)
    K_sl, J_sl = parentindices(V)
    view(bandeddata(parentblock(V)), :, J_sl.block.indices[1])
end

function bandeddata(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) where T
    A = parent(V)
    K_sl, J_sl = parentindices(V)
    bandeddata(parentblock(V))
end

@inline function inbands_getindex(V::BandedBlockBandedBlock, k::Int, j::Int)
    A = parent(V)
    banded_getindex(bandeddata(V), A.λ, A.μ, k, j)
end

@inline function inbands_setindex!(V::BandedBlockBandedBlock, v, k::Int, j::Int)
    A = parent(V)
    banded_setindex!(bandeddata(V), A.λ, A.μ, v, k, j)
end

# we need getindex as we use block-wise get index above
@propagate_inbounds function getindex(V::BandedBlockBandedBlock, k::Int, j::Int)
    @boundscheck checkbounds(V, k, j)
    A = parent(V)
    K,J = parentblocks2Int(V)
    if -A.l ≤ J-K ≤ A.u && -A.λ ≤ j-k ≤ A.μ
        inbands_getindex(V, k, j)
    else
        zero(eltype(V))
    end
end

@propagate_inbounds function setindex!(V::BandedBlockBandedBlock, v, k::Int, j::Int)
    @boundscheck checkbounds(V, k, j)
    A = parent(V)
    K,J = parentblocks2Int(V)
    if -A.l ≤ J-K ≤ A.u && -A.λ ≤ j-k ≤ A.μ
        inbands_setindex!(V, v, k, j)
    elseif iszero(v) # allow setindex for 0 datya
        v
    else
        throw(BandError(V, J-K))
    end
end

BLAS.axpy!(a::T, A::BandedBlockBandedMatrix{T}, B::BandedBlockBandedMatrix{T}) where T =
    B .= a .* A .+ B

lmul!(x::Number, A::BandedBlockBandedMatrix) = (lmul!(x, A.data); A)
rmul!(A::BandedBlockBandedMatrix, x::Number) = (rmul!(A.data, x); A)
*(x::Number, A::BandedBlockBandedMatrix) = _BandedBlockBandedMatrix(x*A.data, axes(A,1), blockbandwidths(A), subblockbandwidths(A))
*(A::BandedBlockBandedMatrix, x::Number) = _BandedBlockBandedMatrix(A.data*x, axes(A,1), blockbandwidths(A), subblockbandwidths(A))


text:
using HTTP

export start_runtime

Base.@kwdef struct Invocation
  body::Any
  aws_request_id::String
  deadline_ms::Int
  invoked_function_arn::String
  trace_id::String
end

function get_endpoint(host::String)::String
  "http://$host/2018-06-01/runtime/invocation/"
end

function lambda_respond(response::String, endpoint::String, aws_request_id::String)
  HTTP.request(
    "POST",
    "$(endpoint)$(aws_request_id)/response",
    [],
    response,
  )
end

function lambda_error(error::String, endpoint::String, aws_request_id::String)
  HTTP.request(
    "POST",
    "$(endpoint)$(aws_request_id)/error",
    [("Lambda-Runtime-Function-Error-Type", "Unhandled")],
    JSON3.write(error),
  )
end

function start_runtime(host::String, func_name::String, param_type::String; single_shot=false)
  param_type = eval(Meta.parse(param_type))
  start_runtime(host, eval(Meta.parse(func_name)), param_type)
end

function start_runtime(host::String, react_function::Function, ::Type{T}; single_shot=false) where {T}
  endpoint = get_endpoint(host)
  println("Starting runtime at $endpoint")

  while true
    http = HTTP.request("GET", "$(endpoint)next"; verbose=3)
    body_raw = String(http.body)
    request_id = string(HTTP.header(http, "Lambda-Runtime-Aws-Request-Id"))

    body = try
      JSON3.read(body_raw, T)
    catch e
      body_sample = length(body_raw) > 50 ? "$(body_raw[start:50])..." : body_raw
      lambda_error("Unable to parse input JSON $body_sample", endpoint, request_id)
      continue
    end

    reaction = try
      react_function(body)
    catch e
      err(msg) = lambda_error(msg, endpoint, request_id)
      if isa(e, MethodError) && e.f == String(Symbol(react_function))
        err("react function is not a valid method for parsed JSON type")
      else
        err("$e")
      end
      continue
    end

    reaction_json = try
      JSON3.write(reaction)
    catch e
      err("Unable to parse function return value $(reaction) to JSON")
      continue
    end

    lambda_respond(reaction_json, endpoint, request_id)
    single_shot && break
  end
end

ext: jl


text:
using Distributions
using Base.Test
import Distributions: PoissonADSampler, PoissonCountSampler

function test_poissonsampler(s, μ::Float64, ns::Int, tol::Float64)
    ub = iceil(μ * 3) + 1
    pv = Distributions.poissonpvec(μ, ub)
    @assert length(pv) == ub + 1

    cnts = zeros(Int, ub+1)
    for i = 1:ns
        x = rand(s)
        @assert x >= 0
        if x <= ub
            cnts[x + 1] += 1
        end
    end

    q = cnts ./ ns
    @test_approx_eq_eps q pv tol
end

# test cases

for μ in [0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 15.0, 20.0, 30.0]
    test_poissonsampler(PoissonCountSampler(μ), μ, 10^5, 0.02)
end

for μ in [5.0, 10.0, 15.0, 20.0, 30.0]
    test_poissonsampler(PoissonADSampler(μ), μ, 10^5, 0.02)
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using HTTP
using JSON
using SHA

# =========================================================================================
"""
    download_file(
        fromUrl::AbstractString,
        toFile::AbstractString
    )::String

Download the the URL and safe it to the file.  Returns the GitHub SHA1 value for the file
that was downloaded.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | fromUrl  |         | The URL to download.
| R | toFile   |         | The file to save the contents to.

__Exceptions:__

| Name         | Reason
| :----------- | :-----
| FileError    | There was an issue writing to the file.
| NetworkError | There was an issue downloading the file.

__Details:__

If the file exists it is overwritten.  The GitHub SHA1 is slightly different that a sha1sum
calculated locally.  First "blob $size\0" is hashed followed by the contents of the file.
"""
function download_file(
            fromUrl::AbstractString,
            toFile::AbstractString
        )::String
    local sha1 = SHA1_CTX()

    try
        # ---------------------------------------------------------------------------------
        # Open the file for writing, and make the HTTP request.
        open(toFile; create = true, write = true) do output

            HTTP.open("GET", fromUrl; redirect = false) do incoming

                # -------------------------------------------------------------------------
                # Start reading the data.
                startread(incoming)

                # -------------------------------------------------------------------------
                # Check for an error.
                if incoming.message.status != 200
                    @warn "Failed to download $fromUrl, HTTP error: $(incoming.message.status)"
                    throw(NetworkError(fromUrl))
                end

                # -------------------------------------------------------------------------
                # Start the hash.
                local size = parse(Int, HTTP.header(incoming, "Content-Length"))
                update!(sha1, Vector{UInt8}("blob $size\0"))

                # -------------------------------------------------------------------------
                # Download the file.
                local buffer = Vector{UInt8}(undef, 8192)

                while eof(incoming) == false
                    local bytes = readbytes!(incoming, buffer)
                    write(output, view(buffer, 1:bytes))
                    update!(sha1, buffer, bytes)
                end
            end
        end
    catch ex
        # ---------------------------------------------------------------------------------
        # We had an error nuke the file.
        if isfile(toFile)
            rm(toFile)
        end

        # ---------------------------------------------------------------------------------
        # SystemError means we had an file error.
        if isa(ex, SystemError)
            throw(FileError(toFile))
        end

        # ---------------------------------------------------------------------------------
        # Unknown error, rethrow it.
        rethrow(ex)
    end

    # -------------------------------------------------------------------------------------
    # Return the DataFile object.
    return bytes2hex(digest!(sha1))
end

# =========================================================================================
"""
    download_json(
        url::AbstractString,
    )::String

Download the the URL and parse the response as JSON.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | url  |         | The URL to download.

__Exceptions:__

| Name         | Reason
| :----------- | :-----
| NetworkError | There was an issue downloading the JSON object.
"""
function download_json(
            url::AbstractString
        )::Union{AbstractDict, AbstractArray}
    local response = HTTP.request("GET",  url)

    if response.status != 200
        @warn "Failed to download $url, HTTP error: $(response.status)"
        throw(NetworkError(url))
    end

    local raw = String(response.body)
    local json

    try
        json = JSON.parse(raw)
    catch ex
        if isa(ex, ErrorException)
            @warn "Unexpected response from the server: $raw"
            throw(NetworkError(url))
        end
        rethrow(ex)
    end

    return json
end


text:
# pick a nice default x range given a distribution
function default_range(dist::Distribution, alpha = 0.0001)
    minval = isfinite(minimum(dist)) ? minimum(dist) : quantile(dist, alpha)
    maxval = isfinite(maximum(dist)) ? maximum(dist) : quantile(dist, 1-alpha)
    minval, maxval
end

function default_range(m::Distributions.MixtureModel, alpha = 0.0001)
    minval = maxval = 0.0
    for c in m.components
        thismin = isfinite(minimum(c)) ? minimum(c) : quantile(c, alpha)
        thismax = isfinite(maximum(c)) ? maximum(c) : quantile(c, 1-alpha)
        if thismin < minval
            minval = thismin
        end
        if thismax > maxval
            maxval = thismax
        end
    end
    minval, maxval
end

yz_args(dist) = default_range(dist)
yz_args(dist::Distribution{N, T}) where N where T<:Discrete = (UnitRange(default_range(dist)...),)

# this "user recipe" adds a default x vector based on the distribution's μ and σ
@recipe function f(dist::Distribution)
    if dist isa Distribution{Univariate,Discrete}
        seriestype --> :scatterpath
    end
    (dist, yz_args(dist)...)
end

@recipe function f(m::Distributions.MixtureModel; components = true)
    if components
        for c in m.components
            @series begin
                (c, yz_args(c)...)
            end
        end
    else
        (m, yz_args(m)...)
    end
end

@recipe function f(distvec::AbstractArray{<:Distribution}, yz...)
    for di in distvec
        @series begin
            seriesargs = isempty(yz) ? yz_args(di) : yz
            if di isa Distribution{Univariate,Discrete}
                seriestype --> :scatterpath
            end
            (di, seriesargs...)
        end
    end
end

# this "type recipe" replaces any instance of a distribution with a function mapping xi to yi
@recipe function f(::Type{T}, dist::T; func = pdf) where T<:Distribution
    xi -> func(dist, xi)
end

#-----------------------------------------------------------------------------
# qqplots

@recipe function f(h::QQPair; qqline = :identity)
    if qqline in (:fit, :quantile, :identity, :R)
        xs = [extrema(h.qx)...]
        if qqline == :identity
            ys = xs
        elseif qqline == :fit
            itc, slp = hcat(fill!(similar(h.qx), 1), h.qx) \ h.qy
            ys = slp .* xs .+ itc
        else # if qqline == :quantile || qqline == :R
            quantx, quanty = quantile(h.qx, [0.25, 0.75]), quantile(h.qy, [0.25, 0.75])
            slp = diff(quanty) ./ diff(quantx)
            ys = quanty .+ slp .* (xs .- quantx)
        end

        @series begin
            primary := false
            seriestype := :path
            xs, ys
        end
    end

    seriestype --> :scatter
    legend --> false
    h.qx, h.qy
end

loc(D::Type{T}, x) where T<:Distribution = fit(D, x), x
loc(D, x) = D, x

@userplot QQPlot
recipetype(::Val{:qqplot}, args...) = QQPlot(args)
@recipe f(h::QQPlot) = qqbuild(loc(h.args[1], h.args[2])...)

@userplot QQNorm
recipetype(::Val{:qqnorm}, args...) = QQNorm(args)
@recipe f(h::QQNorm) = QQPlot((Normal, h.args[1]))


text:
abstract AbstractFilter{T<:Number,D<:AbstractDomain}

immutable Asymmetric1DFilter{T<:Number} <: AbstractFilter{T,SpatialDomain{1}}
    neg::Vector{T}
    pos::Vector{T}
    zero::T
end

immutable Symmetric1DFilter{T<:Number} <: AbstractFilter{T,SpatialDomain{1}}
    leg::Vector{T}
    zero::T
end

# product with scalar number is commutative
Base.:*(b::Number, ψ::AbstractFilter) = ψ * b

# element-wise multiplication operator ".*" with scalar falls back to "*"
Base.:.*(ψ::AbstractFilter, b::Number) = ψ * b

function spin!{T,D}(ψs::Array{AbstractFilter{T,D},3})
    (nΘs, nΧs, nJs) = size(ψs)
    for χ in 0:(nΧs-1), j in 0:(nJs-1)
        ψs[2:end, 1+χ, 1+j] = spin(ψs[1, 1+χ, 1+j])
    end
end


text:
function get_commandline_applescript(exe_dir, executable_path)
    """
    set RootPath to POSIX path of (path to me)
    tell application id "com.apple.terminal"
      do script ("exec '" & RootPath & "$exe_dir/$executable_path'")
      activate
    end tell
    """
end

# Build a wrapper app that opens a terminal and runs the provided binary.
# Returns a new applet name only if binary_name is already "applet" (to prevent collision).
function build_commandline_app_bundle(builddir, binary_name, appname, verbose)
    println("~~~~~~ Creating commandline-app wrapper applet. ~~~~~~~")

    mkpath(builddir)  # Create builddir if it doesn't already exist.
    applescript_file = "$builddir/build_$appname.applescript"
    exe_dir = "Contents/MacOS"  # Put the binaries next to the applet in MacOS.
    verbose && println("    Creating Applescript wrapper file: $applescript_file")
    write(applescript_file, get_commandline_applescript(exe_dir, binary_name))

    app_path = "$builddir/$appname.app"
    run(`osacompile -o $app_path $applescript_file`)
    applet_name = "applet"
    if binary_name == applet_name  # Prevent collisions.
        applet_name = "applet_wrapper"
        Compat.mv("$app_path/Contents/MacOS/applet", "$app_path/Contents/MacOS/$applet_name", force=true)
    end
    # Remove unneeded applet files it creates.
    rm("$app_path/Contents/Resources/applet.icns")

    # Cleaning up temporary file:
    verbose && println("    Cleaning up temporary Applescript file: $applescript_file")
    rm(applescript_file)

    return applet_name
end


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

using Test
import Libdl

# these could fail on an embedded installation
# but for now, we don't handle that case
dlls = Libdl.dllist()
@test !isempty(dlls)
@test length(dlls) > 3 # at a bare minimum, probably have some version of libstdc, libgcc, libjulia, ...
if !Sys.iswindows() || Sys.windows_version() >= Sys.WINDOWS_VISTA_VER
    for dl in dlls
        if isfile(dl) && (Libdl.dlopen(dl; throw_error=false) !== nothing)
            @test Base.samefile(Libdl.dlpath(dl), dl)
        end
    end
end
@test length(filter(dlls) do dl
      if Base.DARWIN_FRAMEWORK
          return occursin(Regex("^$(Base.DARWIN_FRAMEWORK_NAME)(?:_debug)?\$"), basename(dl))
      else
          return occursin(Regex("^libjulia(?:.*)\\.$(Libdl.dlext)(?:\\..+)?\$"), basename(dl))
      end
    end) == 1 # look for something libjulia-like (but only one)

# library handle pointer must not be NULL
@test_throws ArgumentError Libdl.dlsym(C_NULL, :foo)
@test_throws ArgumentError Libdl.dlsym_e(C_NULL, :foo)

cd(@__DIR__) do

# Find the library directory by finding the path of libjulia (or libjulia-debug, as the case may be)
# and then adding on /julia to that directory path to get the private library directory, if we need
# to (where "need to" is defined as private_libdir/julia/libccalltest.dlext exists
private_libdir = if Base.DARWIN_FRAMEWORK
    if ccall(:jl_is_debugbuild, Cint, ()) != 0
        dirname(abspath(Libdl.dlpath(Base.DARWIN_FRAMEWORK_NAME * "_debug")))
    else
        joinpath(dirname(abspath(Libdl.dlpath(Base.DARWIN_FRAMEWORK_NAME))),"Frameworks")
    end
elseif ccall(:jl_is_debugbuild, Cint, ()) != 0
    dirname(abspath(Libdl.dlpath("libjulia-debug")))
else
    dirname(abspath(Libdl.dlpath("libjulia")))
end

if isfile(joinpath(private_libdir,"julia","libccalltest."*Libdl.dlext))
    private_libdir = joinpath(private_libdir, "julia")
end

@test !isempty(Libdl.find_library(["libccalltest"], [private_libdir]))
@test !isempty(Libdl.find_library("libccalltest", [private_libdir]))
@test !isempty(Libdl.find_library(:libccalltest, [private_libdir]))

# dlopen should be able to handle absolute and relative paths, with and without dlext
let dl = C_NULL
    try
        dl = Libdl.dlopen(abspath(joinpath(private_libdir, "libccalltest")); throw_error=false)
        @test dl !== nothing
    finally
        Libdl.dlclose(dl)
    end
end

let dl = C_NULL
    try
        dl = Libdl.dlopen(abspath(joinpath(private_libdir, "libccalltest.$(Libdl.dlext)")); throw_error=false)
        @test dl !== nothing
    finally
        Libdl.dlclose(dl)
    end
end

let dl = C_NULL
    try
        dl = Libdl.dlopen(relpath(joinpath(private_libdir, "libccalltest")); throw_error=false)
        @test dl !== nothing
    finally
        Libdl.dlclose(dl)
    end
end

let dl = C_NULL
    try
        dl = Libdl.dlopen(relpath(joinpath(private_libdir, "libccalltest.$(Libdl.dlext)")); throw_error=false)
        @test dl !== nothing
    finally
        Libdl.dlclose(dl)
    end
end

let dl = C_NULL
    try
        dl = Libdl.dlopen("./foo"; throw_error=false)
        @test dl === nothing
    finally
        Libdl.dlclose(dl)
    end
end

# unqualified names present in DL_LOAD_PATH
let dl = C_NULL
    try
        dl = Libdl.dlopen("libccalltest"; throw_error=false)
        @test dl !== nothing
    finally
        Libdl.dlclose(dl)
    end
end

let dl = C_NULL
    try
        dl = Libdl.dlopen(string("libccalltest",".",Libdl.dlext); throw_error=false)
        @test dl !== nothing
    finally
        Libdl.dlclose(dl)
    end
end

# path with dlopen-able file first in load path
#=
let dl = C_NULL,
    tmpdir = mktempdir(),
    fpath = joinpath(tmpdir,"libccalltest")
    try
        write(open(fpath,"w"))
        push!(Libdl.DL_LOAD_PATH, @__DIR__)
        push!(Libdl.DL_LOAD_PATH, dirname(fpath))
        dl = Libdl.dlopen_e("libccalltest")
        @test dl != C_NULL
    finally
        pop!(Libdl.DL_LOAD_PATH)
        pop!(Libdl.DL_LOAD_PATH)
        rm(tmpdir, recursive=true)
    end
end
=#

# path with dlopen-able file second in load path
#=
let dl = C_NULL,
    tmpdir = mktempdir(),
    fpath = joinpath(tmpdir,"libccalltest")
    try
        write(open(fpath,"w"))
        push!(Libdl.DL_LOAD_PATH, dirname(fpath))
        push!(Libdl.DL_LOAD_PATH, @__DIR__)
        dl = Libdl.dlopen_e("libccalltest")
        @test dl != C_NULL
    finally
        pop!(Libdl.DL_LOAD_PATH)
        pop!(Libdl.DL_LOAD_PATH)
        rm(tmpdir, recursive=true)
    end
end
=#

# test dlpath
let dl = C_NULL
    try
        path = abspath(joinpath(private_libdir, "libccalltest"))
        dl = Libdl.dlopen(path; throw_error=false)
        @test dl !== nothing
        @test Base.samefile(abspath(Libdl.dlpath(dl)),
                            abspath(Libdl.dlpath(path)))
        @test Base.samefile(abspath(Libdl.dlpath(dl)),
                            string(path,".",Libdl.dlext))
    finally
        Libdl.dlclose(dl)
    end
end

# opening a library that does not exist throws an ErrorException
@test_throws ErrorException Libdl.dlopen("./foo")

# opening a versioned library that does not exist does not result in adding extension twice
err = @test_throws ErrorException Libdl.dlopen("./foo.$(Libdl.dlext).0")
@test !occursin("foo.$(Libdl.dlext).0.$(Libdl.dlext)", err.value.msg)
err = @test_throws ErrorException Libdl.dlopen("./foo.$(Libdl.dlext).0.22.1")
@test !occursin("foo.$(Libdl.dlext).0.22.1.$(Libdl.dlext)", err.value.msg)

# test dlsym
let dl = C_NULL
    try
        dl = Libdl.dlopen(abspath(joinpath(private_libdir, "libccalltest")))
        fptr = Libdl.dlsym(dl, :set_verbose)
        @test fptr !== nothing
        @test_throws ErrorException Libdl.dlsym(dl, :foo)

        fptr = Libdl.dlsym_e(dl, :set_verbose)
        @test fptr != C_NULL
        fptr = Libdl.dlsym_e(dl, :foo)
        @test fptr == C_NULL
    finally
        Libdl.dlclose(dl)
    end
end

# test dlclose
# If dl is NULL, jl_dlclose should return -1 and dlclose should return false
# dlclose should return true on success and false on failure
let dl = C_NULL
    @test -1 == ccall(:jl_dlclose, Cint, (Ptr{Cvoid},), dl)
    @test !Libdl.dlclose(dl)

    dl = Libdl.dlopen("libccalltest"; throw_error=false)
    @test dl !== nothing

    @test Libdl.dlclose(dl)
    @test_skip !Libdl.dlclose(dl)   # Syscall doesn't fail on Win32
end

end

max_stars_repo_name: musm/julia

max_issues_repo_name: musm/julia

max_forks_repo_name: musm/julia


text:
function g(x) x = 42 end
function f(x) x = 4; g(x); println(x) end
f(3)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
#=
    Julia is a bit even more snappy when indicies, bitshifts and other register-centric ops
        are performed with Int16 or Int32 values.  Int64 values are very fast, too.
     The functions I use to optimize the type of the indexing value are needed
        when they are being leaned upon to provide time-minimized access.
=#

@inline bitsof(::Type{T}) where T = sizeof(T) * 8
bitsof(x::T) where T = bitsof(T)

function bits_required(posint::T)::T where T<:Integer
    posint > 0 || throw(ErrorException("posint ($(posint)) must be > 0"))
    return bitsof(T) - leading_zeros(posint)%T
end

const TYPES_FOR_INDEXING = [ Int8, Int16, Int32, Int64, Int128 ]

function type_for_indexing(posint::T) where T<:Union{Int8, Int16, Int32, Int64, Int128}
    posint < 0 && throw(ErrorException("posint ($(posint)) must be => 0"))
    posint = max(127%T, posint)
    type_index = nextpow2(bits_required(posint) + 1) >>> 2
    type_index = trailing_zeros(type_index)
    type_index = max(1, min(5, type_index))
    return TYPES_FOR_INDEXING[ type_index ]
end


text:
using PyPlot
include("style_colorscheme.jl")
using3D()

DrWatson._wsave(s, fig::Figure) = fig.savefig(s, dpi = 600, transparent = false)

"""
    rdspl(x, n = 2)
Round `x` to `n` sigdigits for display purposes.
"""
rdspl(x::Real, n = 2) = round(x, digits=n)
rdspl(x::AbstractVector, n = 2) = Tuple((round.(Float64.(x); sigdigits=n)))

PyPlot.rc("lines", lw = 3.0)
PyPlot.rc("errorbar", capsize = 6)
PyPlot.rc("axes", grid = true)
PyPlot.rc("grid", color = "0.75", alpha = 0.75)
PyPlot.rc("axes", axisbelow=true) # makes grid behind other plotting elements

PyPlot.rc("font", size = 24) # set default fontsize
PyPlot.rc("xtick", labelsize = 24)
PyPlot.rc("ytick", labelsize = 24)
PyPlot.rc("axes", labelsize = 28)
PyPlot.rc("legend";
    fontsize = 24, handletextpad = 0.5, handlelength = 2,
    labelspacing = 0.25, columnspacing = 1.0,
)
# PyPlot.rc("font", family = "Times New Roman") # Serif main font
PyPlot.rc("font", family = "DejaVu Sans") # sans main font
# PyPlot.rc("mathtext", rm = "sanserif", fontset="dejavusans") # sans math font
PyPlot.rc("mathtext", rm = "serif", fontset="cm") # serif math font

for z in ("x", "y")
    PyPlot.rc("$(z)tick.major", size = 7, width = 1.2)
    PyPlot.rc("$(z)tick.minor", size = 3, visible = false)
end

figx = 12 # default width correspoding to full text width
figy = 9  # default height corresponding to 1 row of plots
PyPlot.rc("figure", figsize = (figx, figy))
PyPlot.rc("savefig", dpi = 600, transparent = true, format = "png")

# set default color cycle
PyPlot.rc("axes", prop_cycle = matplotlib.cycler(color=COLORSCHEME))

if false # test color scheme
    fig = figure(figsize = (20, 15)) # show colors
    ax1 = subplot(231)
    ax2 = subplot(232)
    ax3 = subplot(233)
    ax4 = subplot(223)
    ax5 = subplot(224)
    lw = 60
    L = length(COLORSCHEME)
    for (i, c) in enumerate(COLORS)
        chsv = matplotlib.colors.rgb_to_hsv(matplotlib.colors.to_rgb(c))
        ax1.plot([0, 1], [0, 0] .+ i, color = c, lw = lw)
        ax1.set_title("color")
        ax2.plot([0, 1], [0, 0] .+ i, color = string(chsv[3]), lw = lw)
        ax2.set_title("brightness")
        ax3.plot([0, 1], [0, 0] .+ i, color = string(chsv[2]), lw = lw)
        ax3.set_title("saturation")
        x = 0:0.05:5π
        ax4.plot(x, cos.(x .+ i/2) .+ rand(length(x))/2; color=c, lw = 2)
        ax5.bar(collect(1:4) .+ (i-1)/L, 0.5rand(4) .+ 0.5, 1/L; color=c)
    end
    fig = tight_layout()
end

bbox = Dict(:boxstyle => "round,pad=0.3", :facecolor=>"white", :alpha => 1.0)

"`add_identifiers!(fig = gcf(), axs = fig.get_axes(); xloc = 0.985, yloc = 0.975)`"
function add_identifiers!(fig = gcf(), axs = fig.get_axes(); xloc = 0.975, yloc = 1)
    bbox = Dict(:boxstyle => "round,pad=0.3", :facecolor=>"white", :alpha => 1.0)
    for (i, ax) in enumerate(axs)
        l = ('a':'z')[i]
        try
            ax.text(xloc, yloc, "$(l)"; transform = ax.transAxes,
            bbox = bbox, zorder = 99, va = "top")
        catch e
            @warn string(e)
        end
    end
end

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
#=
https://www.juliabox.com/notebook/notebooks/tutorials/intro-to-julia/04.%20Loops.ipynb
=#

#=
let's use for loops to create some addition tables,
where the value of every entry is the sum of its row and column indices.
=#

function create_table1(m,n)
    A = fill(0, (m, n))
    for i in 1:m
        for j in 1:n
            A[i,j] = i+j
        end
    end
    return A
end

#=
Here is some syntactic sugar for the same nested `for` loop
=#

function create_table2(m,n)
    B = fill(0, (m, n))
    for i in 1:m, j in 1:n
        B[i,j] = i+j
    end
    return B
end

#=
The more "Julia" way to create this addition table would have been
with an array comprehension.
=#

function create_table3(m,n)
    C = [i+j for i in 1:m, j in 1:n]
    return C
end

@show create_table1(4,5)
@show create_table2(4,5)
@show create_table3(4,5)

using InteractiveUtils
versioninfo()

using BenchmarkTools
m,n=1000,1000
@btime create_table1(m,n)
@btime create_table2(m,n)
@btime create_table3(m,n)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
cd("lesson-packages")
using Pkg
Pkg.activate("TestPackage1")

#%%
using TestPackage1

TestPackage1.greet()
#%%
mySpecialFunction(3)
#%%
using Pkg
Pkg.instantiate()

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples


text:
struct ∇{T}
  operand::T
  string::String
end

struct ∇o{T}
  operand::T
  string::String
end

struct ∇x{T}
  operand::T
  string::String
end

String(x::∇o) = "∇o" * x.string
String(x::∇x) = "∇×" * x.string

(g::∇)(x) = CovariantVector(x->grad(g.operand, x), true)
(c::∇x)(x) = ContravariantVector(x->curl(c.operand, x), true)
(d::∇o)(v::VectorField) = Scalar(dot(d.operand, v))
(d::∇o)(c::∇x) = 0
(c::∇x)(g::∇) = zeros(3)

Base.show(io::IO, x::∇o) = print(io, String(x))
Base.show(io::IO, x::∇x) = print(io, String(x))

max_stars_repo_name: jwscook/DifferentialGeometries

max_issues_repo_name: jwscook/DifferentialGeometries

max_forks_repo_name: jwscook/DifferentialGeometries


text:
function save_features_instance!(db::SQLite.DB, name, scenario, source_type, source_path)
    println("Saving instance features: ($name, $scenario)")
    network = PowerFlowNetwork(source_path, source_type)
    merge_duplicate_branch!(network)
    features = get_features_instance(network)
    query = "UPDATE instances SET "
    for (feature_name, feature_value) in features
        if feature_value == Inf
            query *= "$feature_name = '+Infinity', "
        elseif feature_value == -Inf
            query *= "$feature_name = '-Infinity', "
        elseif isnan(feature_value)
            query *= "$feature_name = 'NaN', "
        else
            query *= "$feature_name = $feature_value, "
        end
    end
    query = query[begin:end - 2]
    query *= " WHERE name = '$name' AND scenario = $scenario"
    execute_query(db, query)
end

function save_features_instance_dfrow!(db::SQLite.DB, row)
    save_features_instance!(db, row[:name], row[:scenario], row[:source_type], row[:source_path])
end

function save_features_instances!(db::SQLite.DB; min_nv=typemin(Int), max_nv=typemax(Int), recompute=false, subset=nothing)
    query = "SELECT name, scenario, source_type, source_path FROM instances WHERE nbus >= $min_nv AND nbus <= $max_nv"
    if !recompute
        query *= " AND nb_edge IS NULL"
    end
    if !isnothing(subset)
        query *= " AND id IN ($(join(subset, ',')))"
    end
    results = DBInterface.execute(db, query) |> DataFrame
    println("Saving instance features config: min_nv=$min_nv, max_nv=$max_nv, recompute=$recompute")
    save_function!(row) = save_features_instance_dfrow!(db, row)
    save_function!.(eachrow(results[!, [:name, :scenario, :source_type, :source_path]]))
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using CUDA
using OrderedCollections: OrderedDict

using Oceananigans.Advection
using Oceananigans: AbstractOutputWriter, AbstractDiagnostic, TimeStepper
using Oceananigans.Architectures: AbstractArchitecture
using Oceananigans.Buoyancy: validate_buoyancy
using Oceananigans.TurbulenceClosures: ν₀, κ₀, with_tracers

mutable struct IncompressibleModel{TS, E, A<:AbstractArchitecture, G, T, B, R, SW, U, C, Φ, F,
                                   V, S, K} <: AbstractModel
       architecture :: A         # Computer `Architecture` on which `Model` is run
               grid :: G         # Grid of physical points on which `Model` is solved
              clock :: Clock{T}  # Tracks iteration number and simulation time of `Model`
          advection :: V         # Advection scheme for velocities _and_ tracers
           buoyancy :: B         # Set of parameters for buoyancy model
           coriolis :: R         # Set of parameters for the background rotation rate of `Model`
      surface_waves :: SW        # Set of parameters for surfaces waves via the Craik-Leibovich approximation
            forcing :: F         # Container for forcing functions defined by the user
            closure :: E         # Diffusive 'turbulence closure' for all model fields
         velocities :: U         # Container for velocity fields `u`, `v`, and `w`
            tracers :: C         # Container for tracer fields
          pressures :: Φ         # Container for hydrostatic and nonhydrostatic pressure
      diffusivities :: K         # Container for turbulent diffusivities
        timestepper :: TS        # Object containing timestepper fields and parameters
    pressure_solver :: S         # Pressure/Poisson solver
end

"""
    IncompressibleModel(;
                   grid,
           architecture = CPU(),
             float_type = Float64,
                  clock = Clock{float_type}(0, 0, 1),
              advection = CenteredSecondOrder(),
               buoyancy = SeawaterBuoyancy(float_type),
               coriolis = nothing,
          surface_waves = nothing,
                forcing = ModelForcing(),
                closure = IsotropicDiffusivity(float_type, ν=ν₀, κ=κ₀),
    boundary_conditions = (u=UVelocityBoundaryConditions(grid),
                           v=VVelocityBoundaryConditions(grid),
                           w=WVelocityBoundaryConditions(grid)),
                tracers = (:T, :S),
            timestepper = :QuasiAdamsBashforth2,
             velocities = VelocityFields(architecture, grid, boundary_conditions),
              pressures = PressureFields(architecture, grid, boundary_conditions),
          diffusivities = DiffusivityFields(architecture, grid, tracernames(tracers), boundary_conditions, closure),
        pressure_solver = PressureSolver(architecture, grid, PressureBoundaryConditions(grid))
    )

Construct an incompressible `Oceananigans.jl` model on `grid`.

Keyword arguments
=================

    - `grid`: (required) The resolution and discrete geometry on which `model` is solved.
    - `architecture`: `CPU()` or `GPU()`. The computer architecture used to time-step `model`.
    - `float_type`: `Float32` or `Float64`. The floating point type used for `model` data.
    - `advection`: The scheme that advects velocities and tracers. See `Oceananigans.Advection`.
    - `buoyancy`: The buoyancy model. See `Oceananigans.Buoyancy`.
    - `closure`: The turbulence closure for `model`. See `Oceananigans.TurbulenceClosures`.
    - `coriolis`: Parameters for the background rotation rate of the model.
    - `forcing`: User-defined forcing functions that contribute to solution tendencies.
    - `boundary_conditions`: `NamedTuple` containing field boundary conditions.
    - `tracers`: A tuple of symbols defining the names of the modeled tracers, or a `NamedTuple` of
                 preallocated `CellField`s.
    - `timestepper`: A symbol that species the time-stepping method. Either `:QuasiAdamsBashforth2` or
                     `:RungeKutta3`.
"""
function IncompressibleModel(;
                   grid,
           architecture = CPU(),
             float_type = Float64,
                  clock = Clock{float_type}(0, 0, 1),
              advection = CenteredSecondOrder(),
               buoyancy = SeawaterBuoyancy(float_type),
               coriolis = nothing,
          surface_waves = nothing,
                forcing = ModelForcing(),
                closure = IsotropicDiffusivity(float_type, ν=ν₀, κ=κ₀),
    boundary_conditions = (u=UVelocityBoundaryConditions(grid),
                           v=VVelocityBoundaryConditions(grid),
                           w=WVelocityBoundaryConditions(grid)),
                tracers = (:T, :S),
            timestepper = :QuasiAdamsBashforth2,
             velocities = VelocityFields(architecture, grid, boundary_conditions),
              pressures = PressureFields(architecture, grid, boundary_conditions),
          diffusivities = DiffusivityFields(architecture, grid, tracernames(tracers), boundary_conditions, closure),
        pressure_solver = PressureSolver(architecture, grid, PressureBoundaryConditions(grid))
    )

    if architecture == GPU() && !has_cuda()
         throw(ArgumentError("Cannot create a GPU model. No CUDA-enabled GPU was detected!"))
    end

    validate_buoyancy(buoyancy, tracernames(tracers))

    # Regularize forcing and closure for given tracer fields.
    forcing = ModelForcing(tracernames(tracers), forcing)
    closure = with_tracers(tracernames(tracers), closure)

    # Instantiate tracer fields if not already instantiated
    tracer_fields = TracerFields(architecture, grid, tracers, boundary_conditions)

    # Instantiate timestepper if not already instantiated
    timestepper = TimeStepper(timestepper, architecture, grid, tracernames(tracers))

    return IncompressibleModel(architecture, grid, clock, advection, buoyancy, coriolis, surface_waves,
                               forcing, closure, velocities, tracer_fields, pressures, diffusivities,
                               timestepper, pressure_solver)
end


text:
#---------------------------------------------------------
# # [SPECTrecon PSF](@id 3-psf)
#---------------------------------------------------------

# This page explains the PSF portion of the Julia package
# [`SPECTrecon.jl`](https://github.com/JeffFessler/SPECTrecon.jl).

# ### Setup

# Packages needed here.

using SPECTrecon
using MIRTjim: jim, prompt
using Plots: scatter, scatter!, plot!, default
default(markerstrokecolor=:auto, markersize=3)

# The following line is helpful when running this example.jl file as a script;
# this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);

# ### Overview

# After rotating the image and the attenuation map,
# second step in SPECT image forward projection
# is to apply depth-dependent point spread function (PSF).
# Each (rotated) image plane
# is a certain distance from the SPECT detector
# and must be convolved with the 2D PSF appropriate
# for that plane.

# Because SPECT has relatively poor spatial resolution,
# the PSF is usually fairly wide,
# so convolution using FFT operations
# is typically more efficient
# than direct spatial convolution.

# Following other libraries like
# [FFTW.jl](https://github.com/JuliaMath/FFTW.jl),
# the PSF operations herein start with a `plan`
# where work arrays are preallocated
# for subsequent use.
# The `plan` is a `Vector` of `PlanPSF` objects:
# one for each thread.
# (Parallelism is across planes for a 3D image volume.)
# The number of threads defaults to `Threads.nthreads()`.

# ### Example

# Start with a 3D image volume.

T = Float32 # work with single precision to save memory
nx = 32
nz = 30
image = zeros(T, nx, nx, nz) # ny = nx required
image[1nx÷4, 1nx÷4, 3nz÷4] = 1
image[2nx÷4, 2nx÷4, 2nz÷4] = 1
image[3nx÷4, 3nx÷4, 1nz÷4] = 1
jim(image, "Original image")

# Create a synthetic gaussian depth-dependent PSF for a single view

px = 11
nview = 1 # for simplicity in this illustration
psf = repeat(psf_gauss( ; ny, px), 1, 1, 1, nview)
jim(psf, "PSF for each of $nx planes")

# Now plan the PSF modeling
# by specifying
# * the image size (must be square)
# * the PSF size: must be `px × pz × ny × nview`
# * the `DataType` used for the work arrays.

plan = plan_psf( ; nx, nz, px, T)

# Here are the internals for the plan for the first thread:

plan[1]

# With this `plan` pre-allocated, now we can apply the depth-dependent PSF
# to the image volume (assumed already rotated here).

result = similar(image) # allocate memory for the result
fft_conv!(result, image, psf[:,:,:,1], plan) # mutates the first argument
jim(result, "After applying PSF")

# ### Adjoint

# To ensure adjoint consistency between SPECT forward- and back-projection,
# there is also an adjoint routine:

adj = similar(result)
fft_conv_adj!(adj, result, psf[:,:,:,1], plan)
jim(adj, "Adjoint of PSF modeling")

# The adjoint is *not* the same as the inverse
# so one does not expect the output here to match the original image!

# ### LinearMap

# One can form a linear map corresponding to PSF modeling using `LinearMapAA`.
# Perhaps the main purpose is simply for verifying adjoint correctness.

using LinearMapsAA: LinearMapAA

nx, nz, px = 10, 7, 5 # small size for illustration
psf3 = psf_gauss( ; ny, px)
plan = plan_psf( ; nx, nz, px, T)
idim = (nx,nx,nz)
odim = (nx,nx,nz)
forw! = (y,x) -> fft_conv!(y, x, psf3, plan)
back! = (x,y) -> fft_conv_adj!(x, y, psf3, plan)
A = LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, odim, idim)

Afull = Matrix(A)
Aadj = Matrix(A')
jim(cat(dims=3, Afull, Aadj'), "Linear map for PSF modeling and its adjoint")

# The following check verifies adjoint consistency:

@assert Afull ≈ Aadj'


text:
# https://adventofcode.com/2020/day/12

#  * Action `N` means to move north by the given value.
#  * Action `S` means to move south by the given value.
#  * Action `E` means to move east by the given value.
#  * Action `W` means to move west by the given value.
#  * Action `L` means to turn left the given number of degrees.
#  * Action `R` means to turn right the given number of degrees.
#  * Action `F` means to move forward by the given value in the
#    direction the ship is currently facing.

example = "
F10
N3
F7
R90
F11
"

const directions = (:North, :East, :South, :West)
const char_to_direction = Dict("N" => :North, "E" => :East, "S" => :South, "W" => :West)

Point = Tuple{Integer,Integer}

mutable struct Ship
    facing::Symbol
    position::Point
    waypoint::Point
end

facing(ship::Ship) = ship.facing
position(ship::Ship) = ship.position
waypoint(ship::Ship) = ship.waypoint

function rotate(facing::Symbol, direction::Symbol, deg::Integer=90)
    steps, reminder = divrem(deg, 90)
    reminder == 0 || error("ups, it can move diagonally!")
    direction = direction == :Right ? +1 : -1
    direction_index = findfirst(directions .== facing)
    direction_index = mod1(direction_index + steps * direction, 4)
    return directions[direction_index]
end

@assert rotate(:East, :Left) == :North
@assert rotate(:North, :Left) == :West
@assert rotate(:North, :Right) == :East
@assert rotate(:East, :Right) == :South
@assert rotate(:North, :Left, 360) == :North

function part1(string::AbstractString; verbose=false)
    ship = Ship(:East, (0, 0), (0, 0))

    for row in split(string, '\n', keepempty=false)
        row = strip(row)
        m = match(r"([NSEWLRF])(\d+)", row)
        if isnothing(m)
            error("invalid input")
        end
        action, value = m.captures
        value = parse(Int, value)

        if action in ("L", "R")
            direction = rotate(facing(ship), action == "L" ? :Left : :Right, value)
            ship.facing = direction
            value = 0
        elseif action == "F"
            direction = facing(ship)
        else
            direction = char_to_direction[action]
        end

        x, y = position(ship)
        if direction == :North
            x += value
        elseif direction == :South
            x -= value
        elseif direction == :East
            y += value
        else
            y -= value
        end
        ship.position = (x, y)

    end
    return position(ship) .|> abs |> sum
end

@assert part1(example) == 25

function rotate(point::Point, deg::Integer)::Point
    θ = deg2rad(deg)
    x, y = point
    x_new = cos(θ) * x - sin(θ) * y
    y_new = sin(θ) * x + cos(θ) * y
    return Int(round(x_new)), Int(round(y_new))
end

@assert rotate((10, 4), 90) == (-4, 10)

function part2(string::AbstractString; verbose=false)
    ship = Ship(:East, (0, 0), (1, 10))

    for row in split(string, '\n', keepempty=false)
        row = strip(row)
        m = match(r"([NSEWLRF])(\d+)", row)
        if isnothing(m)
            error("invalid input")
        end
        action, value = m.captures
        value = parse(Int, value)

        if action in ("N", "S", "W", "E")
            x, y = waypoint(ship)
            if action == "N"
                x += value
            elseif action == "S"
                x -= value
            elseif action == "E"
                y += value
            else
                y -= value
            end
            ship.waypoint = (x, y)
        elseif action in ("L", "R")
            if action == "L"
                value = 360 - value
            end
            ship.waypoint = rotate(waypoint(ship), value)
        else
            wx, wy = waypoint(ship)
            px, py = position(ship)
            ship.position = (px + wx * value, py + wy * value)
        end

        verbose && println("$(row):\t Ship at $(position(ship)), with waypoint $(waypoint(ship))")

    end
    return position(ship) .|> abs |> sum
end

@assert part2(example) == 286

test = read("data/day-12.txt", String)
println("Part 1: $(result1 = part1(test))")
println("Part 2: $(result2 = part2(test))")

@assert result1 == 415
@assert result2 == 29401

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
#
# Correctness Tests
#

using Test
using PowerAnalysis
using HypothesisTests

include("t_test.jl")
include("z_test.jl")


text:
using Test
using MianalyzerBackend
const Stats = MianalyzerBackend.Stats

@testset "Stats" begin
  @testset ".max" begin
    @test Stats.max([1 2 3 4 5]) == 5
    @test Stats.max([10.1 10.2 10.9 10.6]) == 10.9
  end

  @testset ".min" begin
    @test Stats.min([1 2 3 4 5]) == 1
    @test Stats.min([10.1 10.2 10.9 10.6]) == 10.1
  end

  @testset ".mean" begin
    @test Stats.mean([1 2 3 4 5]) == 3
    @test Stats.mean([10.1 10.2 10.9 10.6]) == 10.45
  end

  @testset ".sum" begin
    @test Stats.sum([1 2 3 4 5]) == 15
    @test Stats.sum([10.1 10.2 10.9 10.6]) == 41.8
  end

  @testset ".count" begin
    @test Stats.count([1 2 3 4 5]) == 5
    @test Stats.count([10.1 10.2 10.9 10.6]) == 4
  end

end


text:
include("debug.jl")

# Data Structures

mutable struct SSLConfig
    data::Ptr{Cvoid}
    rng
    chain::CRT
    dbg
    cert
    key
    alpn_protos

    function SSLConfig()
        conf = new()
        conf.data = Libc.malloc(1000)  # 360
        ccall((:mbedtls_ssl_config_init, libmbedtls), Cvoid, (Ptr{Cvoid},), conf.data)
        finalizer(x->begin
            data = x.data
            @async begin
                ccall((:mbedtls_ssl_config_free, libmbedtls),
                      Cvoid, (Ptr{Cvoid},), data)
                Libc.free(data)
            end
        end, conf)
        conf
    end
end

Base.show(io::IO, c::SSLConfig) = print(io, "MbedTLS.SSLConfig()")

mutable struct SSLContext <: IO
    data::Ptr{Cvoid}
    datalock::ReentrantLock
    waitlock::ReentrantLock
    config::SSLConfig
    isreadable::Bool
    bytesavailable::Int
    close_notify_sent::Bool
    async_exception::Union{Nothing,MbedException}
    bio

    function SSLContext()
        ctx = new()
        ctx.data = Libc.malloc(1000)  # 488
        ctx.datalock = ReentrantLock()
        ctx.waitlock = ReentrantLock()
        ctx.isreadable = false
        ctx.bytesavailable = -1
        ctx.close_notify_sent = false
        ctx.async_exception = nothing
        ccall((:mbedtls_ssl_init, libmbedtls), Cvoid, (Ptr{Cvoid},), ctx.data)
        finalizer(x->begin
            data = x.data
            @async begin
                ccall((:mbedtls_ssl_free, libmbedtls),
                      Cvoid, (Ptr{Cvoid},), ctx.data)
                Libc.free(ctx.data)
            end
        end, ctx)
        ctx
    end
end

# Handshake

function handshake(ctx::SSLContext)

    ctx.isreadable && throw(ArgumentError("handshake() already done!"))
                                                                                ;@😬 "🤝 ..."
    while true
        n = ssl_handshake(ctx)
        if n == 0
            break
        elseif n == MBEDTLS_ERR_SSL_WANT_READ                                   ;@😬 "🤝  ⌛"
            if eof(ctx.bio)
                throw(EOFError())                                               ;@💀 "🤝  🛑"
            end
        else
            ssl_abandon(ctx)                                                    ;@💀 "🤝  💥"
            throw(MbedException(n))
        end
    end
                                                                                ;@😬 "🤝  ✅"
    ctx.isreadable = true
    ctx.bytesavailable = 0
    ctx.close_notify_sent = false

    @async try
        while ctx.isreadable
            wait_for_decrypted_data(ctx)
            while ctx.bytesavailable > 0
                sleep(5)
            end
        end
    catch e
        ctx.async_exception = e
    end

    nothing
end

function check_async_exception(ctx::SSLContext)
    if ctx.async_exception !== nothing
        e = ctx.async_exception
        ctx.async_exception = nothing
        throw(e)
    end
end

# Fatal Errors

"""
The documentation for `ssl_read`, `ssl_write` and `ssl_close_notify` all say:

> If this function returns something other than 0 or
> MBEDTLS_ERR_SSL_WANT_READ/WRITE, you must stop using the SSL context
> for reading or writing, and either free it or call

This function ensures that the `SSLContext` is won't be used again.
"""
function ssl_abandon(ctx::SSLContext)                                           ;@💀 "ssl_abandon 💥"
    ctx.isreadable = false
    ctx.bytesavailable = 0
    ctx.close_notify_sent = true
    close(ctx.bio)
    n = ssl_session_reset(ctx)
    n == 0 || throw(MbedException(n))
end

# Base ::IO Connection State Methods

Sockets.getsockname(ctx::SSLContext) = Sockets.getsockname(ctx.bio)

"""
    isreadable(ctx::SSLContext)

True unless:
 - TLS `close_notify` was received, or
 - the peer closed the connection (and the TLS buffer is empty), or
 - an un-handled exception occurred while reading.
"""
function Base.isreadable(ctx::SSLContext)
    ctx.isreadable || return false
    # It's possible we received the shutdown, but didn't process it yet - if so,
    # do that now.
    if bytesavailable(ctx.bio) > 0 || ssl_check_pending(ctx)
        ssl_unsafe_read(ctx, Ptr{UInt8}(C_NULL), UInt(0))
    end
    return ctx.isreadable
end

"""
    iswritable(ctx::SSLContext)

True unless:
 - `close(::SSLContext)` is called, or
 -  the peer closed the connection.
"""
Base.iswritable(ctx::SSLContext) = !ctx.close_notify_sent && isopen(ctx.bio)

"""
    isopen(ctx::SSLContext)

Same as `iswritable(ctx)`.
> "...a closed stream may still have data to read in its buffer,
>  use eof to check for the ability to read data." [?Base.isopen]
"""
Base.isopen(ctx::SSLContext) = iswritable(ctx)

@static if isdefined(Base, :bytesavailable)
"""
    bytesavailable(ctx::SSLContext)

Number of decrypted bytes waiting in the TLS buffer.
"""
Base.bytesavailable(ctx::SSLContext) = ctx.bytesavailable
else
Base.nb_available(ctx::SSLContext) = ctx.bytesavailable
end

"""
    eof(ctx::SSLContext)

True if not `isreadable` and there are no more `bytesavailable` to read.
"""
function Base.eof(ctx::SSLContext)
    if ctx.bytesavailable > 0
        return false
    end                                                                         ;@😬 "eof ⌛"
    wait_for_decrypted_data(ctx)
    @assert ctx.bytesavailable > 0 || !ctx.isreadable
    return ctx.bytesavailable <= 0
end

"""
    close(ctx::SSLContext)

Send a TLS `close_notify` message to the peer.
"""
function Base.close(ctx::SSLContext)                                            ;@💀 "close iswritable=$(iswritable(ctx))"

    if iswritable(ctx)

        n = ssl_close_notify(ctx)
        ctx.close_notify_sent = true                                            ;@💀 "close 🗣"

        if n == MBEDTLS_ERR_SSL_WANT_READ || n == MBEDTLS_ERR_SSL_WANT_WRITE    ;@💀 "close ⌛"
            @assert false "Should not get to here because `f_send` " *
                          "never returns ...WANT_READ/WRITE."
        elseif n != 0
            ssl_abandon(ctx)
            throw(MbedException(n))
        end
    end
    @assert !iswritable(ctx)
    nothing
end

# Sending Data

"""
Copy `nbytes` from `buf` to `ctx`.

The TLS library function `ssl_write` is called as many times as needed to send
all the data. The TLS library encrypts the data and passes it to the `f_send`
function which sends it to the underlying connection (`ctx.bio`).
See `f_send` and `set_bio!` below.
"""
function ssl_unsafe_write(ctx::SSLContext, buf::Ptr{UInt8}, nbytes::UInt)

    Base.check_open(ctx.bio) # Throw error if `bio` is closed.

    iswritable(ctx) ||
    throw(ArgumentError("`unsafe_write` requires `iswritable(::SSLContext)`"))

    nwritten = 0                                                                ;@🤖 "ssl_write ➡  $nbytes"
    while nwritten < nbytes
        n = ssl_write(ctx, buf + nwritten, nbytes - nwritten)
        if n == MBEDTLS_ERR_SSL_WANT_READ || n == MBEDTLS_ERR_SSL_WANT_WRITE    ;@💀 "ssl_write ⌛"
            @assert false "Should not get to here because `f_send` " *
                          "never returns ...WANT_READ/WRITE."
            yield()
            continue
        elseif n == MBEDTLS_ERR_NET_CONN_RESET
            ssl_abandon(ctx)                                                    ;@🤖 "ssl_write 🛑"
            Base.check_open(ctx.bio)
            @assert false
        elseif n < 0
            ssl_abandon(ctx)                                                    ;@🤖 "ssl_write 💥"
            throw(MbedException(n))
        end
        nwritten += n
    end
    return Int(nwritten)
end

# Sending Encrypted Data

"""
Copy `nbytes` of encrypted data from `buf` to the underlying `bio` connection.
"""
function f_send(c_bio, buf, nbytes)                                             ;@🤖 "f_send ➡  $nbytes"
    bio = unsafe_pointer_to_objref(c_bio)
    if !isopen(bio) || bio.status == Base.StatusClosing
        return Cint(MBEDTLS_ERR_NET_CONN_RESET)
    end
    return Cint(unsafe_write(bio, buf, nbytes))
end

"""
Connect `f_send` and `f_recv` callback functions to `SSLContext`.
"""
function set_bio!(ctx::SSLContext, bio::T) where {T<:IO}
    isopen(bio) || throw(ArgumentError("`set_bio!` requires `isopen(bio)`"))
    ctx.bio = bio
    ssl_set_bio(ctx, pointer_from_objref(bio), c_send[], c_recv[])
    nothing
end

# Receiving Data

"""
While there are no decrypted bytes available but the connection is readable:
 - If the TLS buffer has no pending (unprocessed) data wait for
   more encrypted data to arrive on the underlying connection.
 - Run a zero-byte read to allow the library to process its internal buffer,
   and/or read from the underlying connection.
 - `ssl_unsafe_read` updates the `isreadable` and `bytesavailable` state.
"""
function wait_for_decrypted_data(ctx)
    lock(ctx.waitlock)
    try
        check_async_exception(ctx)
        while ctx.isreadable && ctx.bytesavailable <= 0
            if !ssl_check_pending(ctx)                                          ;@🤖 "wait_for_encrypted_data ⌛";
                wait_for_encrypted_data(ctx)
            end
            ssl_unsafe_read(ctx, Ptr{UInt8}(C_NULL), UInt(0))                   ;@🤖 "wait_for_decrypted_data 📥  $(ctx.bytesavailable)"
        end
    finally
        unlock(ctx.waitlock)
    end
end

"""
    ssl_unsafe_read(::SSLContext, buf, nbytes)

Copies at most `nbytes` of decrypted data into `buf`.
Never blocks to wait for more data to arrive.
Returns number of bytes copied into `buf` (`<= nbytes`).
Updates `ctx.bytesavailable` with the number of decrypted bytes remaining in
the TLS buffer.

Stops when:
 - `nbytes` have been copied, or
 - there are no more decrypted bytes available in the TLS buffer, or
 - a TLS `close_notify` message is received.

When TLS `close_notify` is received:
 - `isreadable` is set to false and `bytesavailable` is set to zero.
   [RFC5246 7.2.1]: "Any data received after a closure alert is ignored."
 - the number of bytes read before the `close_notify` is returned as usual.

Throws a `MbedException` if `ssl_read` returns an unhandled error code.

When an unhandled exception occurs `isreadable` is set to false.
"""
function ssl_unsafe_read(ctx::SSLContext, buf::Ptr{UInt8}, nbytes::UInt)

    check_async_exception(ctx)

    ctx.isreadable ||
    throw(ArgumentError("`ssl_unsafe_read` requires `isreadable(::SSLContext)`"))

    nread::UInt = 0
    try
        while true

            n = ssl_read(ctx, buf + nread, nbytes - nread)                      ;@😬 "ssl_read ⬅  $n $(n == MBEDTLS_ERR_SSL_PEER_CLOSE_NOTIFY ? "(CLOSE_NOTIFY)" :
                                                                                n == MBEDTLS_ERR_NET_CONN_RESET        ? "(CONN_RESET)" :
                                                                                n == MBEDTLS_ERR_SSL_WANT_READ         ? "(WANT_READ)" : "")"
            if n == MBEDTLS_ERR_SSL_PEER_CLOSE_NOTIFY ||
               n == MBEDTLS_ERR_NET_CONN_RESET
                ssl_abandon(ctx)
                return nread
            elseif n == MBEDTLS_ERR_SSL_WANT_READ
                ctx.bytesavailable = 0                                          ;@😬 "ssl_read ⌛ $nread"
                return nread
            elseif n < 0
                ssl_abandon(ctx)
                throw(MbedException(n))
            end

            nread += n
            @assert nread <= nbytes

            if nread == nbytes
                ctx.bytesavailable = ssl_get_bytes_avail(ctx)                   ;@🤖 "ssl_read ⬅  $nread, 📥  $(ctx.bytesavailable)"
                return nread
            end
        end
    catch e                                                                     ;@💀 "ssl_read 💥"
        ssl_abandon(ctx)
        rethrow(e)
    end

    @assert false "unreachable"
end

# Receiving Encrypted Data

function wait_for_encrypted_data(ctx)
    try
        eof(ctx.bio)
    catch e
        if !(e isa Base.IOError) || e.code != Base.UV_ECONNRESET
            rethrow(e)
        end
    end
end

"""
Copy at most `nbytes` of encrypted data to `buf` from the `bio` connection.
If no encrypted bytes are available return:
 - `MBEDTLS_ERR_SSL_WANT_READ` if the connection is still open, or
 - `MBEDTLS_ERR_NET_CONN_RESET` if it is closed.
"""
function f_recv(c_bio, buf, nbytes)
    @assert nbytes > 0
    bio = unsafe_pointer_to_objref(c_bio)
    n = bytesavailable(bio)
    if n == 0                                                                   ;@🤖 "f_recv $(isopen(bio) ? "WANT_READ" : "CONN_RESET")"
        return isopen(bio) ? Cint(MBEDTLS_ERR_SSL_WANT_READ) :
                             Cint(MBEDTLS_ERR_NET_CONN_RESET)
    end
    n = min(nbytes, n)                                                          ;@🤖 "f_recv ⬅  $n"
    unsafe_read(bio, buf, n)
    return Cint(n)
end

# Base ::IO Write Methods -- wrappers for `ssl_unsafe_write`

Base.unsafe_write(ctx::SSLContext, msg::Ptr{UInt8}, N::UInt) =
    ssl_unsafe_write(ctx, msg, N)

Base.write(ctx::SSLContext, msg::UInt8) = write(ctx, Ref(msg))

# Base ::IO Read Methods -- wrappers for `ssl_unsafe_read`

"""
    unsafe_read(ctx::SSLContext, buf::Ptr{UInt8}, nbytes::UInt)

Copy `nbytes` of decrypted data from `ctx` into `buf`.
Wait for sufficient decrypted data to be available.
Throw `EOFError` if the peer sends TLS `close_notify` or closes the
connection before `nbytes` have been copied.
"""
function Base.unsafe_read(ctx::SSLContext, buf::Ptr{UInt8}, nbytes::UInt)
    nread = 0
    while nread < nbytes
        if eof(ctx)                                                             ;@💀 "unsafe_read 🛑"
            throw(EOFError())
        end
        nread += ssl_unsafe_read(ctx, buf + nread, nbytes - nread)
    end                                                                         ;@😬 "unsafe_read ⬅ $nread"
    nothing
end

"""
    readbytes!(ctx::SSLContext, buf::Vector{UInt8}, nbytes=length(buf); kw...)

Copy at most `nbytes` of decrypted data from `ctx` into `buf`.
If `all=true`: wait for sufficient decrypted data to be available.
Less than `nbytes` may be copied if the peer sends TLS `close_notify` or closes
the connection.
Returns number of bytes copied into `buf` (`<= nbytes`).
"""
Base.readbytes!(ctx::SSLContext, buf::Vector{UInt8}, nbytes=length(buf); kw...) =
    readbytes!(ctx, buf, UInt(nbytes); kw...)

function Base.readbytes!(ctx::SSLContext, buf::Vector{UInt8}, nbytes::UInt;
                         all::Bool=true)
    nbytes <= length(buf) || throw(ArgumentError("`buf` too small!"))
    nread = 0
    while nread < nbytes
        nread += ssl_unsafe_read(ctx, pointer(buf) + nread, nbytes - nread)
        if (nread == nbytes) || !all || eof(ctx)
            break
        end
    end                                                                         ;@😬 "readbytes! ⬅  $nread"
    return nread
end

"""
    readavailable(ctx::SSLContext)

Read available decrypted data from `ctx`,
but don't wait for more data to arrive.

The amount of decrypted data that can be read at once is limited by
`MBEDTLS_SSL_MAX_CONTENT_LEN`.
"""
function Base.readavailable(ctx::SSLContext)
    n = UInt(MBEDTLS_SSL_MAX_CONTENT_LEN)
    buf = Vector{UInt8}(undef, n)
    n = ssl_unsafe_read(ctx, pointer(buf), n)                                   ;@😬 "readavailable ⬅  $n"
    return resize!(buf, n)
end

# Configuration

function config_defaults!(config::SSLConfig; endpoint=MBEDTLS_SSL_IS_CLIENT,
    transport=MBEDTLS_SSL_TRANSPORT_STREAM, preset=MBEDTLS_SSL_PRESET_DEFAULT)
    @err_check ccall((:mbedtls_ssl_config_defaults, libmbedtls), Cint,
        (Ptr{Cvoid}, Cint, Cint, Cint),
        config.data, endpoint, transport, preset)
end

function authmode!(config::SSLConfig, auth)
    ccall((:mbedtls_ssl_conf_authmode, libmbedtls), Cvoid,
        (Ptr{Cvoid}, Cint),
        config.data, auth)
end

function rng!(config::SSLConfig, f_rng::Ptr{Cvoid}, rng)
    ccall((:mbedtls_ssl_conf_rng, libmbedtls), Cvoid,
        (Ptr{Cvoid}, Ptr{Cvoid}, Any),
        config.data, f_rng, rng)
end

function rng!(config::SSLConfig, rng::AbstractRNG)
    config.rng = rng
    rng!(config, c_rng[], rng)
end

function ca_chain!(config::SSLConfig, chain=crt_parse(DEFAULT_CERT))
    config.chain = chain
    ccall((:mbedtls_ssl_conf_ca_chain, libmbedtls), Cvoid,
        (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        config.data, chain.data, C_NULL)
end

"""
Enable / Disable renegotiation support for connection when initiated by peer
(MBEDTLS_SSL_RENEGOTIATION_ENABLED or MBEDTLS_SSL_RENEGOTIATION_DISABLED).
See: https://tls.mbed.org/api/ssl_8h.html#aad4f50fc1c0a018fd5eb18fd9621d0d3
"""
function ssl_conf_renegotiation!(config::SSLConfig, renegotiation)
    ccall((:mbedtls_ssl_conf_renegotiation, libmbedtls), Cvoid,
        (Ptr{Cvoid}, Cint),
        config.data, renegotiation)
end

function own_cert!(config::SSLConfig, cert::CRT, key::PKContext)
    config.cert = cert
    config.key = key
    @err_check ccall((:mbedtls_ssl_conf_own_cert, libmbedtls), Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        config.data, cert.data, key.data)
end

function setup!(ctx::SSLContext, conf::SSLConfig)
    ctx.config = conf
    ssl_setup(ctx, conf)
end

function dbg!(conf::SSLConfig, f::Ptr{Cvoid}, p)
    ccall((:mbedtls_ssl_conf_dbg, libmbedtls), Cvoid,
        (Ptr{Cvoid}, Ptr{Cvoid}, Any),
        conf.data, f, p)
end

function f_dbg(f, level, filename, number, msg)
    f(level, unsafe_string(filename), number, unsafe_string(msg))
    nothing
end

function dbg!(conf::SSLConfig, f)
    conf.dbg = f
    dbg!(conf, c_dbg[], f)
    nothing
end

@enum(DebugThreshold,
    NONE = 0,
    ERROR,
    STATE_CHANGE,
    INFO,
    VERBOSE)

function set_dbg_level(level)
    ccall((:mbedtls_debug_set_threshold, libmbedtls), Cvoid,
        (Cint,), Cint(level))
    nothing
end

function set_alpn!(conf::SSLConfig, protos)
    conf.alpn_protos = protos
    @err_check ccall((:mbedtls_ssl_conf_alpn_protocols, libmbedtls), Cint,
                     (Ptr{Cvoid}, Ptr{Ptr{Cchar}}), conf.data, protos)
    nothing
end

function alpn_proto(ctx::SSLContext)
    rv = ccall((:mbedtls_ssl_get_alpn_protocol, libmbedtls), Ptr{Cchar},
               (Ptr{Cvoid},), ctx.data)
    unsafe_string(rv)
end

# C API

macro lockdata(ctx, expr)
    esc(quote
        lock($ctx.datalock)
        @assert $ctx.datalock.reentrancy_cnt == 1
        try
            $expr
        finally
            unlock($ctx.datalock)
        end
    end)
end

function get_peer_cert(ctx::SSLContext)
    data = ccall((:mbedtls_ssl_get_peer_cert, libmbedtls), Ptr{Cvoid}, (Ptr{Cvoid},), ctx.data)
    return CRT(data)
end

function get_version(ctx::SSLContext)
    if isdefined(ctx, :config)
        data = ccall((:mbedtls_ssl_get_version, libmbedtls), Ptr{UInt8}, (Ptr{Cvoid},), ctx.data)
        return unsafe_string(data)
    else
        throw(ArgumentError("`ctx` hasn't been initialized with an MbedTLS.SSLConfig; run `MbedTLS.setup!(ctx, conf)`"))
    end
end

function get_ciphersuite(ctx::SSLContext)
    data = ccall((:mbedtls_ssl_get_ciphersuite, libmbedtls), Ptr{UInt8}, (Ptr{Cvoid},), ctx.data)
    return unsafe_string(data)
end

function ssl_setup(ctx::SSLContext, conf::SSLConfig)
    @lockdata ctx begin
        @err_check ccall((:mbedtls_ssl_setup, libmbedtls), Cint,
            (Ptr{Cvoid}, Ptr{Cvoid}),
            ctx.data, conf.data)
    end
end

function ssl_handshake(ctx::SSLContext)
    n = @lockdata ctx begin
        ccall((:mbedtls_ssl_handshake, libmbedtls), Cint,
              (Ptr{Cvoid},), ctx.data)
    end
end

"""
Return the number of application data bytes remaining to be read
from the current record.

https://tls.mbed.org/api/ssl_8h.html#ad43142085f3182e9b0dc967ec582032b:
"""
function ssl_get_bytes_avail(ctx::SSLContext)::Int
    @lockdata ctx begin
        return ccall((:mbedtls_ssl_get_bytes_avail, libmbedtls),
                     Csize_t, (Ptr{Cvoid},), ctx.data)
    end
end

"""
    ssl_check_pending(::SSLContext)::Bool

Check if there is data already read from the underlying transport
but not yet processed.

If the SSL/TLS module successfully returns from an operation - e.g.
a handshake or an application record read - and you're awaiting
incoming data next, you must not immediately idle on the underlying
transport to have data ready, but you need to check the value of
this function first.  The reason is that the desired data might
already be read but not yet processed.  If, in contrast, a previous
call to the SSL/TLS module returned MBEDTLS_ERR_SSL_WANT_READ, it
is not necessary to call this function, as the latter error code
entails that all internal data has been processed.

https://tls.mbed.org/api/ssl_8h.html#a4075f7de9877fd667bcfa2e819e33426
"""
function ssl_check_pending(ctx::SSLContext)::Bool
    @lockdata ctx begin
	return ccall((:mbedtls_ssl_check_pending, libmbedtls),
		     Cint, (Ptr{Cvoid},), ctx.data) > 0
    end
end

function ssl_set_bio(ctx::SSLContext, bio, f_send::Ptr{Cvoid}, f_recv::Ptr{Cvoid})
    @lockdata ctx begin
        ccall((:mbedtls_ssl_set_bio, libmbedtls), Cvoid,
            (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
            ctx.data, bio, f_send, f_recv, C_NULL)
    end
end

"""
    ssl_read(::SSLContext, ptr, n)

Returns One of the following:
0 if the read end of the underlying transport was closed,
the (positive) number of bytes read, or
a negative error code on failure.

If MBEDTLS_ERR_SSL_WANT_READ is returned, no application data is
available from the underlying transport. In this case, the function
needs to be called again at some later stage.

If this function returns something other than a positive value or
MBEDTLS_ERR_SSL_WANT_READ/WRITE or MBEDTLS_ERR_SSL_CLIENT_RECONNECT,
you must stop using the SSL context for reading or writing, and
either free it or call mbedtls_ssl_session_reset() on it before
re-using it for a new connection; the current connection must be
closed.

https://tls.mbed.org/api/ssl_8h.html#aa2c29eeb1deaf5ad9f01a7515006ede5
"""
function ssl_read(ctx::SSLContext, ptr, n)::Int
    @lockdata ctx begin
        return ccall((:mbedtls_ssl_read, libmbedtls), Cint,
                     (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
                     ctx.data, ptr, n)
    end
end

"""
    ssl_write(::SSLContext, ptr, n)

When this function returns MBEDTLS_ERR_SSL_WANT_WRITE/READ, it must
be called later with the same arguments, until it returns a value
greater that or equal to 0. When the function returns
MBEDTLS_ERR_SSL_WANT_WRITE there may be some partial data in the
output buffer, however this is not yet sent.

If this function returns something other than 0, a positive value
or MBEDTLS_ERR_SSL_WANT_READ/WRITE, you must stop using the SSL
context for reading or writing, and either free it or call
mbedtls_ssl_session_reset() on it before re-using it for a new
connection; the current connection must be closed.

https://tls.mbed.org/api/ssl_8h.html#a5bbda87d484de82df730758b475f32e5
"""
function ssl_write(ctx::SSLContext, ptr, n)::Int
    @lockdata ctx begin
        return ccall((:mbedtls_ssl_write, libmbedtls), Cint,
                     (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
                     ctx.data, ptr, n)
    end
end

"""
    ssl_close_notify(::SSLContext)

Notify the peer that the connection is being closed.
0 if successful, or a specific SSL error code.

If this function returns something other than 0 or
MBEDTLS_ERR_SSL_WANT_READ/WRITE, you must stop using the SSL context
for reading or writing, and either free it or call
mbedtls_ssl_session_reset() on it before re-using it for a new
connection; the current connection must be closed.

https://tls.mbed.org/api/ssl_8h.html#ac2c1b17128ead2df3082e27b603deb4c
"""
function ssl_close_notify(ctx::SSLContext)
    @lockdata ctx begin
        return ccall((:mbedtls_ssl_close_notify, libmbedtls),
                     Cint, (Ptr{Cvoid},), ctx.data)
    end
end

"""
Reset an already initialized SSL context for re-use while retaining
application-set variables, function pointers and data.
"""
function ssl_session_reset(ctx::SSLContext)
    @lockdata ctx begin
        return ccall((:mbedtls_ssl_session_reset, libmbedtls),
                     Cint, (Ptr{Cvoid},), ctx.data)
    end
end

function hostname!(ctx::SSLContext, hostname)
    @err_check ccall((:mbedtls_ssl_set_hostname, libmbedtls), Cint,
      (Ptr{Cvoid}, Cstring), ctx.data, hostname)
end

const c_send = Ref{Ptr{Cvoid}}(C_NULL)
const c_recv = Ref{Ptr{Cvoid}}(C_NULL)
const c_dbg = Ref{Ptr{Cvoid}}(C_NULL)
const DEFAULT_CERT = read(joinpath(@__DIR__, "cacert.pem"), String)

function __sslinit__()
    c_send[] = @cfunction(f_send, Cint, (Ptr{Cvoid}, Ptr{UInt8}, Csize_t))
    c_recv[] = @cfunction(f_recv, Cint, (Ptr{Cvoid}, Ptr{UInt8}, Csize_t))
    c_dbg[] = @cfunction(f_dbg, Cvoid, (Any, Cint, Ptr{UInt8}, Cint, Ptr{UInt8}))
end