text:
#######################################################################
# Hamiltonian
#######################################################################
struct HamiltonianHarmonic{L,M,A<:Union{AbstractMatrix{M},SparseMatrixBuilder{M}}}
    dn::SVector{L,Int}
    h::A
end

HamiltonianHarmonic{L,M,A}(dn::SVector{L,Int}, n::Int, m::Int) where {L,M,A<:SparseMatrixCSC{M}} =
    HamiltonianHarmonic(dn, sparse(Int[], Int[], M[], n, m))

HamiltonianHarmonic{L,M,A}(dn::SVector{L,Int}, n::Int, m::Int) where {L,M,A<:Matrix{M}} =
    HamiltonianHarmonic(dn, zeros(M, n, m))

struct Hamiltonian{LA<:AbstractLattice,L,M,A<:AbstractMatrix,
                   H<:HamiltonianHarmonic{L,M,A},
                   O<:Tuple{Vararg{Tuple{Vararg{NameType}}}}} <: AbstractMatrix{M}
    lattice::LA
    harmonics::Vector{H}
    orbitals::O
end

function Hamiltonian(lat, hs::Vector{H}, orbs, n::Int, m::Int) where {L,M,H<:HamiltonianHarmonic{L,M}}
    sort!(hs, by = h -> abs.(h.dn))
    if isempty(hs) || !iszero(first(hs).dn)
        pushfirst!(hs, H(zero(SVector{L,Int}), empty_sparse(M, n, m)))
    end
    return Hamiltonian(lat, hs, orbs)
end

Base.show(io::IO, ham::Hamiltonian) = show(io, MIME("text/plain"), ham)
function Base.show(io::IO, ::MIME"text/plain", ham::Hamiltonian)
    i = get(io, :indent, "")
    print(io, i, summary(ham), "\n",
"$i  Bloch harmonics  : $(length(ham.harmonics)) ($(displaymatrixtype(ham)))
$i  Harmonic size    : $((n -> "$n × $n")(nsites(ham)))
$i  Orbitals         : $(displayorbitals(ham))
$i  Element type     : $(displayelements(ham))
$i  Onsites          : $(nonsites(ham))
$i  Hoppings         : $(nhoppings(ham))
$i  Coordination     : $(nhoppings(ham) / nsites(ham))")
    ioindent = IOContext(io, :indent => string("  "))
    issuperlattice(ham.lattice) && print(ioindent, "\n", ham.lattice.supercell)
end

Base.summary(h::Hamiltonian{LA}) where {E,L,LA<:Lattice{E,L}} =
    "Hamiltonian{<:Lattice} : Hamiltonian on a $(L)D Lattice in $(E)D space"

Base.summary(::Hamiltonian{LA}) where {E,L,T,L´,LA<:Superlattice{E,L,T,L´}} =
    "Hamiltonian{<:Superlattice} : $(L)D Hamiltonian on a $(L´)D Superlattice in $(E)D space"

Base.eltype(::Hamiltonian{<:Any,<:Any,M}) where {M} = M

displaymatrixtype(h::Hamiltonian) = displaymatrixtype(matrixtype(h))
displaymatrixtype(::Type{<:SparseMatrixCSC}) = "SparseMatrixCSC, sparse"
displaymatrixtype(::Type{<:Array}) = "Matrix, dense"
displaymatrixtype(A::Type{<:AbstractArray}) = string(A)
displayelements(h::Hamiltonian) = displayelements(blocktype(h))
displayelements(::Type{S}) where {N,T,S<:SMatrix{N,N,T}} = "$N × $N blocks ($T)"
displayelements(::Type{T}) where {T} = "scalar ($T)"
displayorbitals(h::Hamiltonian) =
    replace(replace(string(h.orbitals), "Symbol(\"" => ":"), "\")" => "")

SparseArrays.issparse(h::Hamiltonian{LA,L,M,A}) where {LA,L,M,A<:AbstractSparseMatrix} = true
SparseArrays.issparse(h::Hamiltonian{LA,L,M,A}) where {LA,L,M,A} = false

# Internal API #

latdim(h::Hamiltonian{LA}) where {E,L,LA<:AbstractLattice{E,L}} = L

matrixtype(::Hamiltonian{LA,L,M,A}) where {LA,L,M,A} = A
realtype(::Hamiltonian{<:Any,<:Any,M}) where {M} = real(eltype(M))

# find SVector type that can hold all orbital amplitudes in any lattice sites
orbitaltype(orbs, type::Type{Tv} = Complex{T}) where {T,Tv} =
    _orbitaltype(SVector{1,Tv}, orbs...)
_orbitaltype(::Type{S}, ::NTuple{D,NameType}, os...) where {N,Tv,D,S<:SVector{N,Tv}} =
    (M = max(N,D); _orbitaltype(SVector{M,Tv}, os...))
_orbitaltype(t::Type{SVector{N,Tv}}) where {N,Tv} = t
_orbitaltype(t::Type{SVector{1,Tv}}) where {Tv} = Tv

# find SMatrix type that can hold all matrix elements between lattice sites
blocktype(orbs, type::Type{Tv}) where {Tv} =
    _blocktype(orbitaltype(orbs, Tv))
_blocktype(::Type{S}) where {N,Tv,S<:SVector{N,Tv}} = SMatrix{N,N,Tv,N*N}
_blocktype(::Type{S}) where {S<:Number} = S

blocktype(h::Hamiltonian{LA,L,M}) where {LA,L,M} = M

promote_blocktype(hs::Hamiltonian...) = promote_blocktype(blocktype.(hs)...)
promote_blocktype(s1::Type, s2::Type, ss::Type...) =
    promote_blocktype(promote_blocktype(s1, s2), ss...)
promote_blocktype(::Type{SMatrix{N1,N1,T1,NN1}}, ::Type{SMatrix{N2,N2,T2,NN2}}) where {N1,NN1,T1,N2,NN2,T2} =
    SMatrix{max(N1, N2), max(N1, N2), promote_type(T1, T2), max(NN1,NN2)}
promote_blocktype(T1::Type{<:Number}, T2::Type{<:Number}) = promote_type(T1, T2)
promote_blocktype(T::Type) = T

blockdim(h::Hamiltonian) = blockdim(blocktype(h))
blockdim(::Type{S}) where {N,S<:SMatrix{N,N}} = N
blockdim(::Type{T}) where {T<:Number} = 1

checkfinitedim(h::Hamiltonian{LA,L}) where {LA,L} =
    L == 0 && throw(ArgumentError("A finite-dimensional Hamiltonian is required, not zero-dimensional"))

function nhoppings(ham::Hamiltonian)
    count = 0
    for h in ham.harmonics
        count += iszero(h.dn) ? (_nnz(h.h) - _nnzdiag(h.h)) : _nnz(h.h)
    end
    return count
end

function nonsites(ham::Hamiltonian)
    count = 0
    for h in ham.harmonics
        iszero(h.dn) && (count += _nnzdiag(h.h))
    end
    return count
end

_nnz(h::AbstractSparseMatrix) = nnz(h)
_nnz(h::DenseMatrix) = count(!iszero, h)

function _nnzdiag(s::SparseMatrixCSC)
    count = 0
    rowptrs = rowvals(s)
    for col in 1:size(s,2)
        for ptr in nzrange(s, col)
            rowptrs[ptr] == col && (count += 1; break)
        end
    end
    return count
end
_nnzdiag(s::Matrix) = count(!iszero, s[i,i] for i in 1:minimum(size(s)))

# Iteration tools #

struct EachIndexNonzeros{H}
    h::H
    rowrange::UnitRange{Int}
    colrange::UnitRange{Int}
end

eachindex_nz(h, rowrange = 1:size(h, 1), colrange = 1:size(h, 2)) =
    EachIndexNonzeros(h, rclamp(rowrange, 1:size(h, 1)), rclamp(colrange, 1:size(h, 2)))

function firststate(itr::EachIndexNonzeros{<:Hamiltonian}, nhar)
    m = itr.h.harmonics[nhar].h
    row, col = nextnonzero_row_col(m, itr)
    return (row, col, nhar)
end

function nextnonzero_row_col(m::DenseMatrix, itr, col = first(itr.colrange))
    for col´ in col:last(itr.colrange), row in itr.rowrange
        iszero(m[row, col´]) || return (row, col´)
    end
    # (0, 0) is sentinel for "no non-zero row for col´ >= col
    return (0, 0)
end

function nextnonzero_row_col(m::AbstractSparseMatrix, itr, col = first(itr.colrange))
    rows = rowvals(m)
    for col´ in col:last(itr.colrange)
        ptridx = findfirst(p -> isvalidrowcol(rows[p], col´, m, itr), nzrange(m, col´))
        ptridx === nothing || return (ptridx, col´)
    end
    # (0, 0) is sentinel for "no non-zero row for col´ >= col
    return (0, 0)
end

isvalidrowcol(row, col, m, itr) = row in itr.rowrange && !iszero(m[row, col])

function Base.iterate(itr::EachIndexNonzeros{<:Hamiltonian}, (ptridx, col, nhar) = firststate(itr, 1))
    nhar > length(itr.h.harmonics) && return nothing
    har = itr.h.harmonics[nhar]
    i = _iterate(har.h, itr, ptridx, col)
    if i === nothing
        nhar´ = nhar + 1
        return nhar´ > length(itr.h.harmonics) ? nothing : iterate(itr, firststate(itr, nhar´))
    else
        ((row, col), (ptridx, col)) = i
        return (row, col, har.dn), (ptridx, col, nhar)
    end
end

firststate(itr::EachIndexNonzeros{<:HamiltonianHarmonic}) = nextnonzero_row_col(itr.h.h, itr)

function Base.iterate(itr::EachIndexNonzeros{<:HamiltonianHarmonic}, (row, col) = firststate(itr))
    _iterate(itr.h.h, itr, row, col)
end

# Returns nothing or ((row, col), (nextrow, nextcol)), where row and nextrow can be a ptridx
function _iterate(m::AbstractSparseMatrix, itr, ptridx, col)
    col in itr.colrange || return nothing  # will also return nothing if col == 0 (sentinel)
    ptrs = nzrange(m, col)
    rows = rowvals(m)
    if ptridx <= length(ptrs)
        row = rows[ptrs[ptridx]]
        row in itr.rowrange && return (row, col), (ptridx + 1, col)
    end
    ptridx´, col´ = nextnonzero_row_col(m, itr, col + 1)
    return _iterate(m, itr, ptridx´, col´)
end

function _iterate(m::DenseMatrix, itr, row, col)
    col in itr.colrange || return nothing
    for row´ in row:last(itr.rowrange)
        iszero(m[row´, col]) || return (row´, col), (row´ + 1, col)
    end
    row´, col´ = nextnonzero_row_col(m, itr, col + 1)
    return _iterate(m, itr, row´, col´)
end

Base.IteratorSize(::EachIndexNonzeros) = Base.SizeUnknown()
Base.IteratorEltype(::EachIndexNonzeros) = Base.HasEltype()
Base.eltype(s::EachIndexNonzeros{<:Hamiltonian}) = Tuple{Int, Int, typeof(first(s.h.harmonics).dn)}
Base.eltype(s::EachIndexNonzeros{<:HamiltonianHarmonic}) = Tuple{Int, Int}

# stored_indices(h::Hamiltonian) = ((har.dn, rowvals(har.h)[ptr], col) for har in h.harmonics
#                                   for col in 1:size(har.h, 2) for ptr in nzrange(har.h, col))

# External API #
"""
    hamiltonian(lat[, model]; orbitals, type)

Create a `Hamiltonian` by additively applying `model::TighbindingModel` to the lattice `lat`
(see `hopping` and `onsite` for details on building tightbinding models).

The number of orbitals on each sublattice can be specified by the keyword `orbitals`
(otherwise all sublattices have one orbital by default). The following, and obvious
combinations, are possible formats for the `orbitals` keyword:

    orbitals = :a                # all sublattices have 1 orbital named :a
    orbitals = (:a,)             # same as above
    orbitals = (:a, :b, 3)       # all sublattices have 3 orbitals named :a and :b and :3
    orbitals = ((:a, :b), (:c,)) # first sublattice has 2 orbitals, second has one
    orbitals = ((:a, :b), :c)    # same as above
    orbitals = (Val(2), Val(1))  # same as above, with automatic names
    orbitals = (:A => (:a, :b), :D => :c) # sublattice :A has two orbitals, :D and rest have one
    orbitals = :D => Val(4)      # sublattice :D has four orbitals, rest have one

The matrix sizes of tightbinding `model` must match the orbitals specified. Internally, we
define a block size `N = max(num_orbitals)`. If `N = 1` (all sublattices with one orbital)
the the Hamiltonian element type is `type`. Otherwise it is `SMatrix{N,N,type}` blocks,
padded with the necessary zeros as required. Keyword `type` is `Complex{T}` by default,
where `T` is the number type of `lat`.

    h(ϕ₁, ϕ₂, ...)
    h((ϕ₁, ϕ₂, ...))

Build the Bloch Hamiltonian matrix `bloch(h, (ϕ₁, ϕ₂, ...))` of a `h::Hamiltonian` on an
`L`D lattice. (See also `bloch!` for a non-allocating version of `bloch`.)

    lat |> hamiltonian(model[, funcmodel]; kw...)

Functional `hamiltonian` form equivalent to `hamiltonian(lat, model[, funcmodel]; kw...)`.

# Indexing

Indexing into a Hamiltonian `h` works as follows. Access the `HamiltonianHarmonic` matrix at
a given `dn::NTuple{L,Int}` with `h[dn]`. Assign `v` into element `(i,j)` of said matrix
with `h[dn][i,j] = v` or `h[dn, i, j] = v`. Broadcasting with vectors of indices `is` and
`js` is supported, `h[dn][is, js] = v_matrix`.

To add an empty harmonic with a given `dn::NTuple{L,Int}`, do `push!(h, dn)`. To delete it,
do `deleteat!(h, dn)`.

# Examples
```jldoctest
julia> h = hamiltonian(LatticePresets.honeycomb(), hopping(@SMatrix[1 2; 3 4], range = 1/√3), orbitals = Val(2))
Hamiltonian{<:Lattice} : 2D Hamiltonian on a 2D Lattice in 2D space
  Bloch harmonics  : 5 (SparseMatrixCSC, sparse)
  Harmonic size    : 2 × 2
  Orbitals         : ((:a, :a), (:a, :a))
  Element type     : 2 × 2 blocks (Complex{Float64})
  Onsites          : 0
  Hoppings         : 6
  Coordination     : 3.0

julia> push!(h, (3,3)) # Adding a new Hamiltonian harmonic (if not already present)
Hamiltonian{<:Lattice} : 2D Hamiltonian on a 2D Lattice in 2D space
  Bloch harmonics  : 6 (SparseMatrixCSC, sparse)
  Harmonic size    : 2 × 2
  Orbitals         : ((:a, :a), (:a, :a))
  Element type     : 2 × 2 blocks (Complex{Float64})
  Onsites          : 0
  Hoppings         : 6
  Coordination     : 3.0

julia> h[3,3][1,1] = @SMatrix[1 2; 2 1]; h[3,3] # element assignment
2×2 SparseArrays.SparseMatrixCSC{StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4},Int64} with 1 stored entry:
  [1, 1]  =  [1.0+0.0im 2.0+0.0im; 2.0+0.0im 1.0+0.0im]

julia> h[3,3][[1,2],[1,2]] .= rand(SMatrix{2,2,Float64}, 2, 2) # Broadcast assignment
2×2 view(::SparseArrays.SparseMatrixCSC{StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4},Int64}, [1, 2], [1, 2]) with eltype StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4}:
 [0.271152+0.0im 0.921417+0.0im; 0.138212+0.0im 0.525911+0.0im]  [0.444284+0.0im 0.280035+0.0im; 0.565106+0.0im 0.121869+0.0im]
 [0.201126+0.0im 0.912446+0.0im; 0.372099+0.0im 0.931358+0.0im]  [0.883422+0.0im 0.874016+0.0im; 0.296095+0.0im 0.995861+0.0im]

julia> hopfunc(;k = 0) = hopping(k); hamiltonian(LatticePresets.square(), onsite(1) + hopping(2), hopfunc) # Parametric Hamiltonian
Parametric Hamiltonian{<:Lattice} : 2D Hamiltonian on a 2D Lattice in 2D space
  Bloch harmonics  : 5 (SparseMatrixCSC, sparse)
  Harmonic size    : 1 × 1
  Orbitals         : ((:a,),)
  Elements         : scalars (Complex{Float64})
  Onsites          : 1
  Hoppings         : 4
  Coordination     : 4.0

```
"""
hamiltonian(lat, ts...; orbitals = missing, kw...) =
    _hamiltonian(lat, sanitize_orbs(orbitals, lat.unitcell.names), ts...; kw...)
_hamiltonian(lat::AbstractLattice, orbs; kw...) = _hamiltonian(lat, orbs, TightbindingModel(); kw...)
_hamiltonian(lat::AbstractLattice, orbs, f::Function; kw...) = _hamiltonian(lat, orbs, TightbindingModel(), f; kw...)
_hamiltonian(lat::AbstractLattice, orbs, m::TightbindingModel; type::Type = Complex{numbertype(lat)}, kw...) =
    hamiltonian_sparse(blocktype(orbs, type), lat, orbs, m; kw...)

hamiltonian(t::TightbindingModel...; kw...) =
    z -> hamiltonian(z, t...; kw...)
hamiltonian(f::Function, t::TightbindingModel...; kw...) =
    z -> hamiltonian(z, f, t...; kw...)
hamiltonian(h::Hamiltonian) =
    z -> hamiltonian(z, h)

(h::Hamiltonian)(phases...) = bloch(h, phases...)

sanitize_orbs(o::Union{Val,NameType,Integer}, names::NTuple{N}) where {N} =
    ntuple(n -> sanitize_orbs(o), Val(N))
sanitize_orbs(o::NTuple{M,Union{NameType,Integer}}, names::NTuple{N}) where {M,N} =
    (ont = nametype.(o); ntuple(n -> ont, Val(N)))
sanitize_orbs(o::Missing, names) = sanitize_orbs((:a,), names)
sanitize_orbs(o::Pair, names) = sanitize_orbs((o,), names)
sanitize_orbs(os::NTuple{M,Pair}, names::NTuple{N}) where {N,M} =
    ntuple(Val(N)) do n
        for m in 1:M
            first(os[m]) == names[n] && return sanitize_orbs(os[m])
        end
        return (:a,)
    end
sanitize_orbs(os::NTuple{M,Any}, names::NTuple{N}) where {N,M} =
    ntuple(n -> n > M ? (:a,) : sanitize_orbs(os[n]), Val(N))

sanitize_orbs(p::Pair) = sanitize_orbs(last(p))
sanitize_orbs(o::Integer) = (nametype(o),)
sanitize_orbs(o::NameType) = (o,)
sanitize_orbs(o::Val{N}) where {N} = ntuple(_ -> :a, Val(N))
sanitize_orbs(o::NTuple{N,Union{Integer,NameType}}) where {N} = nametype.(o)
sanitize_orbs(p) = throw(ArgumentError("Wrong format for orbitals, see `hamiltonian`"))

Base.Matrix(h::Hamiltonian) = Hamiltonian(h.lattice, Matrix.(h.harmonics), h.orbitals)
Base.Matrix(h::HamiltonianHarmonic) = HamiltonianHarmonic(h.dn, Matrix(h.h))

Base.copy(h::Hamiltonian) = Hamiltonian(copy(h.lattice), copy.(h.harmonics), h.orbitals)
Base.copy(h::HamiltonianHarmonic) = HamiltonianHarmonic(h.dn, copy(h.h))

Base.size(h::Hamiltonian, n) = size(first(h.harmonics).h, n)
Base.size(h::Hamiltonian) = size(first(h.harmonics).h)
Base.size(h::HamiltonianHarmonic, n) = size(h.h, n)
Base.size(h::HamiltonianHarmonic) = size(h.h)

function LinearAlgebra.ishermitian(h::Hamiltonian)
    for hh in h.harmonics
        isnonnegative(hh.dn) || continue
        isassigned(h, -hh.dn) || return false
        hh.h == h[-hh.dn]' || return false
    end
    return true
end

bravais(h::Hamiltonian) = bravais(h.lattice)

issemibounded(h::Hamiltonian) = issemibounded(h.lattice)

nsites(h::Hamiltonian) = isempty(h.harmonics) ? 0 : nsites(first(h.harmonics))
nsites(h::HamiltonianHarmonic) = size(h.h, 1)

nsublats(h::Hamiltonian) = nsublats(h.lattice)

norbitals(h::Hamiltonian) = length.(h.orbitals)

# External API #

"""
    transform!(h::Hamiltonian, f::Function)

Transform the site positions of the Hamiltonian's lattice in place without modifying the
Hamiltonian harmonics.
"""
function transform!(h::Hamiltonian, f::Function)
    transform!(h.lattice, f)
    return h
end

# Indexing #

Base.push!(h::Hamiltonian{<:Any,L}, dn::NTuple{L,Int}) where {L} = push!(h, SVector(dn...))
Base.push!(h::Hamiltonian{<:Any,L}, dn::Vararg{Int,L}) where {L} = push!(h, SVector(dn...))
function Base.push!(h::Hamiltonian{<:Any,L,M,A}, dn::SVector{L,Int}) where {L,M,A}
    for hh in h.harmonics
        hh.dn == dn && return hh
    end
    hh = HamiltonianHarmonic{L,M,A}(dn, size(h)...)
    push!(h.harmonics, hh)
    return h
end

Base.getindex(h::Hamiltonian, dn::NTuple) = getindex(h, SVector(dn))
@inline function Base.getindex(h::Hamiltonian{<:Any,L}, dn::SVector{L,Int}) where {L}
    nh = findfirst(hh -> hh.dn == dn, h.harmonics)
    nh === nothing && throw(BoundsError(h, dn))
    return h.harmonics[nh].h
end
Base.getindex(h::Hamiltonian, dn::NTuple, i::Vararg{Int}) = h[dn][i...]
Base.getindex(h::Hamiltonian{LA, L}, i::Vararg{Int}) where {LA,L} = h[zero(SVector{L,Int})][i...]

Base.deleteat!(h::Hamiltonian{<:Any,L}, dn::Vararg{Int,L}) where {L} =
    deleteat!(h, toSVector(dn))
Base.deleteat!(h::Hamiltonian{<:Any,L}, dn::NTuple{L,Int}) where {L} =
    deleteat!(h, toSVector(dn))
function Base.deleteat!(h::Hamiltonian{<:Any,L}, dn::SVector{L,Int}) where {L}
    nh = findfirst(hh -> hh.dn == SVector(dn...), h.harmonics)
    nh === nothing || deleteat!(h.harmonics, nh)
    return h
end

Base.isassigned(h::Hamiltonian, dn::Vararg{Int}) = isassigned(h, SVector(dn))
Base.isassigned(h::Hamiltonian, dn::NTuple) = isassigned(h, SVector(dn))
Base.isassigned(h::Hamiltonian{<:Any,L}, dn::SVector{L,Int}) where {L} =
    findfirst(hh -> hh.dn == dn, h.harmonics) != nothing

#######################################################################
# auxiliary types
#######################################################################
struct IJV{L,M}
    dn::SVector{L,Int}
    i::Vector{Int}
    j::Vector{Int}
    v::Vector{M}
end

struct IJVBuilder{L,M,E,T,O,LA<:AbstractLattice{E,L,T}}
    lat::LA
    orbs::O
    ijvs::Vector{IJV{L,M}}
    kdtrees::Vector{KDTree{SVector{E,T},Euclidean,T}}
end

IJV{L,M}(dn::SVector{L} = zero(SVector{L,Int})) where {L,M} =
    IJV(dn, Int[], Int[], M[])

function IJVBuilder(lat::AbstractLattice{E,L,T}, orbs, ijvs::Vector{IJV{L,M}}) where {E,L,T,M}
    kdtrees = Vector{KDTree{SVector{E,T},Euclidean,T}}(undef, nsublats(lat))
    return IJVBuilder(lat, orbs, ijvs, kdtrees)
end

IJVBuilder(lat::AbstractLattice{E,L}, orbs, ::Type{M}) where {E,L,M} =
    IJVBuilder(lat, orbs, IJV{L,M}[])

function IJVBuilder(lat::AbstractLattice{E,L}, orbs, hs::Hamiltonian...) where {E,L}
    M = promote_blocktype(hs...)
    ijvs = IJV{L,M}[]
    builder = IJVBuilder(lat, orbs, ijvs)
    offset = 0
    for h in hs
        for har in h.harmonics
            ijv = builder[har.dn]
            push_block!(ijv, har, offset)
        end
        offset += size(h, 1)
    end
    return builder
end

function Base.getindex(b::IJVBuilder{L,M}, dn::SVector{L2,Int}) where {L,L2,M}
    L == L2 || throw(error("Tried to apply an $L2-dimensional model to an $L-dimensional lattice"))
    for e in b.ijvs
        e.dn == dn && return e
    end
    e = IJV{L,M}(dn)
    push!(b.ijvs, e)
    return e
end

Base.length(h::IJV) = length(h.i)
Base.isempty(h::IJV) = length(h) == 0
Base.copy(h::IJV) = IJV(h.dn, copy(h.i), copy(h.j), copy(h.v))

function Base.resize!(h::IJV, n)
    resize!(h.i, n)
    resize!(h.j, n)
    resize!(h.v, n)
    return h
end

Base.push!(ijv::IJV, (i, j, v)::Tuple) = (push!(ijv.i, i); push!(ijv.j, j); push!(ijv.v, v))

function push_block!(ijv::IJV{L,M}, h::HamiltonianHarmonic, offset) where {L,M}
    I, J, V = findnz(h.h)
    for (i, j, v) in zip(I, J, V)
        push!(ijv, (i + offset, j + offset, padtotype(v, M)))
    end
    return ijv
end

#######################################################################
# hamiltonian_sparse
#######################################################################
function hamiltonian_sparse(Mtype, lat, orbs, model)
    builder = IJVBuilder(lat, orbs, Mtype)
    return hamiltonian_sparse!(builder, lat, orbs, model)
end

function hamiltonian_sparse!(builder::IJVBuilder{L,M}, lat::AbstractLattice{E,L}, orbs, model) where {E,L,M}
    checkmodelorbs(model, orbs, lat)
    applyterms!(builder, terms(model)...)
    n = nsites(lat)
    HT = HamiltonianHarmonic{L,M,SparseMatrixCSC{M,Int}}
    harmonics = HT[HT(e.dn, sparse(e.i, e.j, e.v, n, n)) for e in builder.ijvs if !isempty(e)]
    return Hamiltonian(lat, harmonics, orbs, n, n)
end

applyterms!(builder, terms...) = foreach(term -> applyterm!(builder, term), terms)

applyterm!(builder::IJVBuilder, term::Union{OnsiteTerm, HoppingTerm}) =
    applyterm!(builder, term, sublats(term, builder.lat))

function applyterm!(builder::IJVBuilder{L,M}, term::OnsiteTerm, termsublats) where {L,M}
    selector = term.selector
    lat = builder.lat
    for s in termsublats
        is = siterange(lat, s)
        dn0 = zero(SVector{L,Int})
        ijv = builder[dn0]
        offset = lat.unitcell.offsets[s]
        for i in is
            isinregion(i, dn0, selector.region, lat) || continue
            r = lat.unitcell.sites[i]
            vs = orbsized(term(r,r), builder.orbs[s])
            v = padtotype(vs, M)
            term.forcehermitian ? push!(ijv, (i, i, 0.5 * (v + v'))) : push!(ijv, (i, i, v))
        end
    end
    return nothing
end

function applyterm!(builder::IJVBuilder{L,M}, term::HoppingTerm, termsublats) where {L,M}
    selector = term.selector
    checkinfinite(selector)
    lat = builder.lat
    for (s1, s2) in termsublats
        is, js = siterange(lat, s1), siterange(lat, s2)
        dns = dniter(selector.dns, Val(L))
        for dn in dns
            addadjoint = term.forcehermitian
            foundlink = false
            ijv = builder[dn]
            addadjoint && (ijvc = builder[negative(dn)])
            for j in js
                sitej = lat.unitcell.sites[j]
                rsource = sitej - lat.bravais.matrix * dn
                itargets = targets(builder, selector.range, rsource, s1)
                for i in itargets
                    isselfhopping((i, j), (s1, s2), dn) && continue
                    isinregion((i, j), (dn, zero(dn)), selector.region, lat) || continue
                    foundlink = true
                    rtarget = lat.unitcell.sites[i]
                    r, dr = _rdr(rsource, rtarget)
                    vs = orbsized(term(r, dr), builder.orbs[s1], builder.orbs[s2])
                    v = padtotype(vs, M)
                    if addadjoint
                        v *= redundancyfactor(dn, (s1, s2), selector)
                        push!(ijv, (i, j, v))
                        push!(ijvc, (j, i, v'))
                    else
                        push!(ijv, (i, j, v))
                    end
                end
            end
            foundlink && acceptcell!(dns, dn)
        end
    end
    return nothing
end

orbsized(m, orbs) = orbsized(m, orbs, orbs)
orbsized(m, o1::NTuple{D1}, o2::NTuple{D2}) where {D1,D2} = padtotype(m, SMatrix{D1,D2})

dniter(dns::Missing, ::Val{L}) where {L} = BoxIterator(zero(SVector{L,Int}))
dniter(dns, ::Val) = dns

function targets(builder, range::Real, rsource, s1)
    if !isassigned(builder.kdtrees, s1)
        sites = view(builder.lat.unitcell.sites, siterange(builder.lat, s1))
        (builder.kdtrees[s1] = KDTree(sites))
    end
    targets = inrange(builder.kdtrees[s1], rsource, range)
    targets .+= builder.lat.unitcell.offsets[s1]
    return targets
end

targets(builder, range::Missing, rsource, s1) = eachindex(builder.lat.sublats[s1].sites)

checkinfinite(selector) =
    selector.dns === missing && (selector.range === missing || !isfinite(selector.range)) &&
    throw(ErrorException("Tried to implement an infinite-range hopping on an unbounded lattice"))

isselfhopping((i, j), (s1, s2), dn) = i == j && s1 == s2 && iszero(dn)

# Avoid double-counting hoppings when adding adjoint
redundancyfactor(dn, ss, selector) =
    isnotredundant(dn, selector) || isnotredundant(ss, selector) ? 1.0 : 0.5
isnotredundant(dn::SVector, selector) = selector.dns !== missing && !iszero(dn)
isnotredundant((s1, s2)::Tuple{Int,Int}, selector) = selector.sublats !== missing && s1 != s2

#######################################################################
# unitcell/supercell for Hamiltonians
#######################################################################
function supercell(ham::Hamiltonian, args...; kw...)
    slat = supercell(ham.lattice, args...; kw...)
    return Hamiltonian(slat, ham.harmonics, ham.orbitals)
end

function unitcell(ham::Hamiltonian{<:Lattice}, args...; modifiers = (), kw...)
    sham = supercell(ham, args...; kw...)
    return unitcell(sham; modifiers = modifiers)
end

function unitcell(ham::Hamiltonian{LA,L}; modifiers = ()) where {E,L,T,L´,LA<:Superlattice{E,L,T,L´}}
    lat = ham.lattice
    sc = lat.supercell
    modifiers´ = resolve.(ensuretuple(modifiers), Ref(lat))
    mapping = OffsetArray{Int}(undef, sc.sites, sc.cells.indices...) # store supersite indices newi
    mapping .= 0
    foreach_supersite((s, oldi, olddn, newi) -> mapping[oldi, Tuple(olddn)...] = newi, lat)
    dim = nsites(sc)
    B = blocktype(ham)
    S = typeof(SparseMatrixBuilder{B}(dim, dim))
    harmonic_builders = HamiltonianHarmonic{L´,B,S}[]
    pinvint = pinvmultiple(sc.matrix)
    foreach_supersite(lat) do s, source_i, source_dn, newcol
        for oldh in ham.harmonics
            rows = rowvals(oldh.h)
            vals = nonzeros(oldh.h)
            target_dn = source_dn + oldh.dn
            super_dn = new_dn(target_dn, pinvint)
            wrapped_dn = wrap_dn(target_dn, super_dn, sc.matrix)
            newh = get_or_push!(harmonic_builders, super_dn, dim, newcol)
            for p in nzrange(oldh.h, source_i)
                target_i = rows[p]
                # check: wrapped_dn could exit bounding box along non-periodic direction
                checkbounds(Bool, mapping, target_i, Tuple(wrapped_dn)...) || continue
                newrow = mapping[target_i, Tuple(wrapped_dn)...]
                val = applymodifiers(vals[p], lat, (source_i, target_i), (source_dn, target_dn), modifiers´...)
                iszero(newrow) || pushtocolumn!(newh.h, newrow, val)
            end
        end
        foreach(h -> finalizecolumn!(h.h), harmonic_builders)
    end
    harmonics = [HamiltonianHarmonic(h.dn, sparse(h.h)) for h in harmonic_builders]
    unitlat = unitcell(lat)
    orbs = ham.orbitals
    return Hamiltonian(unitlat, harmonics, orbs)
end

function get_or_push!(hs::Vector{<:HamiltonianHarmonic{L,B,<:SparseMatrixBuilder}}, dn, dim, currentcol) where {L,B}
    for h in hs
        h.dn == dn && return h
    end
    newh = HamiltonianHarmonic(dn, SparseMatrixBuilder{B}(dim, dim))
    currentcol > 1 && finalizecolumn!(newh.h, currentcol - 1) # for columns that have been already processed
    push!(hs, newh)
    return newh
end

wrap_dn(olddn::SVector, newdn::SVector, supercell::SMatrix) = olddn - supercell * newdn

applymodifiers(val, lat, inds, dns) = val

function applymodifiers(val, lat, inds, dns, m::ElementModifier{Val{false}}, ms...)
    selected = m.selector(lat, inds, dns)
    val´ = selected ? m.f(val) : val
    return applymodifiers(val´, lat, inds, dns, ms...)
end

function applymodifiers(val, lat, (row, col), (dnrow, dncol), m::Onsite!{Val{true}}, ms...)
    selected = m.selector(lat, (row, col), (dnrow, dncol))
    if selected
        r = sites(lat)[col] + bravais(lat) * dncol
        val´ = selected ? m.f(val, r) : val
    else
        val´ = val
    end
    return applymodifiers(val´, lat, (row, col), (dnrow, dncol), ms...)
end

function applymodifiers(val, lat, (row, col), (dnrow, dncol), m::Hopping!{Val{true}}, ms...)
    selected = m.selector(lat, (row, col), (dnrow, dncol))
    if selected
        br = bravais(lat)
        r, dr = _rdr(sites(lat)[col] + br * dncol, sites(lat)[row] + br * dnrow)
        val´ = selected ? m.f(val, r, dr) : val
    else
        val´ = val
    end
    return applymodifiers(val´, lat, (row, col), (dnrow, dncol), ms...)
end

#######################################################################
# wrap
#######################################################################
"""
    wrap(h::Hamiltonian, axis::Int; factor = 1)

Build a new Hamiltonian wherein the Bravais `axis` is wrapped into a loop. If a `factor` is
given, the wrapped hoppings will be multiplied by said factor. This is useful to represent a
flux Φ through the loop, if `factor = exp(im * 2π * Φ/Φ₀)`.

    h |> wrap(axis; kw...)

Functional form equivalent to `wrap(h, axis; kw...)`.

# Examples
```
julia> LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = 1/√3)) |> unitcell((1,-1), (10, 10)) |> wrap(2)
Hamiltonian{<:Lattice} : 1D Hamiltonian on a 1D Lattice in 2D space
  Bloch harmonics  : 3 (SparseMatrixCSC, sparse)
  Harmonic size    : 40 × 40
  Orbitals         : ((:a,), (:a,))
  Element type     : scalar (Complex{Float64})
  Onsites          : 0
  Hoppings         : 120
  Coordination     : 3.0
```
"""
function wrap(h::Hamiltonian{<:Lattice,L}, axis; factor = 1) where {L}
    1 <= axis <= L || throw(ArgumentError("wrap axis should be between 1 and the lattice dimension $L"))
    lattice´ = _wrap(h.lattice, axis)
    harmonics´ = _wrap(h.harmonics, axis, factor)
    return Hamiltonian(lattice´, harmonics´, h.orbitals)
end

wrap(axis; kw...) = h -> wrap(h, axis; kw...)

_wrap(lat::Lattice, axis) = Lattice(_wrap(lat.bravais, axis), lat.unitcell)

function _wrap(br::Bravais{E,L}, axis) where {E,L}
    mask = deleteat(SVector{L}(1:L), axis)
    return Bravais(br.matrix[:, mask], br.semibounded[mask])
end

function _wrap(harmonics::Vector{HamiltonianHarmonic{L,M,A}}, axis, factor) where {L,M,A}
    harmonics´ = HamiltonianHarmonic{L-1,M,A}[]
    for har in harmonics
        dn = har.dn
        dn´ = deleteat(dn, axis)
        factor´ = iszero(dn[axis]) ? 1 : factor
        add_or_push!(harmonics´, dn´, har.h, factor´)
    end
    return harmonics´
end

function add_or_push!(hs::Vector{<:HamiltonianHarmonic}, dn, matrix::AbstractMatrix, factor)
    for h in hs
        if h.dn == dn
            h.h .+= factor .* matrix
            return h
        end
    end
    newh = HamiltonianHarmonic(dn, matrix)
    push!(hs, newh)
    return newh
end

#######################################################################
# combine
#######################################################################
"""
    combine(hams::Hamiltonian...; coupling = missing)

Build a new Hamiltonian `h` that combines all `hams` as diagonal blocks, and applies
`coupling::Model`, if provided, to build the off-diagonal couplings. Note that the diagonal
blocks are not modified by the coupling model.
"""
combine(hams::Hamiltonian...; coupling = missing) = _combine(coupling, hams...)

_combine(::Missing, hams...) = _combine(TightbindingModel(), hams...)

function _combine(model::TightbindingModel, hams::Hamiltonian...)
    lat = combine((h -> h.lattice).(hams)...)
    orbs = tuplejoin((h -> h.orbitals).(hams)...)
    builder = IJVBuilder(lat, orbs, hams...)
    model´ = offdiagonal(model, lat, nsublats.(hams))
    ham = hamiltonian_sparse!(builder, lat, orbs, model´)
    return ham
end

#######################################################################
# Bloch routines
#######################################################################
struct SupercellBloch{L,T,H<:Hamiltonian{<:Superlattice}}
    hamiltonian::H
    phases::SVector{L,T}
    axis::Int
end
SupercellBloch(h, ϕs) = SupercellBloch(h, ϕs, 0)

Base.summary(h::SupercellBloch{L,T}) where {L,T} =
    "SupercellBloch{$L)}: Bloch Hamiltonian matrix lazily defined on an $(L)D supercell"

function Base.show(io::IO, sb::SupercellBloch)
    ioindent = IOContext(io, :indent => string("  "))
    print(io, summary(sb), "
  Phases          : $(Tuple(sb.phases))
  Axis            : $(iszero(sb.axis) ? "none" : sb.axis)\n")
    print(ioindent, sb.hamiltonian.lattice.supercell)
end

"""
    bloch!(matrix, h::Hamiltonian, ...)

In-place version of `bloch`. Overwrite `matrix` with the Bloch Hamiltonian matrix of `h` for
the specified Bloch phases `ϕs` (see `bloch` for definition and API).  A conventient way to
obtain a `matrix` is to use `similarmatrix(h)`, which will return an `AbstractMatrix` of the
same type as the Hamiltonian's. Note, however, that matrix need not be of the same type
(e.g. it can be dense, while `h` is sparse).

# Examples
```
julia> h = LatticePresets.honeycomb() |> hamiltonian(onsite(1), hopping(2));

julia> bloch!(similarmatrix(h), h, (.2,.3), 2)  # velocity operator along `bravais(h)[2]`
2×2 SparseArrays.SparseMatrixCSC{Complex{Float64},Int64} with 4 stored entries:
  [1, 1]  =  1.99001-0.199667im
  [2, 1]  =  1.96013-0.397339im
  [1, 2]  =  1.96013+0.397339im
  [2, 2]  =  1.99001-0.199667im
```

# See also:
    bloch, optimize!, similarmatrix
"""
bloch!(matrix, h, ϕs::Number...) = _bloch!(matrix, h, toSVector(ϕs), 0)
bloch!(matrix, h, ϕs::Tuple, axis = 0) = _bloch!(matrix, h, toSVector(ϕs), axis)
bloch!(matrix, h, ϕs::SVector, axis = 0) = _bloch!(matrix, h, ϕs, axis)
bloch!(matrix, ϕs::Number...) = h -> bloch!(matrix, h, ϕs...)
bloch!(matrix, ϕs::Tuple, axis = 0) = h -> bloch!(matrix, h, ϕs, axis)

function _bloch!(matrix::AbstractMatrix, h::Hamiltonian{<:Lattice,L,M}, ϕs, axis::Number) where {L,M}
    rawmatrix = parent(matrix)
    if iszero(axis)
        _copy!(rawmatrix, first(h.harmonics).h) # faster copy!(dense, sparse) specialization
        add_harmonics!(rawmatrix, h, ϕs, dn -> 1)
    else
        fill!(rawmatrix, zero(M)) # There is no guarantee of same structure
        add_harmonics!(rawmatrix, h, ϕs, dn -> -im * dn[axis])
    end
    return matrix
end

function _bloch!(matrix::AbstractMatrix, h::Hamiltonian{<:Lattice,L,M}, ϕs, dnfunc::Function) where {L,M}
    prefactor0 = dnfunc(zero(ϕs))
    rawmatrix = parent(matrix)
    if iszero(prefactor0)
        fill!(rawmatrix, zero(M))
    else
        _copy!(rawmatrix, first(h.harmonics).h)
        rmul!(rawmatrix, prefactor0)
    end
    add_harmonics!(rawmatrix, h, ϕs, dnfunc)
    return matrix
end

add_harmonics!(zerobloch, h::Hamiltonian{<:Lattice}, ϕs::SVector{0}, _) = zerobloch

function add_harmonics!(zerobloch, h::Hamiltonian{<:Lattice,L}, ϕs::SVector{L}, dnfunc) where {L}
    ϕs´ = ϕs'
    for ns in 2:length(h.harmonics)
        hh = h.harmonics[ns]
        hhmatrix = hh.h
        prefactor = dnfunc(hh.dn)
        iszero(prefactor) && continue
        ephi = prefactor * cis(-ϕs´ * hh.dn)
        _add!(zerobloch, hhmatrix, ephi)
    end
    return zerobloch
end

"""
    bloch(h::Hamiltonian{<:Lattice}, ϕs::Real...)

Build the Bloch Hamiltonian matrix of `h`, for the specified Bloch phases `ϕs`. In terms of
Bloch wavevector `k`, `ϕs = k * bravais(h)`, it is defined as Overwrite `matrix` with the
Bloch Hamiltonian matrix of `h`, for the specified Bloch phases `ϕs`, defined as
`H(ϕs) = ∑exp(-im * ϕs' * dn) h_dn` where `h_dn` are Bloch harmonics connecting unit cells
at a distance `dR = bravais(h) * dn`.

    bloch(h::Hamiltonian{<:Lattice})

Build the intra-cell Hamiltonian matrix of `h`, without adding any Bloch harmonics.

    bloch(h::Hamiltonian{<:Lattice}, ϕs::NTuple{L,Real}[, axis::Int = 0])

A nonzero `axis` produces the derivative of the Bloch matrix respect to `ϕs[axis]` (i.e. the
velocity operator along this axis), `∂H(ϕs) = ∑ -im * dn[axis] * exp(-im * ϕs' * dn) h_dn`

    bloch(matrix, h::Hamiltonian{<:Lattice}, ϕs::NTuple{L,Real}, dnfunc::Function)

Generalization that applies a prefactor `dnfunc(dn) * exp(im * ϕs' * dn)` to the `dn`
harmonic.

    h |> bloch(ϕs...)
    h(ϕs...)

Functional forms of `bloch`, equivalent to `bloch(h, ϕs...)`

    bloch(h::Hamiltonian{<:Superlattice}, ϕs::Number...)
    bloch(h::Hamiltonian{<:Superlattice}, ϕs::Tuple[, axis = 0])

Build a `SupercellBloch` object that lazily implements the Bloch Hamiltonian in the
`Superlattice` without actually building the matrix (e.g. for matrix-free diagonalization).

# Notes

`bloch` allocates a new matrix on each call. For a non-allocating version of `bloch`, see
`bloch!`.

# Examples
```
julia> h = LatticePresets.honeycomb() |> hamiltonian(onsite(1), hopping(2)) |> bloch(.2,.3)
2×2 SparseArrays.SparseMatrixCSC{Complex{Float64},Int64} with 4 stored entries:
  [1, 1]  =  1.99001-0.199667im
  [2, 1]  =  1.96013-0.397339im
  [1, 2]  =  1.96013+0.397339im
  [2, 2]  =  1.99001-0.199667im
```

# See also:
    bloch!, optimize!, similarmatrix
"""
bloch(ϕs::Number...) = h -> bloch(h, ϕs...)
bloch(ϕs::Tuple, axis = 0) = h -> bloch(h, ϕs, axis)
bloch(h::Hamiltonian{<:Lattice}, args...) = bloch!(similarmatrix(h), h, args...)
bloch(h::Hamiltonian{<:Superlattice}, ϕs::Number...) =
    SupercellBloch(h, toSVector(ϕs))
bloch(h::Hamiltonian{<:Superlattice}, ϕs::Tuple, axis = 0) =
    SupercellBloch(h, toSVector(ϕs), axis)

"""
    similarmatrix(h::Hamiltonian; optimize = true)

Create an uninitialized matrix of the same type of the Hamiltonian's matrix, calling
`optimize!(h)` first if `optimize = true` to produce an optimal work matrix in the sparse
case.
"""
function similarmatrix(h::Hamiltonian; optimize = true)
    optimize && optimize!(h)
    sm = size(h)
    T = eltype(h)
    matrix = similar(h.harmonics[1].h, T, sm[1], sm[2])
    return matrix
end

"""
    optimize!(h::Hamiltonian)

Prepare a sparse Hamiltonian `h` to increase the performance of subsequent calls to
`bloch(h, ϕs...)` and `bloch!(matrix, h, ϕs...)` by minimizing memory reshufflings. It also
adds missing structural zeros to the diagonal to enable shifts by `α*I` (for
shift-and-invert methods).

No optimization will be performed on non-sparse Hamiltonians, or those defined on
`Superlattice`s, for which Bloch Hamiltonians are lazily evaluated.

Note that when calling `similarmatrix(h)` on a sparse `h`, `optimize!` is called first.

# See also:
    bloch, bloch!
"""
function optimize!(ham::Hamiltonian{<:Lattice,L,M,A}) where {L,M,A<:SparseMatrixCSC}
    h0 = first(ham.harmonics)
    n, m = size(h0.h)
    iszero(h0.dn) || throw(ArgumentError("First Hamiltonian harmonic is not the fundamental"))
    nh = length(ham.harmonics)
    builder = SparseMatrixBuilder{M}(n, m)
    for col in 1:m
        for i in eachindex(ham.harmonics)
            h = ham.harmonics[i].h
            for p in nzrange(h, col)
                v = i == 1 ? nonzeros(h)[p] : zero(M)
                row = rowvals(h)[p]
                pushtocolumn!(builder, row, v, false) # skips repeated rows
            end
        end
        pushtocolumn!(builder, col, zero(M), false) # if not present already, add structural zeros to diagonal
        finalizecolumn!(builder)
    end
    ho = sparse(builder)
    copy!(h0.h, ho) # Inject new structural zeros into zero harmonics
    return ham
end
# IDEA: could sum and subtract all harmonics instead
# Tested, it is slower

function optimize!(ham::Hamiltonian{<:Lattice,L,M,A}) where {L,M,A<:AbstractMatrix}
    # @warn "Hamiltonian is not sparse. Nothing changed."
    return ham
end

function optimize!(ham::Hamiltonian{<:Superlattice})
    # @warn "Hamiltonian is defined on a Superlattice. Nothing changed."
    return ham
end

"""
    flatten(h::Hamiltonian)

Flatten a multiorbital Hamiltonian `h` into one with a single orbital per site. The
associated lattice is flattened also, so that there is one site per orbital for each initial
site (all at the same position). Note that in the case of sparse Hamiltonians, zeros in
hopping/onsite matrices are preserved as structural zeros upon flattening.

    h |> flatten()

Functional form equivalent to `flatten(h)` of `h |> flatten` (included for consistency with
the rest of the API).

# Examples
```
julia> h = LatticePresets.honeycomb() |> hamiltonian(hopping(@SMatrix[1 2], range = 1/√3, sublats = (:A,:B)), orbitals = (Val(1), Val(2)))
Hamiltonian{<:Lattice} : 2D Hamiltonian on a 2D Lattice in 2D space
  Bloch harmonics  : 5 (SparseMatrixCSC, sparse)
  Harmonic size    : 2 × 2
  Orbitals         : ((:a,), (:a, :a))
  Element type     : 2 × 2 blocks (Complex{Float64})
  Onsites          : 0
  Hoppings         : 6
  Coordination     : 3.0

julia> flatten(h)
Hamiltonian{<:Lattice} : 2D Hamiltonian on a 2D Lattice in 2D space
  Bloch harmonics  : 5 (SparseMatrixCSC, sparse)
  Harmonic size    : 3 × 3
  Orbitals         : ((:flat,), (:flat,))
  Element type     : scalar (Complex{Float64})
  Onsites          : 0
  Hoppings         : 12
  Coordination     : 4.0
```
"""
flatten() = h -> flatten(h)

function flatten(h::Hamiltonian)
    all(isequal(1), norbitals(h)) && return copy(h)
    harmonics´ = [flatten(har, h.orbitals, h.lattice) for har in h.harmonics]
    lattice´ = flatten(h.lattice, h.orbitals)
    orbitals´ = (_ -> (:flat, )).(h.orbitals)
    return Hamiltonian(lattice´, harmonics´, orbitals´)
end

flatten(h::HamiltonianHarmonic, orbs, lat) =
    HamiltonianHarmonic(h.dn, _flatten(h.h, length.(orbs), lat))

function _flatten(src::SparseMatrixCSC{<:SMatrix{N,N,T}}, norbs::NTuple{S,<:Any}, lat) where {N,T,S}
    offsets´ = flattenoffsets(lat.unitcell.offsets, norbs)
    dim´ = last(offsets´)

    builder = SparseMatrixBuilder{T}(dim´, dim´, nnz(src) * N * N)

    for col in 1:size(src, 2)
        scol = sublat(lat, col)
        for j in 1:norbs[scol]
            for p in nzrange(src, col)
                row = rowvals(src)[p]
                srow = sublat(lat, row)
                rowoffset´ = flattenind(row, lat, norbs, offsets´) - 1
                val = nonzeros(src)[p]
                for i in 1:norbs[srow]
                    pushtocolumn!(builder, rowoffset´ + i, val[i, j])
                end
            end
            finalizecolumn!(builder, false)
        end
    end
    matrix = sparse(builder)
    return matrix
end

function _flatten(src::DenseMatrix{<:SMatrix{N,N,T}}, norbs::NTuple{S,<:Any}, lat) where {N,T,S}
    offsets´ = flattenoffsets(lat.unitcell.offsets, norbs)
    dim´ = last(offsets´)
    matrix = similar(src, T, dim´, dim´)

    for col in 1:size(src, 2), row in 1:size(src, 1)
        srow, scol = sublat(lat, row), sublat(lat, col)
        nrow, ncol = norbs[srow], norbs[scol]
        val = src[row, col]
        rowoffset´ = flattenind(row, lat, norbs, offsets´) - 1
        coloffset´ = flattenind(col, lat, norbs, offsets´) - 1
        for j in 1:ncol, i in 1:nrow
            matrix[rowoffset´ + i, coloffset´ + j] = val[i, j]
        end
    end
    return matrix
end

# sublat offsets after flattening (without padding zeros)
function flattenoffsets(offsets, norbs)
    offsets´ = similar(offsets)
    i = 0
    for s in eachindex(norbs)
        offsets´[s] = i
        ns = offsets[s + 1] - offsets[s]
        no = norbs[s]
        i += ns * no
    end
    offsets´[end] = i
    return offsets´
end

# index of first orbital of site i after flattening
function flattenind(i, lat, norbs, offsets´)
    s = sublat(lat, i)
    offset = lat.unitcell.offsets[s]
    Δi = i - offset
    i´ = offsets´[s] + (Δi - 1) * norbs[s] + 1
    return i´
end

function flatten(lat::Lattice, orbs)
    length(orbs) == nsublats(lat) || throw(ArgumentError("Msmatch between sublattices and orbitals"))
    unitcell´ = flatten(lat.unitcell, length.(orbs))
    bravais´ = lat.bravais
    lat´ = Lattice(bravais´, unitcell´)
end

function flatten(unitcell::Unitcell, norbs::NTuple{S,Int}) where {S}
    offsets´ = flattenoffsets(unitcell.offsets, norbs)
    ns´ = last(offsets´)
    sites´ = similar(unitcell.sites, ns´)
    i = 1
    for sl in 1:S, site in sites(unitcell, sl), rep in 1:norbs[sl]
        sites´[i] = site
        i += 1
    end
    names´ = unitcell.names
    unitcell´ = Unitcell(sites´, names´, offsets´)
    return unitcell´
end

max_issues_repo_name: pablosanjose/QBox

max_forks_repo_name: pablosanjose/QBox


text:
using PEGParser
using Compat

immutable RangeType
  first
  last
  step
end

@grammar lang begin
  start = list(stmt, NEWLINE)
  stmt = if_stmt | for_stmt | decl_stmt |
    func_stmt | func_call_stmt | return_stmt | block
  if_stmt = -IF + condition + stmt + *(else_stmt)
  else_stmt = ELSE + stmt
  block = (LBRACE + block_list + RBRACE){ _2 }
  block_list = list(stmt, NEWLINE) #{ children }
  atom = string | range | number | name
  value = expr | atom

  condition = value + comparison + value
  comparison = EQ | NE | LTE | GTE | LT | GT

  string = (-SPACE + r"\"[^\"]*\""){ _1.value }

  number = float_number | integer_number
  integer_number = (-SPACE + integer){ parse(Int, _1.value) }
  float_number = (-SPACE + float){ parse(Float64, _1.value) }

  bool = TRUE | FALSE
  range = (number + COLON + number){ RangeType(_1, _3, 1) }
  name = (-SPACE + r"[a-zA-Z_][0-9a-zA-Z_]*"){ symbol(_1.value) }

  expr_op = term + op1 + expr
  expr = expr_op | term
  term_op = factor + op2 + term

  term = term_op | factor
  factor = atom | pfactor
  pfactor = (LPAREN + expr + RPAREN){ _2 }
  op1 = (ADD | SUB){ symbol(_1.value) }
  op2 = (MULT | DIV){ symbol(_1.value) }

  for_stmt = -FOR + name + -IN + value + stmt

  decl_stmt = var_stmt | const_stmt

  var_type = name
  var_name_type = name + -":" + var_type
  var_stmt = -VAR + (var_name_type | name) + -EQUALS + value
  const_stmt = -CONST + (var_name_type | name) + -EQUALS + value

  func_stmt = -FUNCTION + name + -LPAREN + param_list + -RPAREN + ?(-RARROW + var_type) + block
  param_list = list(var_name_type, COMMA, 0)

  func_call_stmt = name + -LPAREN + param_call_list + -RPAREN
  param_call_list = list(value, COMMA, 0)

  return_stmt = -RETURN + value

  ADD = (-SPACE + "+"){ _1 }
  SUB = (-SPACE + "-"){ _1 }
  MULT = (-SPACE + "*"){ _1 }
  DIV = (-SPACE + "/"){ _1 }
  FUNCTION = SPACE + "function"
  IF = SPACE + "if"
  FOR = SPACE + "for"
  IN = SPACE + "in"
  VAR = SPACE + "var"
  CONST = SPACE + "const"
  TRUE = SPACE + "true"
  FALSE = SPACE + "false"
  RETURN = SPACE + "return"
  EQUALS = SPACE + "="
  EQ = SPACE + "=="
  NEQ = SPACE + "!="
  LT = SPACE + "<"
  GT = SPACE + ">"
  LTE = SPACE + "<="
  GTE = SPACE + ">="
  ELSE = SPACE + "else"
  LBRACE = SPACE + "{"
  RBRACE = SPACE + "}"
  LPAREN = SPACE + "("
  RPAREN = SPACE + ")"
  RARROW = SPACE + "->"
  SPACE = r"[ \n\r\t]*"
  NEWLINE = r"[\n\r]+" | SEMICOLON
  COLON = SPACE + ":"
  SEMICOLON = SPACE + ";"
  COMMA = SPACE + ","
end

toexpr(node, children, ::MatchRule{:default}) = node

toexpr(node, children, ::MatchRule{:start}) = Expr(:block, children...)

toexpr(node, children, ::MatchRule{:term_op}) =
  Expr(:call, children[2], children[1], children[3])

toexpr(node, children, ::MatchRule{:expr_op}) =
  Expr(:call, children[2], children[1], children[3])

toexpr(node, children, ::MatchRule{:var_stmt}) =
  Expr(:(=), children[1], children[2])

toexpr(node, children, ::MatchRule{:const_stmt}) =
  Expr(:const, Expr(:(=), children[1], children[2]))

toexpr(node, children, ::MatchRule{:func_call_stmt}) =
  Expr(:call, children[1], children[2]...)

toexpr(node, children, ::MatchRule{:block_list}) =
  Expr(:block, children...)

toexpr(node, children, ::MatchRule{:param_call_list}) = children

toexpr(node, children, ::MatchRule{:for_stmt}) =
  Expr(:for,
    Expr(:(=), children[1],
      Expr(:(:), children[2].first, children[2].last)),
    children[3])

data = """
const k = 20
var j = (k+5) / 2.0
for i in 0:10 {
  println(i/2.0)
}
"""

(ast, pos, error) = parse(lang, data)
println(ast)
println(error)

code = transform(toexpr, ast)
println(eval(code))

ext: jl


text:
# Autogenerated wrapper script for libpcap_jll for armv7l-linux-musleabihf
export libpcap

JLLWrappers.@generate_wrapper_header("libpcap")
JLLWrappers.@declare_library_product(libpcap, "libpcap.so.1")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libpcap,
        "lib/libpcap.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Revise

include("Interface.jl")

module Attractors

using ..Interface: @ODE

export lorenz, TSUCS1, TSUCS2, yuwang, aizawa,
       rucklidge, genesio_tesi, finance, shimizu_morioka,
       noose_hoover, liu_chen, arneodo, bouali,
       burke_shaw, chen_celikovsky, chen_lee, dequan_li,
       hadley, halvorsen, newton_leipnik

lorenz = @ODE quote
    dx = σ * (y - x)
    dy = x * (ρ - z) - y
    dz = x * y - β * z

    σ = 10
    ρ = 28
    β = 8/3
end

TSUCS1 = @ODE quote
    dx = α * (y - x) + δ * x * z
    dy = ζ * y - x * z
    dz = β * z + x * y - ϵ * x * x

    α = 40
    β = 0.833
    δ = 0.5
    ϵ = 0.65
    ζ = 20
end

TSUCS2 = @ODE quote
    dx = α * (y - x) + δ * x * z
    dy = (ς * x) - (x * z) + ζ * y
    dz = (β * z) + (x * y) - (ϵ * x * x)

    α = 40
    δ = 0.16
    ς = 55
    ζ = 20
    β = 1.833
    ϵ = 0.65
end

yuwang = @ODE quote
    dx = α * (y - x)
    dy = β * x - σ * x * z
    dz = exp(x * y) - δ * z

    α = 10
    β = 40
    σ = 2
    δ = 2.5
end

aizawa = @ODE quote
    dx = (z - β) * x - δ * y
    dy = δ * x + (z - β) * y
    dz = σ + α * z - (z^3)/3 - (x^2 + y^2)*(1 + ϵ * z) + ζ * z * x^3

    α = 0.95
    β = 0.7
    σ = 0.6
    δ = 3.5
    ϵ = 0.25
    ζ = 0.1
end

lorenz_mod_2 = @ODE quote
    dx = -α*x + y^2 - z^2 + α*σ
    dy = x*(y-β*z) + δ
    dz = -z + x*(β * y + z)

    α = 0.9
    β = 5
    σ = 9.9
    δ = 1
end

thomas = @ODE quote
    dx = -β * x + sin(y)
    dy = -β * y + sin(z)
    dz = -β * z + sin(x)

    β = 0.19
end

rucklidge = @ODE quote
    dx = -κ*x + α*y - y* z
    dy = x
    dz = -z + y^2

    κ = 2
    α = 6.7
end

genesio_tesi = @ODE quote
    dx = y
    dy = z
    dz = -σ * x - β * y - α * z + x^2

    α = 0.44
    β = 1.1
    σ = 1.0
end

finance = @ODE quote
    dx = ((1/β) - α)*x + z + x*y
    dy = -β*y - x^2
    dz = -x - σ*z

    α = 0.001
    β = 0.2
    σ = 1.1
end

shimizu_morioka = @ODE quote
    dx = y
    dy = (1-z)*x - α*y
    dz = x^2 - β*z

    α = 0.75
    β = 0.45
end

noose_hoover = @ODE quote
    dx = y
    dy = -x + y*z
    dz = α - y^2

    α = 1.5
end

liu_chen = @ODE quote
    dx = α*y + β*x + σ*y*z
    dy = δ*y - z + ϵ*x*z
    dz = ζ*z + η*x*y

    α = 2.4
    β = -3.78
    σ = 14
    δ = -11
    ϵ = 4
    ζ = 5.58
    η = 1
end

arneodo = @ODE quote
    dx = y
    dy = z
    dz = -α*x - β*y - z + σ*x^3

    α = -5.5
    β = 3.5
    σ = -1.0
end

bouali = @ODE quote
    dx = x*(4-y) + α*z
    dy = -y*(1-x^2)
    dz = -x*(1.5 - σ*z) - 0.05*z

    α = 0.3
    σ = 1.0
end

burke_shaw = @ODE quote
    dx = -σ*(x+y)
    dy = -y -σ*x*z
    dz = σ*x*y + ν

    σ = 10
    ν = 4.272
end

chen_celikovsky = @ODE quote
    dx = α*(y-x)
    dy = -x*z + σ*y
    dz = x*y - β*z

    α = 36
    β = 3
    σ = 20
end

chen_lee = @ODE quote
    dx = α*x - y*z
    dy = β*y + x*z
    dz = σ*z + x*(y/3)

    α = 5
    β = -10
    σ = -0.38
end

dequan_li = TSUCS2

hadley = @ODE quote
    dx = -y^2 -z^2 -α*x + α*ζ
    dy = x*y - β*x*z - y + η
    dz = β*x*y + x*z - z

    α = 0.2
    β = 4
    ζ = 8
    η = 1
end

halvorsen = @ODE quote
    dx = -α*x -4*y -4*z - y^2
    dy = -α*y - 4*z - 4*x - z^2
    dz = -α*z - 4*x - 4*y - x^2

    α = 1.4
end

newton_leipnik = @ODE quote
    dx = -α*x + y + 10*y*z
    dy = -x - 0.4*y + 5*x*z
    dz = β*z -5*x*y

    α = 0.4
    β = 0.175
end

end


text:
"Partitioned Differential Algebraic Equation"
immutable PDAE{T} <: Equation{T}
    m::Int
    n::Int
    f::Function
    g::Function
    u::Function
    v::Function
    ϕ::Function
    t₀::T
    q₀::Array{T,1}
    p₀::Array{T,1}
    λ₀::Array{T,1}

    function PDAE(m, n, f, g, u, v, ϕ, t₀, q₀, p₀, λ₀)
        @assert m == length(q₀)
        @assert m == length(p₀)
        @assert n == length(λ₀)
        @assert 2m ≥ n
        @assert T == eltype(q₀)
        @assert T == eltype(p₀)
        @assert T == eltype(λ₀)

        new(m, n, f, g, u, v, ϕ, t₀, q₀, p₀, λ₀)
    end
end

function PDAE{T}(f::Function, g::Function, u::Function, v::Function, ϕ::Function,
                 t₀::Real, q₀::Vector{T}, p₀::Vector{T}, λ₀::Vector{T})
    @assert length(q₀) == length(p₀)
    PDAE{T}(length(q₀), length(λ₀), f, g, u, v, ϕ, t₀, q₀, p₀, λ₀)
end

function PDAE{T}(f::Function, g::Function, u::Function, v::Function, ϕ::Function,
                 q₀::Vector{T}, p₀::Vector{T}, λ₀::Vector{T})
    PDAE(f, g, u, v, ϕ, 0, q₀, p₀, λ₀)
end


text:
#Wrappers for functions directly in the rospy namespace
export init_node, is_shutdown, spin,
       get_param, has_param, set_param, delete_param,
       logdebug, loginfo, logwarn, logerr, logfatal

"""
    init_node(name; args...)

Initialize this node, registering it with the ROS master. All arguments are passed on directly to
the rospy init_node function.
"""
init_node(name::AbstractString; args...) =
    __rospy__.init_node(ascii(name); args...)

"""
    is_shutdown()

Return the shutdown status of the node.
"""
is_shutdown()          = __rospy__.is_shutdown()

get_published_topics() = __rospy__.get_published_topics()
get_ros_root()         = __rospy__.get_ros_root()

"""
    spin()

Block execution and process callbacks/service calls until the node is shut down.
"""
function spin()
    #Have to make sure both Julia tasks and python threads can wake up so
    #can't just call rospy's spin
    while ! is_shutdown()
        rossleep(Duration(0.001))
    end
end

#Parameter server API
"""
    get_param(param_name, default=nothing)

Request the value of a parameter from the parameter server, with optional default value. If no
default is given, throws a `KeyError` if the parameter cannot be found.
"""
function get_param(param_name::AbstractString, def=nothing)
    try
        if def == nothing
            __rospy__.get_param(ascii(param_name))
        else
            __rospy__.get_param(ascii(param_name), def)
        end
    catch ex
        throw(KeyError(pycall(pybuiltin("str"), PyAny, ex.val)[2:end-1]))
    end
end

"""
    set_param(param_name, val)

Set the value of a parameter on the parameter server.
"""
set_param(param_name::AbstractString, val) =
    __rospy__.set_param(ascii(param_name), val)

"""
    has_param(param_name)

Return a boolean specifying if a parameter exists on the parameter server.
"""
has_param(param_name::AbstractString) =
    __rospy__.has_param(ascii(param_name))

"""
    delete_param(param_name)

Delete a parameter from the parameter server. Throws a `KeyError` if no such parameter exists.
"""
function delete_param(param_name::AbstractString)
    try
        __rospy__.delete_param(ascii(param_name))
    catch ex
        throw(KeyError(pycall(pybuiltin("str"), PyAny, ex.val)[2:end-1]))
    end
end

#Doesn't work for some reason
#rospy_search_param(param_name::AbstractString) =
#    __rospy__.rospy_search_param(ascii(param_name))
get_param_names() = __rospy__.get_param_names()

#Logging API
logdebug(msg, args...) = __rospy__.logdebug(msg, args...)
loginfo(msg, args...)  = __rospy__.loginfo(msg, args...)
logwarn(msg, args...)  = __rospy__.logwarn(msg, args...)
logerr(msg, args...)   = __rospy__.logerr(msg, args...)
logfatal(msg, args...) = __rospy__.logfatal(msg, args...)

"""
    logdebug, loginfo, logwarn, logerr, logfatal

Call the rospy logging system at the corresponding message level, passing a message and other
arguments directly.
"""
logdebug, loginfo, logwarn, logerr, logfatal

#Node information
get_name()             = __rospy__.get_name()
get_namespace()        = __rospy__.get_namespace()
get_node_uri()         = __rospy__.get_node_uri()
get_caller_id()        = __rospy__.get_caller_id()


text:
# obs_weight is a shortcut function for getting the relative likelihood of an observation without having to construct the observation distribution. Useful for particle filtering
# maintained by @zsunberg

"""
    obs_weight(pomdp, sp, o)
    obs_weight(pomdp, a, sp, o)
    obs_weight(pomdp, s, a, sp, o)

Return a weight proportional to the likelihood of receiving observation o from state sp (and a and s if they are present).

This is a useful shortcut for particle filtering so that the observation distribution does not have to be represented.
"""
function obs_weight end

@generated function obs_weight(p, s, a, sp, o)
    ow_impl = :(obs_weight(p, a, sp, o))
    o_impl = :(pdf(observation(p, s, a, sp), o))
    if implemented(obs_weight, Tuple{p, a, sp, o})
        return ow_impl
    elseif implemented(observation, Tuple{p, s, a, sp})
        return o_impl
    else
        return quote
            try # trick to get the compiler to put the right backedges in
                $ow_impl
                $o_impl
            catch
                throw(MethodError(obs_weight, (p,s,a,sp,o)))
            end
        end
    end
end

@generated function obs_weight(p, a, sp, o)
    ow_impl = :(obs_weight(p, sp, o))
    o_impl = :(pdf(observation(p, a, sp), o))
    if implemented(obs_weight, Tuple{p, sp, o})
        return ow_impl
    elseif implemented(observation, Tuple{p, a, sp})
        return o_impl
    else
        return quote
            try # trick to get the compiler to put the right backedges in
                $ow_impl
                $o_impl
            catch
                throw(MethodError(obs_weight, (p, a, sp, o)))
            end
        end
    end
end

@generated function obs_weight(p, sp, o)
    impl = :(pdf(observation(p, sp), o))
    if implemented(observation, Tuple{p, sp})
        return impl
    else
         return quote
            try # trick to get the compiler to put the right backedges in
                $impl
            catch
                return :(throw(MethodError(obs_weight, (p, sp, o))))
            end
        end
    end
end

function implemented(f::typeof(obs_weight), TT::Type)
    m = which(f, TT)
    if length(TT.parameters) == 5
        P, S, A, _, O = TT.parameters
        reqs_met = implemented(observation, Tuple{P,S,A,S}) || implemented(obs_weight, Tuple{P,A,S,O})
    elseif length(TT.parameters) == 4
        P, A, S, O = TT.parameters
        reqs_met = implemented(observation, Tuple{P,A,S}) || implemented(obs_weight, Tuple{P,S,O})
    elseif length(TT.parameters) == 3
        P, S, O = TT.parameters
        reqs_met = implemented(observation, Tuple{P,S})
    else
        return method_exists(f, TT)
    end
    if m.module == POMDPModelTools && !reqs_met
        return false
    else
        true
    end
end


text:
mutable struct Printer
    inputpath::Channel{String}
    errorpath::Channel{String}
    inkremaining::Int32
    reserve::Printer
    name::String
    function Printer(ch1, ch2, ink, name)
        this = new()
        this.inputpath = ch1
        this.errorpath = ch2
        this.inkremaining = ink
        this.name = name
        this.reserve = this
        this
    end
end

function lineprintertask(printer)
    while true
        line = take!(printer.inputpath)
        linesprinted = 0
        if(printer.inkremaining < 1)
            if(printer.reserve == printer)
                put!(printer.errorpath, "Error: printer $(printer.name) out of ink")
            else
                put!(printer.reserve.inputpath, line)
            end
        else
            println(line)
            printer.inkremaining -= 1
        end
    end
end

function schedulework(poems)
    printerclose(printer) = (close(printer.inputpath); close(printer.errorpath))
    reserveprinter = Printer(Channel{String}(1), Channel{String}(10), 5, "Reserve")
    mainprinter = Printer(Channel{String}(1), Channel{String}(10), 5, "Main")
    mainprinter.reserve = reserveprinter
    @async(lineprintertask(mainprinter))
    @async(lineprintertask(reserveprinter))
    printers = [mainprinter, reserveprinter]
    activeprinter = 1
    @sync(
    for poem in poems
        activeprinter = (activeprinter % length(printers)) + 1
        @async(
        for line in poem
            put!(printers[activeprinter].inputpath, line)
        end)
    end)
    for p in printers
        while isready(p.errorpath)
            println(take!(p.errorpath))
        end
        printerclose(p)
    end
end

const humptydumpty = ["Humpty Dumpty sat on a wall.",
                      "Humpty Dumpty had a great fall.",
                      "All the king's horses and all the king's men,",
                      "Couldn't put Humpty together again."]

const oldmothergoose = ["Old Mother Goose,",
                        "When she wanted to wander,",
                        "Would ride through the air,",
                        "On a very fine gander.",
                        "Jack's mother came in,",
                        "And caught the goose soon,",
                        "And mounting its back,",
                        "Flew up to the moon."]

schedulework([humptydumpty, oldmothergoose])

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAdmissionregistrationV1beta1RuleWithOperations <: SwaggerModel
    apiGroups::Any # spec type: Union{ Nothing, Vector{String} } # spec name: apiGroups
    apiVersions::Any # spec type: Union{ Nothing, Vector{String} } # spec name: apiVersions
    operations::Any # spec type: Union{ Nothing, Vector{String} } # spec name: operations
    resources::Any # spec type: Union{ Nothing, Vector{String} } # spec name: resources
    scope::Any # spec type: Union{ Nothing, String } # spec name: scope

    function IoK8sApiAdmissionregistrationV1beta1RuleWithOperations(;apiGroups=nothing, apiVersions=nothing, operations=nothing, resources=nothing, scope=nothing)
        o = new()
        validate_property(IoK8sApiAdmissionregistrationV1beta1RuleWithOperations, Symbol("apiGroups"), apiGroups)
        setfield!(o, Symbol("apiGroups"), apiGroups)
        validate_property(IoK8sApiAdmissionregistrationV1beta1RuleWithOperations, Symbol("apiVersions"), apiVersions)
        setfield!(o, Symbol("apiVersions"), apiVersions)
        validate_property(IoK8sApiAdmissionregistrationV1beta1RuleWithOperations, Symbol("operations"), operations)
        setfield!(o, Symbol("operations"), operations)
        validate_property(IoK8sApiAdmissionregistrationV1beta1RuleWithOperations, Symbol("resources"), resources)
        setfield!(o, Symbol("resources"), resources)
        validate_property(IoK8sApiAdmissionregistrationV1beta1RuleWithOperations, Symbol("scope"), scope)
        setfield!(o, Symbol("scope"), scope)
        o
    end
end # type IoK8sApiAdmissionregistrationV1beta1RuleWithOperations

const _property_map_IoK8sApiAdmissionregistrationV1beta1RuleWithOperations = Dict{Symbol,Symbol}(Symbol("apiGroups")=>Symbol("apiGroups"), Symbol("apiVersions")=>Symbol("apiVersions"), Symbol("operations")=>Symbol("operations"), Symbol("resources")=>Symbol("resources"), Symbol("scope")=>Symbol("scope"))
const _property_types_IoK8sApiAdmissionregistrationV1beta1RuleWithOperations = Dict{Symbol,String}(Symbol("apiGroups")=>"Vector{String}", Symbol("apiVersions")=>"Vector{String}", Symbol("operations")=>"Vector{String}", Symbol("resources")=>"Vector{String}", Symbol("scope")=>"String")
Base.propertynames(::Type{ IoK8sApiAdmissionregistrationV1beta1RuleWithOperations }) = collect(keys(_property_map_IoK8sApiAdmissionregistrationV1beta1RuleWithOperations))
Swagger.property_type(::Type{ IoK8sApiAdmissionregistrationV1beta1RuleWithOperations }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAdmissionregistrationV1beta1RuleWithOperations[name]))}
Swagger.field_name(::Type{ IoK8sApiAdmissionregistrationV1beta1RuleWithOperations }, property_name::Symbol) =  _property_map_IoK8sApiAdmissionregistrationV1beta1RuleWithOperations[property_name]

function check_required(o::IoK8sApiAdmissionregistrationV1beta1RuleWithOperations)
    true
end

function validate_property(::Type{ IoK8sApiAdmissionregistrationV1beta1RuleWithOperations }, name::Symbol, val)
end


text:
# Module for The King is Dead bot
# This houses all of the functions

function random_pick_from_supply(board::Board)
    supply_array = Array{String}(undef, 0)
    for (f_name, f_num) in board.supply
        append!(supply_array, repeat([f_name], f_num))
    end
    pick = supply_array[rand(1:length(supply_array))]
    return pick
end

function initialise_player_courts!(players::Array{Player}, board::Board)
    for player in players
        while player.court_size < court_size
            pick = random_pick_from_supply(board)
            player.court[pick] += 1
            player.court_size += 1
            board.supply[pick] -= 1
        end
    end
end

function initialise_board!(board::Board)
    # First make the players' courts
    for player in board.players
        while player.court_size < court_size
            pick = random_pick_from_supply(board)
            player.court[pick] += 1
            player.court_size += 1
            board.supply[pick] -= 1
        end
    end
    # Then put followers in home bases
    for region in board.regions
        for (f, starting_region) in faction_names
            if region.name == starting_region
                region.followers[f] += 2
                region.size = 2
                board.supply[f] -= 2 # Reduce the supply
                for r in board.unresolved_regions
                    if r.name == starting_region
                        r.controller = f
                    end
                end
            end
        end
    end
    # Then randomly pick followers from the supply
    for region in board.regions
        while region.size < region_size
            pick = random_pick_from_supply(board)
            region.followers[pick] += 1
            region.size += 1
            board.supply[pick] -= 1
        end
    end
end

function resolve_region!(board::Board, region::String)
    index = 0
    for (i, r) in enumerate(board.unresolved_regions)
        if r.name == region
            # Check which faction controls this region
            max_followers = maximum([v for (k,v) in r.followers])
            controlling_factions = [k for (k,v) in r.followers if v == max_followers]
            if length(controlling_factions) == 1
                r.controller = controlling_factions[1]
            else
                r.controller = "France"
                board.invasions += 1
            end
            r.resolved = true
            # Return all followers to the supply
            for (f_name, f_num) in r.followers
                board.supply[f_name] += f_num
            end
            index = i
        end
    end
    if index == 0
        error("The region to be resolved does not exist in the unresolved regions\n",
              "Tried to resolve ", region)
    end
    push!(board.resolved_regions, board.unresolved_regions[index])
    deleteat!(board.unresolved_regions, index)
end

function resolve_board!(board::Board)
    while length(board.unresolved_regions) > 0
        r = board.unresolved_regions[1]
        resolve_region!(board, r.name)
    end
end

function evaluate_current_board(board::Board)
    evaluation = 0
    player = board.players[board.turn]
    for r in board.resolved_regions
        evaluation += get(player.court, r.controller,
                        minimum([v for (k,v) in player.court]))
    end
    evaluation
end

function evaluate_resolved_board(board::Board)
    board_copy = deepcopy(board)
    resolve_board!(board_copy)
    evaluate_current_board(board_copy)
end

function take_follower!(board::Board, region::String,
                        faction::String)
    for r in board.unresolved_regions
        if r.name == region
            if r.followers[faction] == 0
                error("Trying to take a follower from a region where that faction is not present\n",
                "Tried to take a(n) ", faction, " follower from ", region)
            end
            r.followers[faction] -= 1
            r.size -= 1
            board.players[board.turn].court[faction] += 1
            board.players[board.turn].court_size += 1
        end
    end
end

function take_follower!(board::Board, player::Player,
    p::Pair{String, String})
    take_follower!(board, player, p.first, p.second)
end

function place_follower!(board::Board, region::String, faction::String)
    for r in board.unresolved_regions
        if r.name == region
            if board.supply[faction] == 0
                prinln("Trying to take a follower from the supply when that faction is not present\n",
                "Tried to take a(n) ", faction, " follower from the supply")
                break
            end
            r.followers[faction] += 1
            r.size += 1
            board.supply[faction] -= 1
        end
    end
end

function place_follower!(board::Board, d::Pair{String, String})
    place_follower!(board, d.second, d.first)
end

max_stars_repo_name: danielpicone/the-king-is-dead

max_issues_repo_name: danielpicone/the-king-is-dead

max_forks_repo_name: danielpicone/the-king-is-dead


text:
# high-level interface for CURAND
#
# the interface is split in two levels:
# - functions that extend the Random standard library, and take an RNG as first argument,
#   will only ever dispatch to CURAND and as a result are limited in the types they support.
# - functions that take an array will dispatch to either CURAND or GPUArrays
# - `cu`-prefixed functions are provided for constructing GPU arrays from only an eltype

## seeding

seed!(rng::RNG=generator()) = generate_seeds(rng)

## in-place

# uniform
Random.rand!(rng::RNG, A::CuArray{Float32}) = generate_uniform(rng, A)
Random.rand!(rng::RNG, A::CuArray{Float64}) = generate_uniform_double(rng, A)

# normal
Random.randn!(rng::RNG, A::CuArray{Float32}; mean=0, stddev=1) = generate_normal(rng, A, mean, stddev)
Random.randn!(rng::RNG, A::CuArray{Float64}; mean=0, stddev=1) = generate_normal_double(rng, A, mean, stddev)

# log-normal
rand_logn!(rng::RNG, A::CuArray{Float32}; mean=0, stddev=1) = generate_log_normal(rng, A, mean, stddev)
rand_logn!(rng::RNG, A::CuArray{Float64}; mean=0, stddev=1) = generate_log_normal_double(rng, A, mean, stddev)

# log-normal
rand_poisson!(rng::RNG, A::CuArray{Cuint}; lambda=1) = generate_poisson(rng, A, lambda)

## out of place

Random.rand(rng::RNG, ::Type{X}, dims::Dims) where {X} = rand!(rng, CuArray{X}(undef, dims))
Random.randn(rng::RNG, ::Type{X}, dims::Dims; kwargs...) where {X} = randn!(rng, CuArray{X}(undef, dims); kwargs...)
rand_logn(rng::RNG, ::Type{X}, dims::Dims; kwargs...) where {X} = rand_logn!(rng, CuArray{X}(undef, dims); kwargs...)
rand_poisson(rng::RNG, ::Type{X}, dims::Dims; kwargs...) where {X} = rand_poisson!(rng, CuArray{X}(undef, dims); kwargs...)

# specify default types
Random.rand(rng::RNG, dims::Integer...; kwargs...) = rand(rng, Float32, dims...; kwargs...)
Random.randn(rng::RNG, dims::Integer...; kwargs...) = randn(rng, Float32, dims...; kwargs...)
rand_logn(rng::RNG, dims::Integer...; kwargs...) = rand_logn(rng, Float32, dims...; kwargs...)
rand_poisson(rng::RNG, dims::Integer...; kwargs...) = rand_poisson(rng, Cuint, dims...; kwargs...)

# convenience
Random.randn(rng::RNG, ::Type{X}, dim1::Integer, dims::Integer...; kwargs...) where {X} =
    randn(rng, X, Dims((dim1, dims...)); kwargs...)
rand_logn(rng::RNG, ::Type{X}, dim1::Integer, dims::Integer...; kwargs...) where {X} =
    rand_logn(rng, X, Dims((dim1, dims...)); kwargs...)
rand_poisson(rng::RNG, ::Type{X}, dim1::Integer, dims::Integer...; kwargs...) where {X} =
    rand_poisson(rng, X, Dims((dim1, dims...)); kwargs...)

## functions that dispatch to either CURAND or GPUArrays

uniform_rng(::CuArray{<:Union{Float32,Float64}}) = generator()
uniform_rng(A::CuArray) = GPUArrays.global_rng(A)

normal_rng(::CuArray{<:Union{Float32,Float64}}) = generator()
normal_rng(::CuArray{T}) where {T} =
    error("CuArrays does not support generating normally distributed numbers of type $T")

logn_rng(::CuArray{<:Union{Float32,Float64}}) = generator()
logn_rng(::CuArray{T}) where {T} =
    error("CuArrays does not support generating lognormally distributed numbers of type $T")

poisson_rng(::CuArray{Cuint}) = generator()
poisson_rng(::CuArray{T}) where {T} =
    error("CuArrays does not support generating Poisson distributed numbers of type $T")

Random.rand!(A::CuArray; kwargs...) = rand!(uniform_rng(A), A; kwargs...)
Random.randn!(A::CuArray; kwargs...) = randn!(normal_rng(A), A; kwargs...)
rand_logn!(A::CuArray; kwargs...) = rand_logn!(logn_rng(A), A; kwargs...)
rand_poisson!(A::CuArray; kwargs...) = rand_poisson!(poisson_rng(A), A; kwargs...)
rand_logn(A::CuArray; kwargs...) = rand_logn!(logn_rng(A), A; kwargs...)
rand_poisson(A::CuArray; kwargs...) = rand_poisson!(poisson_rng(A), A; kwargs...)

rand(::Type{X}, args...; kwargs...) where {X} = rand!(CuArray{X}(undef, args...); kwargs...)
randn(::Type{X}, args...; kwargs...) where {X} = randn!(CuArray{X}(undef, args...); kwargs...)
rand_logn(::Type{X}, args...; kwargs...) where {X} = rand_logn!(CuArray{X}(undef, args...); kwargs...)
rand_poisson(::Type{X}, args...; kwargs...) where {X} = rand_poisson!(CuArray{X}(undef, args...); kwargs...)

# specify default types
rand(args...; kwargs...) where {X} = rand(Float32, args...; kwargs...)
randn(args...; kwargs...) where {X} = randn(Float32, args...; kwargs...)
rand_logn(args...; kwargs...) where {X} = rand_logn(Float32, args...; kwargs...)
rand_poisson(args...; kwargs...) where {X} = rand_poisson(Cuint, args...; kwargs...)


text:
"""
    tri_quadrature(deg; vertices, transform)

Calculate quadrature points in a triangle

@arg deg: polynomial degree
@arg vertices: vertex coordinates

- for a equilateral triangle, the vertices are ([-1.0, -1 / √3], [1.0, -1 / √3], [0.0, 2 / √3])
- for a right triangle, the vertices take ([-1.0, -1.0], [1.0, -1.0], [-1.0, 1.0])

@arg if transform equals to true, then we transform the coordinates from x-y to r-s plane
"""
function tri_quadrature(
    deg;
    vertices = ([-1.0, -1 / √3], [1.0, -1 / √3], [0.0, 2 / √3]),
    transform = true,
)

    py"""
    import quadpy.t2._williams_shunn_jameson as qwsj
    def create_quadrature(n):
        if n == 1:
            scheme = qwsj.williams_shunn_jameson_1()
        elif n == 2:
            scheme = qwsj.williams_shunn_jameson_2()
        elif n == 3:
            scheme = qwsj.williams_shunn_jameson_3()
        elif n == 4:
            scheme = qwsj.williams_shunn_jameson_4()
        elif n == 5:
            scheme = qwsj.williams_shunn_jameson_5()
        elif n == 6:
            scheme = qwsj.williams_shunn_jameson_6()
        elif n == 7:
            scheme = qwsj.williams_shunn_jameson_7()
        elif n == 8:
            scheme = qwsj.williams_shunn_jameson_8()
        else:
            print("Not a valid polynomial degree")

        return scheme.points, scheme.weights
    """

    # trilinear coordinates (三线坐标)
    points0, weights = py"create_quadrature"(deg + 1)

    # cartesian coordinates of vertices
    p1, p2, p3 = vertices

    a = norm(p2 .- p3)
    b = norm(p3 .- p1)
    c = norm(p2 .- p1)

    points = similar(points0, size(points0, 2), 2)
    for i in axes(points, 1)
        x, y, z = points0[:, i]

        # trilinear -> cartesian
        points[i, :] .= (a * x .* p1 + b * y .* p2 + c * z .* p3) / (a * x + b * y + c * z)
    end

    if transform == true
        r, s = xy_rs(points)
        points .= hcat(r, s)
    end

    return points, weights
end

"""
    triface_quadrature(N)

Calculate quadrature points along a face of right triangle

- face 1: 1 -> 2
- face 2: 2 -> 3
- face 3: 3 -> 1

Face 2 is assumed to be hypotenuse

@arg N: polynomial degree
"""
function triface_quadrature(N)
    Δf = [1.0, √2, 1.0]

    pf = Array{Float64}(undef, 3, N + 1, 2)
    wf = Array{Float64}(undef, 3, N + 1)

    p0, w0 = gausslegendre(N + 1)

    pf[1, :, 1] .= p0
    pf[2, :, 1] .= p0[end:-1:1]
    pf[3, :, 1] .= -1.0
    pf[1, :, 2] .= -1.0
    pf[2, :, 2] .= p0
    pf[3, :, 2] .= p0[end:-1:1]

    wf[1, :] .= w0 .* Δf[1]
    wf[2, :] .= w0 .* Δf[2]
    wf[3, :] .= w0 .* Δf[3]

    return pf, wf
end


text:
# MGH problem 7 - Helical valley function
#
#  Source:
#  J. J. Moré, B. S. Garbow and K. E. Hillstrom
#  Testing Unconstrained Optimization Software
#  ACM Transactions on Mathematical Software, 7(1):17-41, 1981
#
# A. Montoison, Montreal, 05/2018.

export mgh07

"Helical valley function"
function mgh07()

  nls  = Model()
  x0   = [-1.0,  0.0,  0.0]
  @variable(nls, x[i=1:3], start = x0[i])

  θ_aux(t) = (t > 0 ? 0.0 : 0.5)
  JuMP.register(nls, :θ_aux, 1, θ_aux, autodiff=true)

  @NLexpression(nls, F1, 10*(x[3] - 10*(atan(x[2]/x[1])/(2*π) + θ_aux(x[1]))))
  @NLexpression(nls, F2, 10*(sqrt(x[1]^2 + x[2]^2) - 1.0))
  @NLexpression(nls, F3, 1*x[3])

  return MathProgNLSModel(nls, [F1, F2, F3], name="mgh07")
end


text:
# Lets check normalization
using GeometryTypes
using Makie
using Colors
using JuliennedArrays

export spherical_plot, complex_spherical_plot

function sampleSphere(n::Int)
    map(normalize, Slices(randn(3,n),1))
end

function complexImage(vals::Array{Complex{Float64}})
    max = maximum(abs.(vals))
    map(x -> HSV(angle(x)*360/(2π),(abs(x)/max)^.5,1),vals)
end

function spherical_plot(fun::Function, n::Int)
    # n sets the resolution
    θ = [0;(0.5:n-0.5)/n;1]
    φ = [(0:2n-2)*2/(2n-1);2]
    x = [cospi(φ)*sinpi(θ) for θ in θ, φ in φ]
    y = [sinpi(φ)*sinpi(θ) for θ in θ, φ in φ]
    z = [cospi(θ) for θ in θ, φ in φ]
    vals = map(((x,y,z) -> fun( [x,y,z])), x,y,z);
    s = Makie.surface(x, y, z, color = complexImage(vals))
end

function comptosphere(z::Complex)
    # solve for the
    x = real(z)
    y = imag(z)
    t = 2/(1+abs(z)^2)
    [t*x,t*y,1-t]
end

function complex_spherical_plot(fun::Function, n::Int)
    # plots a function of the complex plane on the riemann sphere
    function spheretocomp(r::Array)
        # solve for the complex plane point
        x = r[1]
        y = r[2]
        t = 1-r[3]
        z = (x + im*y)/t
    end
    spherical_plot((x->fun(spheretocomp(x))), n)
end


text:
using InkLocalizer
using Test
using JSON

data = joinpath(@__DIR__, "data/TheIntercept.json")

a, b = InkLocalizer.localize(data)

output = joinpath(@__DIR__, "data/TheIntercept_localized.json")
localkeys = joinpath(@__DIR__, "data/TheIntercept_localkeys.json")
write(output, JSON.json(a))
write(localkeys, JSON.json(b, 2))


text:
# Julia wrapper for header: /usr/include/scip/pub_prop.h
# Automatically generated using Clang.jl wrap_c

function SCIPpropComp(elem1, elem2)
    ccall((:SCIPpropComp, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPpropCompPresol(elem1, elem2)
    ccall((:SCIPpropCompPresol, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPpropCompName(elem1, elem2)
    ccall((:SCIPpropCompName, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPpropGetData(prop)
    ccall((:SCIPpropGetData, libscip), Ptr{SCIP_PROPDATA}, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropSetData(prop, propdata)
    ccall((:SCIPpropSetData, libscip), Cvoid, (Ptr{SCIP_PROP}, Ptr{SCIP_PROPDATA}), prop, propdata)
end

function SCIPpropGetName(prop)
    ccall((:SCIPpropGetName, libscip), Cstring, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetDesc(prop)
    ccall((:SCIPpropGetDesc, libscip), Cstring, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetPriority(prop)
    ccall((:SCIPpropGetPriority, libscip), Cint, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetPresolPriority(prop)
    ccall((:SCIPpropGetPresolPriority, libscip), Cint, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetFreq(prop)
    ccall((:SCIPpropGetFreq, libscip), Cint, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetSetupTime(prop)
    ccall((:SCIPpropGetSetupTime, libscip), Cdouble, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropSetFreq(prop, freq)
    ccall((:SCIPpropSetFreq, libscip), Cvoid, (Ptr{SCIP_PROP}, Cint), prop, freq)
end

function SCIPpropGetTime(prop)
    ccall((:SCIPpropGetTime, libscip), Cdouble, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetStrongBranchPropTime(prop)
    ccall((:SCIPpropGetStrongBranchPropTime, libscip), Cdouble, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetRespropTime(prop)
    ccall((:SCIPpropGetRespropTime, libscip), Cdouble, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetPresolTime(prop)
    ccall((:SCIPpropGetPresolTime, libscip), Cdouble, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetNCalls(prop)
    ccall((:SCIPpropGetNCalls, libscip), Clonglong, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetNRespropCalls(prop)
    ccall((:SCIPpropGetNRespropCalls, libscip), Clonglong, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetNCutoffs(prop)
    ccall((:SCIPpropGetNCutoffs, libscip), Clonglong, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetNDomredsFound(prop)
    ccall((:SCIPpropGetNDomredsFound, libscip), Clonglong, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropIsDelayed(prop)
    ccall((:SCIPpropIsDelayed, libscip), UInt32, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropWasDelayed(prop)
    ccall((:SCIPpropWasDelayed, libscip), UInt32, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropIsInitialized(prop)
    ccall((:SCIPpropIsInitialized, libscip), UInt32, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetNFixedVars(prop)
    ccall((:SCIPpropGetNFixedVars, libscip), Cint, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetNAggrVars(prop)
    ccall((:SCIPpropGetNAggrVars, libscip), Cint, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetNChgVarTypes(prop)
    ccall((:SCIPpropGetNChgVarTypes, libscip), Cint, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetNChgBds(prop)
    ccall((:SCIPpropGetNChgBds, libscip), Cint, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetNAddHoles(prop)
    ccall((:SCIPpropGetNAddHoles, libscip), Cint, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetNDelConss(prop)
    ccall((:SCIPpropGetNDelConss, libscip), Cint, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetNAddConss(prop)
    ccall((:SCIPpropGetNAddConss, libscip), Cint, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetNUpgdConss(prop)
    ccall((:SCIPpropGetNUpgdConss, libscip), Cint, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetNChgCoefs(prop)
    ccall((:SCIPpropGetNChgCoefs, libscip), Cint, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetNChgSides(prop)
    ccall((:SCIPpropGetNChgSides, libscip), Cint, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetNPresolCalls(prop)
    ccall((:SCIPpropGetNPresolCalls, libscip), Cint, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetTimingmask(prop)
    ccall((:SCIPpropGetTimingmask, libscip), SCIP_PROPTIMING, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropDoesPresolve(prop)
    ccall((:SCIPpropDoesPresolve, libscip), UInt32, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropGetPresolTiming(prop)
    ccall((:SCIPpropGetPresolTiming, libscip), SCIP_PRESOLTIMING, (Ptr{SCIP_PROP},), prop)
end

function SCIPpropSetPresolTiming(prop, presoltiming)
    ccall((:SCIPpropSetPresolTiming, libscip), Cvoid, (Ptr{SCIP_PROP}, SCIP_PRESOLTIMING), prop, presoltiming)
end


text:
@testset "Global Quantities" begin

    # Test Residual
    T = Float64
    N = 10
    dt = 0.1
    p = 3
    i = 2
    model = UnicycleGame(p=p)
    probsize = ProblemSize(N,model)
    x0 = rand(SVector{model.n,T})
    opts = Options()

    Q = [Diagonal(rand(SVector{model.ni[i],T})) for i=1:p]
    R = [Diagonal(rand(SVector{model.mi[i],T})) for i=1:p]
    xf = [i*ones(SVector{model.ni[i],T}) for i=1:p]
    uf = [2i*ones(SVector{model.mi[i],T}) for i=1:p]
    game_obj = GameObjective(Q,R,xf,uf,N,model)

    game_con = GameConstraintValues(probsize)
    prob = GameProblem(N, dt, x0, model, opts, game_obj, game_con)
    @test typeof(prob) <: GameProblem
    residual!(prob)
    ibr_residual!(prob, i)

    # Test Residual Jacobian
    residual_jacobian!(prob)
    ibr_residual_jacobian!(prob, i)

    # Test scn
    @test scn(1234.0) == " 1.2e+3"
    @test scn(-1234.0) == "-1.2e+3"
    @test scn(-0.1234) == "-1.2e-1"
    @test scn( 0.1234) == " 1.2e-1"
    @test scn(0) == " 0.0e+0"
    @test scn(-0) == " 0.0e+0"
    @test scn(0, digits=3) == " 0.000e+0"
    @test scn(1234, digits=3) == " 1.234e+3"
    @test scn(1234, digits=0) == " 1e+3"
    @test_throws AssertionError scn(1234, digits=-1) == " 1e+3"

end


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8,6))
cm = countmap(purchase)
bar([1, 2], [cm["No"], cm["Yes"]])
xticks([1, 2], ["No", "Yes"], fontsize=12)
yticks(fontsize=12)
ylabel("Number of occurences", fontsize=14)

savefig(joinpath(@OUTPUT, "ISL-lab-4-bal2.svg")) # hide


text:
"""
    run([env, cicfg::CIService]; kwargs...)

Run the `RegistryCI.AutoMerge` service.

# Arguments

- `env`: an `AbstractDictionary` used to read environmental variables from.
   Defaults to `ENV` but a plain `Dict` can be passed to mimic an alternate environment.
- `ciccfg`: Configuration struct describing the continuous integration (CI) environment in which AutoMerge is being run.

# Keyword Arguments

- `merge_new_packages`: should AutoMerge merge registration PRs for new packages
- `merge_new_versions`: should AutoMerge merge registration PRs for new versions of packages
- `new_package_waiting_period`: new package waiting period, e.g `Day(3)`.
- `new_jll_package_waiting_period`: new JLL package waiting period, e.g `Minute(20)`.
- `new_version_waiting_period`: new package version waiting period, e.g `Minute(10)`.
- `new_jll_version_waiting_period`: new JLL package version waiting period, e.g `Minute(10)`.
- `registry`: the registry name you want to run AutoMerge on.
- `tagbot_enabled`: if tagbot is enabled.
- `authorized_authors`: list of who can submit registration, e.g `String["JuliaRegistrator"]`.
- `authorized_authors_special_jll_exceptions`: a list of users who can submit JLL packages (which have strict rules about allowed dependencies and are subject to `new_jll_*_waiting_period`s instead of `new_*_waiting_period`s).
- `additional_statuses`: list of additional commit statuses that must pass before AutoMerge will merge a PR
- `additional_check_runs`: list of additional check runs that must pass before AutoMerge will merge a PR
- `error_exit_if_automerge_not_applicable`: if `false`, AutoMerge will not error on PRs made by non-AutoMerge-authorized users
- `master_branch`: name of `master_branch`, e.g you may want to specify this to `"main"` for new GitHub repositories.
- `master_branch_is_default_branch`: if `master_branch` specified above is the default branch.
- `suggest_onepointzero`: should the AutoMerge comment include a suggestion to tag a 1.0 release for v0.x.y packages.
- `registry_deps`: list of registry dependencies, e.g your packages may depend on `General`.
- `api_url`: the registry host API URL, default is `"https://api.github.com"`.
- `check_license`: check package has a valid license, default is `false`.
- `public_registries`: If a new package registration has a UUID that matches
   that of a package already registered in one of these registries supplied here
   (and has either a different name or different URL) then an error will be thrown.
   This to prevent AutoMerge from being used for "dependency confusion"
   attacks on those registries.
- `read_only`: run in read only mode, default is `false`.

# Example

Here is an example of how `General` registry is configured

```julia
using RegistryCI
using Dates

RegistryCI.AutoMerge.run(
    merge_new_packages = ENV["MERGE_NEW_PACKAGES"] == "true",
    merge_new_versions = ENV["MERGE_NEW_VERSIONS"] == "true",
    new_package_waiting_period = Day(3),
    new_jll_package_waiting_period = Minute(20),
    new_version_waiting_period = Minute(10),
    new_jll_version_waiting_period = Minute(10),
    registry = "JuliaLang/General",
    tagbot_enabled = true,
    authorized_authors = String["JuliaRegistrator"],
    authorized_authors_special_jll_exceptions = String["jlbuild"],
    suggest_onepointzero = false,
    additional_statuses = String[],
    additional_check_runs = String[],
    check_license = true,
    public_registries = String["https://github.com/HolyLab/HolyLabRegistry"],
)
```
"""
function run(;
    env=ENV,
    cicfg::CIService=auto_detect_ci_service(; env=env),
    merge_new_packages::Bool,
    merge_new_versions::Bool,
    new_package_waiting_period,
    new_jll_package_waiting_period,
    new_version_waiting_period,
    new_jll_version_waiting_period,
    registry::String,
    #
    tagbot_enabled::Bool=false,
    #
    authorized_authors::Vector{String},
    authorized_authors_special_jll_exceptions::Vector{String},
    #
    additional_statuses::AbstractVector{<:AbstractString}=String[],
    additional_check_runs::AbstractVector{<:AbstractString}=String[],
    #
    error_exit_if_automerge_not_applicable::Bool=false,
    #
    master_branch::String="master",
    master_branch_is_default_branch::Bool=true,
    suggest_onepointzero::Bool=true,
    #
    registry_deps::Vector{<:AbstractString}=String[],
    api_url::String="https://api.github.com",
    check_license::Bool=false,
    # A list of public Julia registries (repository URLs)
    # which will be checked for UUID collisions in order to
    # mitigate the dependency confusion vulnerability. See
    # the `dependency_confusion.jl` file for details.
    public_registries::Vector{<:AbstractString}=String[],
    read_only::Bool=false,
)::Nothing
    all_statuses = deepcopy(additional_statuses)
    all_check_runs = deepcopy(additional_check_runs)
    push!(all_statuses, "automerge/decision")
    unique!(all_statuses)
    unique!(all_check_runs)
    api = GitHub.GitHubWebAPI(HTTP.URI(api_url))

    registry_head = directory_of_cloned_registry(cicfg; env=env)

    # Figure out what type of build this is
    run_pr_build = conditions_met_for_pr_build(cicfg; env=env, master_branch=master_branch)
    run_merge_build = conditions_met_for_merge_build(
        cicfg; env=env, master_branch=master_branch
    )

    if !(run_pr_build || run_merge_build)
        throw_not_automerge_applicable(
            AutoMergeWrongBuildType,
            "Build not determined to be either a PR build or a merge build. Exiting.";
            error_exit_if_automerge_not_applicable=error_exit_if_automerge_not_applicable,
        )
        return nothing
    end

    # Authentication
    key = if run_pr_build || !tagbot_enabled
        "AUTOMERGE_GITHUB_TOKEN"
    else
        "AUTOMERGE_TAGBOT_TOKEN"
    end
    auth = my_retry(() -> GitHub.authenticate(api, env[key]))
    whoami = my_retry(() -> username(api, cicfg; auth=auth))
    @info("Authenticated to GitHub as \"$(whoami)\"")
    registry_repo = my_retry(() -> GitHub.repo(api, registry; auth=auth))

    if run_pr_build
        pr_number = pull_request_number(cicfg; env=env)
        pr_head_commit_sha = current_pr_head_commit_sha(cicfg; env=env)
        pull_request_build(
            api,
            pr_number,
            pr_head_commit_sha,
            registry_repo,
            registry_head;
            auth=auth,
            authorized_authors=authorized_authors,
            authorized_authors_special_jll_exceptions=authorized_authors_special_jll_exceptions,
            error_exit_if_automerge_not_applicable=error_exit_if_automerge_not_applicable,
            master_branch=master_branch,
            master_branch_is_default_branch=master_branch_is_default_branch,
            suggest_onepointzero=suggest_onepointzero,
            whoami=whoami,
            registry_deps=registry_deps,
            check_license=check_license,
            public_registries=public_registries,
            read_only=read_only,
        )
    else
        always_assert(run_merge_build)
        cron_or_api_build(
            api,
            registry_repo;
            auth=auth,
            authorized_authors=authorized_authors,
            authorized_authors_special_jll_exceptions=authorized_authors_special_jll_exceptions,
            merge_new_packages=merge_new_packages,
            merge_new_versions=merge_new_versions,
            new_package_waiting_period=new_package_waiting_period,
            new_jll_package_waiting_period=new_jll_package_waiting_period,
            new_version_waiting_period=new_version_waiting_period,
            new_jll_version_waiting_period=new_jll_version_waiting_period,
            whoami=whoami,
            all_statuses=all_statuses,
            all_check_runs=all_check_runs,
            read_only=read_only,
        )
    end
    return nothing
end