text:
#######################################################################
# OrbitalStructure
#######################################################################
struct OrbitalStructure{T,N,O<:NTuple{N,Tuple{Vararg{NameType}}}}
    orbtype::Type{T}    # T<:Union{Number,SVector} is Hamiltonian's orbitaltype
    orbitals::O
    offsets::Tuple{Int,Vararg{Int,N}}
    flatoffsets::Tuple{Int,Vararg{Int,N}}
end

function OrbitalStructure(::Type{T}, orbs::O, lat::AbstractLattice) where {T,N,O<:NTuple{N,Any}}
    offsets = ntuple(i -> lat.unitcell.offsets[i], Val(N+1))
    offsets´ = flatoffsets(offsets, length.(orbs))
    return OrbitalStructure(T, orbs, offsets, offsets´)
end

OrbitalStructure(::Type{T}, n) where {T} = OrbitalStructure(T, ((:A,),), (0, n), (0, n))

# Equality does not need equal T, or equal orbital names
Base.:(==)(o1::OrbitalStructure, o2::OrbitalStructure) =
    length.(o1.orbitals) == length.(o2.orbitals) && o1.offsets == o2.offsets

function Base.show(io::IO, o::OrbitalStructure{T,N,O}) where {T,N,O}
    i = get(io, :indent, "")
    ioindent = IOContext(io, :indent => string(i, "  "))
    print(io,
"$(i)OrbitalStructure:
$i  Orbital Type  : $T
$i  Orbitals      : $(orbitals(o))
$i  Sublattices   : $(length(sublats(o)))
$i  Dimensions    : $(dimh(o))")
end

#######################################################################
# 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<:OrbitalStructure} # <: AbstractMatrix{M}
    lattice::LA
    harmonics::Vector{H}
    orbstruct::O
    # Enforce sorted-dns-starting-from-zero invariant onto harmonics
    function Hamiltonian{LA,L,M,A,H,O}(lattice, harmonics, orbstruct) where
        {LA<:AbstractLattice,L,M,A<:AbstractMatrix, H<:HamiltonianHarmonic{L,M,A},O<:OrbitalStructure}
        dimh = nsites(lattice)
        all(har -> size(har.h) == (dimh, dimh), harmonics) || throw(DimensionMismatch("Harmonics don't match lattice dimensions"))
        length(harmonics) > 0 && iszero(first(harmonics).dn) ||
            pushfirst!(harmonics, H(zero(SVector{L,Int}), dimh, dimh))
        sort!(harmonics, by = h -> abs.(h.dn))
        return new(lattice, harmonics, orbstruct)
    end
end

Hamiltonian(lat::LA, hs::Vector{H}, orb::O) where {LA<:AbstractLattice,L,M,A<:AbstractMatrix, H<:HamiltonianHarmonic{L,M,A},O<:OrbitalStructure} =
    Hamiltonian{LA,L,M,A,H,O}(lat, hs, orb)

Hamiltonian(lat, hs::Vector{H}, orbs::Tuple) where {L,M,H<:HamiltonianHarmonic{L,M}} =
    Hamiltonian(lat, hs, OrbitalStructure(orbitaltype(M), orbs, lat))

orbitals(h::Hamiltonian) = orbitalstructure(h).orbitals

offsets(h::Hamiltonian) = orbitalstructure(h).offsets

flatoffsets(h::Hamiltonian) = orbitalstructure(h).flatoffsets

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

Base.isequal(h1::HamiltonianHarmonic, h2::HamiltonianHarmonic) =
    h1.dn == h2.dn && h1.h == h2.h
Base.isequal(o1::OrbitalStructure, o2::OrbitalStructure) =
    o1.orbitals == o1.orbitals && o1.offsets == o2.offsets && o1.flatoffsets == o2.flatoffsets

sublat_site(siteidx, o::OrbitalStructure) = sublat_site(siteidx, o.offsets)

dimh(o::OrbitalStructure) = last(o.offsets)

sublats(o::OrbitalStructure) = 1:length(o.orbitals)

orbitals(o::OrbitalStructure) = o.orbitals

siterange(o::OrbitalStructure, sublat) = (1+o.offsets[sublat]):o.offsets[sublat+1]

# sublat offsets after flattening (without padding zeros)
flatoffsets(offsets, norbs) = _flatoffsets((0,), offsets, norbs...)
_flatoffsets(flatoffsets´::NTuple{N,Any}, offsets, n, ns...) where {N} =
    _flatoffsets((flatoffsets´..., flatoffsets´[end] + n * (offsets[N+1] - offsets[N])), offsets, ns...)
_flatoffsets(flatoffsets´, offsets) = flatoffsets´

# offset of site i after flattening
@inline flatoffset_site(i, orbstruct) = first(flatoffsetorbs_site(i, orbstruct))

function flatoffsetorbs_site(i, orbstruct)
    s = sublat_site(i, orbstruct)
    N = length(orbitals(orbstruct)[s])
    offset = orbstruct.offsets[s]
    offset´ = orbstruct.flatoffsets[s]
    Δi = i - offset
    i´ = offset´ + (Δi - 1) * N
    return i´, N
end

function flatindices(i, orbstruct)
    o, N = flatoffsetorbs_site(i, orbstruct)
    return o+1:o+N
end

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) = displayorbitals(orbitalstructure(h))
displayorbitals(o::OrbitalStructure) =
    replace(replace(string(orbitals(o)), "Symbol(\"" => ":"), "\")" => "")

function Base.show(io::IO, 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     : $(coordination(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"

#######################################################################
# flatten
#######################################################################
"""
    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.

    flatten(k::Ket)

Flattens a multiorbital `Ket` to have a scalar `eltype`, instead of `SVector`.

    flatten(s::Subspace)

Rebuild `s` by flattening its basis to have a scalar eltype.

    x |> flatten()

Curried form equivalent to `flatten(x)` or `x |> flatten` (included for consistency with
the rest of the API).

    flatten(x, o::OrbitalStructure)

Flatten object x, if applicable, using the orbital structure `o`, as obtained from a
Hamiltonian `h` with `orbitalstructure(h)`. `x` here is typically an `AbstractArray` of
non-scalar `eltype`.

# Examples

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

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

# See also
    `unflatten`
"""
flatten() = h -> flatten(h)

function flatten(h::Hamiltonian)
    isflat(h) && return copy(h)
    harmonics´ = [flatten(har, orbitalstructure(h)) for har in h.harmonics]
    lattice´ = flatten(h.lattice, orbitalstructure(h))
    orbs´ = (_ -> (:flat, )).(orbitals(h))
    return Hamiltonian(lattice´, harmonics´, orbs´)
end

# special method: if already flat, don't make a copy
maybeflatten(h, args...) = isflat(h) ? h : flatten(h, args...)

isflat(h::Hamiltonian) = all(isequal(1), norbitals(h))
isflat(m::AbstractArray{<:Number}) = true
isflat(m) = false

flatten(h::HamiltonianHarmonic, orbstruct) =
    HamiltonianHarmonic(h.dn, flatten(h.h, orbstruct))

function flatten(o::OrbitalStructure{T}) where {T}
    T´ = eltype(T)
    orbs = (_ -> (:flat, )).(orbitals(o))
    offsets = o.flatoffsets
    flatoffsets = o.flatoffsets
    return OrbitalStructure(T´, orbs, offsets, flatoffsets)
end

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

function flatten(unitcell::Unitcell, orbstruct) # orbs::NTuple{S,Any}) where {S}
    norbs = length.(orbitals(orbstruct))
    nsublats = length(sublats(orbstruct))
    offsets´ = orbstruct.flatoffsets
    ns´ = last(offsets´)
    sites´ = similar(unitcell.sites, ns´)
    i = 1
    for sl in 1:nsublats, site in sitepositions(unitcell, sl), rep in 1:norbs[sl]
        sites´[i] = site
        i += 1
    end
    names´ = unitcell.names
    unitcell´ = Unitcell(sites´, names´, offsets´)
    return unitcell´
end

flatten(src::AbstractArray{<:Number}, orbstruct) = src
flatten(src::AbstractArray{T}, orbstruct, ::Type{T}) where {T<:Number} = src
flatten(src::AbstractArray{T}, orbstruct, ::Type{T´}) where {T<:Number,T´<:Number} = T´.(src)

function flatten(src::SparseMatrixCSC{<:SMatrix{N,N,T}}, orbstruct, ::Type{T´} = T) where {N,T,T´}
    norbs = length.(orbitals(orbstruct))
    offsets´ = orbstruct.flatoffsets
    dim´ = last(offsets´)

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

    for col in 1:size(src, 2)
        scol = sublat_site(col, orbstruct)
        for j in 1:norbs[scol]
            for p in nzrange(src, col)
                row = rowvals(src)[p]
                srow = sublat_site(row, orbstruct)
                rowoffset´ = flatoffset_site(row, orbstruct)
                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::StridedMatrix{<:SMatrix{N,N,T}}, orbstruct, ::Type{T´} = T) where {N,T,T´}
    dim´ = last(orbstruct.flatoffsets)
    dst = similar(src, T´, dim´, dim´)
    return flatten!(dst, src, orbstruct)
end

function flatten(src::StridedMatrix{<:SVector{N,T}}, orbstruct, ::Type{T´} = T) where {N,T,T´}
    dim´ = last(orbstruct.flatoffsets)
    dst = similar(src, T´, dim´, size(src, 2))
    return flatten!(dst, src, orbstruct)
end

function flatten!(dst, src::StridedMatrix{<:SMatrix{N,N,T}}, orbstruct, ::Type{T´} = T) where {N,T,T´}
    size(src) == (dimh(orbstruct), dimh(orbstruct)) || throw(ArgumentError("Dimension mismatch between source array and orbital structure"))
    norbs = length.(orbitals(orbstruct))
    for col in 1:size(src, 2), row in 1:size(src, 1)
        srow, scol = sublat_site(row, orbstruct), sublat_site(col, orbstruct)
        nrow, ncol = norbs[srow], norbs[scol]
        val = src[row, col]
        rowoffset´ = flatoffset_site(row, orbstruct)
        coloffset´ = flatoffset_site(col, orbstruct)
        for j in 1:ncol, i in 1:nrow
            dst[rowoffset´ + i, coloffset´ + j] = val[i, j]
        end
    end
    return dst
end

# for Subspace bases
function flatten!(dst, src::StridedMatrix{<:SVector{N,T}}, orbstruct, ::Type{T´} = T) where {N,T,T´}
    norbs = length.(orbitals(orbstruct))
    for col in 1:size(src, 2), row in 1:size(src, 1)
        srow = sublat_site(row, orbstruct)
        nrow = norbs[srow]
        val = src[row, col]
        rowoffset´ = flatoffset_site(row, orbstruct)
        for i in 1:nrow
            dst[rowoffset´ + i, col] = val[i]
        end
    end
    return dst
end

function flatten_sparse_copy!(dst, src, o::OrbitalStructure)
    fill!(dst, zero(eltype(dst)))
    norbs = length.(orbitals(o))
    coloffset = 0
    for s´ in sublats(o)
        N´ = norbs[s´]
        for col in siterange(o, s´)
            for p in nzrange(src, col)
                val = nonzeros(src)[p]
                row = rowvals(src)[p]
                rowoffset, M´ = flatoffsetorbs_site(row, o)
                for j in 1:N´, i in 1:M´
                    dst[i + rowoffset, j + coloffset] = val[i, j]
                end
            end
            coloffset += N´
        end
    end
    return dst
end

# Specialized flattening copy! and muladd!
function flatten_sparse_muladd!(dst, src, o::OrbitalStructure, α = I)
    norbs = length.(orbitals(o))
    coloffset = 0
    for s´ in sublats(o)
        N´ = norbs[s´]
        for col in siterange(o, s´)
            for p in nzrange(src, col)
                val = α * nonzeros(src)[p]
                row = rowvals(src)[p]
                rowoffset, M´ = flatoffsetorbs_site(row, o)
                for j in 1:N´, i in 1:M´
                    dst[i + rowoffset, j + coloffset] += val[i, j]
                end
            end
            coloffset += N´
        end
    end
    return dst
end

function flatten_dense_muladd!(dst, src, o::OrbitalStructure, α = I)
    norbs = length.(orbitals(o))
    coloffset = 0
    for s´ in sublats(o)
        N´ = norbs[s´]
        for col in siterange(o, s´)
            rowoffset = 0
            for s in sublats(o)
                M´ = norbs[s]
                for row in siterange(o, s)
                    val = α * src[row, col]
                    for j in 1:N´, i in 1:M´
                        dst[i + rowoffset, j + coloffset] += val[i, j]
                    end
                    rowoffset += M´
                end
            end
            coloffset += N´
        end
    end
    return dst
end

function flatten_dense_copy!(dst, src, o::OrbitalStructure)
    fill!(dst, zero(eltype(dst)))
    return flatten_dense_muladd!(dst, src, o, I)
end

## unflatten ##
"""
    unflatten(x, o::OrbitalStructure)

Rebuild object `x` performing the inverse of `flatten(x)` or `flatten(x, o)`. The target `o`
is required.

    x |> unflatten(o::OrbitalStructure)

Curried form equivalent to `unflatten(x, o)` (included for consistency with the
rest of the API).

# Examples
```jldoctest
julia> h = LP.honeycomb() |> hamiltonian(hopping(I), orbitals = (:up,:down)) |> unitcell;

julia> psi = spectrum(h)[around = -1]
Subspace{0}: eigenenergy subspace on a 0D manifold
  Energy       : -0.9999999999999989
  Degeneracy   : 2
  Bloch/params : Float64[]
  Basis eltype : SVector{2, ComplexF64}

julia> psiflat = flatten(psi)
Subspace{0}: eigenenergy subspace on a 0D manifold
  Energy       : -0.9999999999999989
  Degeneracy   : 2
  Bloch/params : Float64[]
  Basis eltype : ComplexF64

julia> unflatten(psiflat, orbitalstructure(psi))
Subspace{0}: eigenenergy subspace on a 0D manifold
  Energy       : -0.9999999999999989
  Degeneracy   : 2
  Bloch/params : Float64[]
  Basis eltype : SVector{2, ComplexF64}

julia> k = ket(ketmodel(1, sublats = :up), h) |> flatten
Ket{ComplexF64}: ket with a 4 × 1 amplitude matrix
  OrbitalStructure:
    Orbital Type  : ComplexF64
    Orbitals      : ((:flat,), (:flat,))
    Sublattices   : 2
    Dimensions    : 4

julia> unflatten(k, orbitalstructure(h))
Ket{SVector{2, ComplexF64}}: ket with a 2 × 1 amplitude matrix
  OrbitalStructure:
    Orbital Type  : SVector{2, ComplexF64}
    Orbitals      : ((:up, :down), (:up, :down))
    Sublattices   : 2
    Dimensions    : 2
```

# See also
    `flatten`, `orbitalstructure`
"""
unflatten(o::OrbitalStructure) = x -> unflatten(x, o)

# Pending implementation for Hamiltonian, Lattice, Unitcell

unflatten_orbitals(v::AbstractVector, o::OrbitalStructure) =
    isunflat_orbitals(v, o) ? copy(v) : unflatten!(similar(v, orbitaltype(o), dimh(o)), v, o)
unflatten_orbitals(m::AbstractMatrix, o::OrbitalStructure) =
    isunflat_orbitals(m, o) ? copy(m) : unflatten!(similar(m, orbitaltype(o), dimh(o), size(m, 2)), m, o)
unflatten_blocks(m::AbstractMatrix, o::OrbitalStructure) =
    isunflat_blocks(m, o) ? copy(m) : unflatten!(similar(m, blocktype(o), dimh(o), dimh(o)), m, o)

maybe_unflatten_orbitals(x, o) = isunflat_orbitals(x, o) ? x : unflatten_orbitals(x, o)
maybe_unflatten_blocks(x, o) = isunflat_blocks(x, o) ? x : unflatten_blocks(x, o)

isunflat_orbitals(m, o) = orbitaltype(o) == eltype(m) && size(m, 1) == dimh(o)
isunflat_blocks(m, o) = blocktype(o) == eltype(m) && size(m, 1) == dimh(o)

function unflatten!(v::AbstractArray, vflat::AbstractArray, o::OrbitalStructure)
    dimflat = last(o.flatoffsets)
    check_unflatten_dst_dims(v, dimh(o))
    check_unflatten_src_dims(vflat, dimflat)
    check_unflatten_eltypes(v, o)
    v = _unflatten!(v, vflat, o)
    return v
end

# unflatten into SMatrix blocks
function _unflatten!(v::AbstractArray{T}, vflat::AbstractArray, o::OrbitalStructure) where {T<:SMatrix}
    norbs = length.(orbitals(o))
    flatoffsets = o.flatoffsets
    col = 0
    for scol in sublats(o)
        M = colstride(norbs, scol, T)
        for j in flatoffsets[scol]+1:M:flatoffsets[scol+1]
            col +=1
            row = 0
            for srow in sublats(o)
                N = rowstride(norbs, srow)
                for i in flatoffsets[srow]+1:N:flatoffsets[srow+1]
                    row += 1
                    val = view(vflat, i:i+N-1, j:j+M-1)
                    v[row, col] = padtotype(val, T)
                end
            end
        end
    end
    return v
end

# unflatten into SVector orbitals
function _unflatten!(v::AbstractArray{T}, vflat::AbstractArray, o::OrbitalStructure) where {T<:SVector}
    norbs = length.(orbitals(o))
    flatoffsets = o.flatoffsets
    col = 0
    for col in 1:size(v, 2)
        row = 0
        for srow in sublats(o)
            N = rowstride(norbs, srow)
            for i in flatoffsets[srow]+1:N:flatoffsets[srow+1]
                row += 1
                val = view(vflat, i:i+N-1, col:col)
                v[row, col] = padtotype(val, T)
            end
        end
    end
    return v
end

rowstride(norbs, s) = norbs[s]
colstride(norbs, s, ::Type{<:SVector}) = 1
colstride(norbs, s, ::Type{<:SMatrix}) = rowstride(norbs, s)

# dest v should be a vector or matrix such that H*v is possible
check_unflatten_dst_dims(v, dimh) =
    size(v, 1) == dimh ||
        throw(ArgumentError("Dimension of destination array is inconsistent with orbital structure"))

# dest v should be a square matrix like the Hamiltonian
check_unflatten_dst_dims(v::AbstractArray{<:SMatrix}, dimh) =
    size(v, 1) == dimh && size(v, 2) == dimh ||
        throw(ArgumentError("Dimension of destination array is inconsistent with orbital structure"))

check_unflatten_src_dims(vflat, dimflat) =
    size(vflat, 1) == dimflat ||
        throw(ArgumentError("Dimension of source array is inconsistent with orbital structure"))

check_unflatten_eltypes(::AbstractArray{T}, o::OrbitalStructure) where {T} =
    T === orbitaltype(o) || T === blocktype(o) ||
        throw(ArgumentError("Eltype of desination array is inconsistent with orbital structure"))

# unflatten_orbitals_or_reinterpret: call unflatten_orbitals but only if we cannot unflatten via reinterpret
unflatten_orbitals_or_reinterpret(vflat, ::Missing) = vflat
# source is already of the correct orbitaltype(h)
function unflatten_orbitals_or_reinterpret(vflat::AbstractArray{<:SArray{N}}, o::OrbitalStructure{<:SArray{N}}) where {N}
    check_unflatten_dst_dims(vflat, dimh(o))
    return vflat
end

function unflatten_orbitals_or_reinterpret(vflat::AbstractArray{<:Number}, o::OrbitalStructure{<:Number})
    check_unflatten_dst_dims(vflat, dimh(o))
    return vflat
end

# source can be reinterpreted, because the number of orbitals is the same M for all N sublattices
unflatten_orbitals_or_reinterpret(vflat::AbstractArray{T}, o::OrbitalStructure{S,N,<:NTuple{N,NTuple{M}}}) where {T,S,N,M} =
    reinterpret(SVector{M,T}, vflat)
# otherwise call unflatten_orbitals
unflatten_orbitals_or_reinterpret(vflat, o) = unflatten_orbitals(vflat, o)

#######################################################################
# similarmatrix
#######################################################################
"""
    similarmatrix(h::Hamiltonian)

Create an uninitialized matrix of the same type and size of the Hamiltonian's matrix.

    similarmatrix(h::Hamiltonian, T::Type{<:AbstractMatrix})

Make the matrix of type `B<:T`. Can be used to specify a different eltype from `h`'s, (e.g.
`T=SparseMatrixCSC{Float64}` with a multiorbital `h`)

    similarmatrix(h::Hamiltonian, flatten)

Create an unitialized matrix of the same type as a flattened version of `h`'s, i.e. with a
scalar eltype as in the example above.

    similarmatrix(h::Hamiltonian, T::AbstractDiagonalizeMethod)

Adapts the type of the matrix (e.g. dense/sparse) to the specified `method`

    similarmatrix(x::Union{ParametricHamiltonian, GreensFunction}, ...)

Equivalent to the above, but adapted to the more general type of `x`.

# Examples

```jldoctest
julia> h = LatticePresets.honeycomb() |> hamiltonian(hopping(I), orbitals = Val(2));

julia> similarmatrix(h) |> summary
"2×2 SparseMatrixCSC{SMatrix{2, 2, ComplexF64, 4}, Int64}"

julia> similarmatrix(h, Matrix{Int}) |> summary
"4×4 Matrix{Int64}"

julia> similarmatrix(h, flatten) |> summary
"4×4 SparseMatrixCSC{ComplexF64, Int64}"

julia> similarmatrix(h, LinearAlgebraPackage()) |> summary
"4×4 Matrix{ComplexF64}"
```

# See also
    `bloch!`
"""
similarmatrix(h, dest_type = missing) = _similarmatrix(dest_type, matrixtype(h), h)

_similarmatrix(::Missing, src_type, h) =
    similar_merged(h.harmonics)
_similarmatrix(::Type{A}, ::Type{A´}, h) where {T<:Number,A<:AbstractSparseMatrix{T},T´<:Number,A´<:AbstractSparseMatrix{T´}} =
    similar_merged(h.harmonics, T)
_similarmatrix(::Type{A}, ::Type{A´}, h) where {N,T<:SMatrix{N,N},A<:AbstractSparseMatrix{T},T´<:SMatrix{N,N},A´<:AbstractSparseMatrix{T´}} =
    similar_merged(h.harmonics, T)
_similarmatrix(::Type{A}, ::Type{A´}, h) where {N,T<:Number,A<:AbstractSparseMatrix{T},T´<:SMatrix{N,N},A´<:AbstractSparseMatrix{T´}} =
    flatten(similar_merged(h.harmonics), orbitalstructure(h), T)
_similarmatrix(::Type{A}, ::Type{A´}, h) where {T<:Number,A<:Matrix{T},T´<:Number,A´<:AbstractMatrix{T´}} =
    similar(A, size(h))
_similarmatrix(::Type{A}, ::Type{A´}, h) where {N,T<:SMatrix{N,N},A<:Matrix{T},T´<:SMatrix{N,N},A´<:AbstractMatrix{T´}} =
    similar(A, size(h))
_similarmatrix(::Type{A}, ::Type{A´}, h) where {N,T<:Number,A<:Matrix{T},T´<:SMatrix{N,N},A´<:AbstractMatrix{T´}} =
    similar(A, flatsize(h))

_similarmatrix(::typeof(flatten), ::Type{A´}, h) where {N,T,S<:SMatrix{N,N,T},A´<:AbstractSparseMatrix{S}} =
    _similarmatrix(AbstractSparseMatrix{T}, A´, h)
_similarmatrix(::typeof(flatten), ::Type{A´}, h) where {N,T,S<:SMatrix{N,N,T},A´<:StridedMatrix{S}} =
    _similarmatrix(Matrix{T}, A´, h)
_similarmatrix(::typeof(flatten), ::Type{A´}, h) where {T<:Number,A´<:AbstractArray{T}} =
    _similarmatrix(A´, A´, h)
_similarmatrix(::Type{SparseMatrixCSC}, ::Type{A´}, h) where {T´,A´<:AbstractMatrix{T´}} =
    _similarmatrix(SparseMatrixCSC{T´}, A´, h)
_similarmatrix(::Type{Matrix}, ::Type{A´}, h) where {T´,A´<:AbstractMatrix{T´}} =
    _similarmatrix(Matrix{T´}, A´, h)

_similarmatrix(dest_type, src_type, h) = throw(ArgumentError("Unexpected `blochtype` ($src_type => $dest_type)"))

# ensure we have stored entries for all harmonics in the sparse case
similar_merged(hs::AbstractVector{H}, ::Type{M´} = M) where {M´,M,H<:HamiltonianHarmonic{<:Any,M,<:AbstractSparseMatrix}} =
    _similar_merged_sparse(hs, M´)
similar_merged(hs::AbstractVector{H}, ::Type{M´} = M) where {M´,M,H<:HamiltonianHarmonic{<:Any,M,<:StridedMatrix}} =
    similar(first(hs).h, M´)

function _similar_merged_sparse(hs::Vector{<:HamiltonianHarmonic{L,M}},::Type{M´} = M) where {L,M,M´}
    h0 = first(hs)
    n, m = size(h0.h)
    iszero(h0.dn) || throw(ArgumentError("First Hamiltonian harmonic is not the fundamental"))
    nh = length(hs)
    builder = SparseMatrixBuilder{M´}(n, m)
    for col in 1:m
        for i in eachindex(hs)
            h = hs[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)
    return ho
end
# IDEA: could sum and subtract all harmonics instead: sum(h->h.h, hs)
# Tested, it is slower

#######################################################################
# Hamiltonian API
#######################################################################
"""
    hamiltonian(lat, model; orbitals, orbtype)

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

    lat |> hamiltonian(model; kw...)

Curried form of `hamiltonian` equivalent to `hamiltonian(lat, model; kw...)`.

# Keywords

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 Hamiltonian element type is `orbtype`. Otherwise it is `SMatrix{N,N,orbtype}` blocks,
padded with the necessary zeros as required. Keyword `orbtype` is `Complex{T}` by default,
where `T` is the number type of `lat`.

# Indexing

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

A slicing syntax `h[rows, cols]` (without specifying `dn`) is also available, that creates a
special `hs::Slice{<:Hamiltonian}` object that represents a slice of the Hamiltonian matrix
restricted to `rows` and `cols`. Here, `rows` and `cols` are collections of site indices, or
alternatively `SiteSelector`s (see `siteselector`s for details). If `rows::Integer` or
`cols::Integer`, they will be converted to a single-element range (to preserve always a
matrix-like slice, unlike for `AbstractArray` indexing). Slices support `bloch` and `bloch!`
to produce the corresponding matrices, and can also be indexed as `hs[dn::Tuple]` that
produces the equivalent to `h[dn][rows, cols]`.

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} : 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 (ComplexF64)
  Onsites          : 0
  Hoppings         : 6
  Coordination     : 3.0

julia> push!(h, (3,3)) # Adding a new Hamiltonian harmonic (if not already present)
Hamiltonian{<:Lattice} : 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 (ComplexF64)
  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 SparseMatrixCSC{SMatrix{2, 2, ComplexF64, 4}, Int64} with 1 stored entry:
 [1.0+0.0im 2.0+0.0im; 2.0+0.0im 1.0+0.0im]                      ⋅
                     ⋅                                           ⋅

julia> h[(3,3)][[1,2],[1,2]] .= Ref(@SMatrix[1 2; 2 1])
2×2 SparseMatrixCSC{SMatrix{2, 2, ComplexF64, 4}, Int64} with 4 stored entries:
 [1.0+0.0im 2.0+0.0im; 2.0+0.0im 1.0+0.0im]  [1.0+0.0im 2.0+0.0im; 2.0+0.0im 1.0+0.0im]
 [1.0+0.0im 2.0+0.0im; 2.0+0.0im 1.0+0.0im]  [1.0+0.0im 2.0+0.0im; 2.0+0.0im 1.0+0.0im]

 julia> h = unitcell(h); h[]
2×2 SparseMatrixCSC{SMatrix{2, 2, ComplexF64, 4}, Int64} with 2 stored entries:
                     ⋅                       [1.0+0.0im 2.0+0.0im; 3.0+0.0im 4.0+0.0im]
 [1.0+0.0im 2.0+0.0im; 3.0+0.0im 4.0+0.0im]                      ⋅

```

# See also
    `onsite`, `hopping`, `bloch`, `bloch!`
"""
hamiltonian(lat::AbstractLattice, 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, m::TightbindingModel; orbtype::Type = Complex{numbertype(lat)}, kw...) =
    hamiltonian_sparse(blocktype(orbs, orbtype), lat, orbs, m; kw...)

hamiltonian(t::TightbindingModel...; kw...) = lat -> hamiltonian(lat, t...; kw...)

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`"))

"""
    orbitalstructure(x::Union{Hamiltonian,ParametricHamiltonian})
    orbitalstructure(x::Subspace)

Return an `OrbitalStructure` containing information about the orbital structure of `x`

# Examples

```jldoctest
julia> sp = spectrum(LP.honeycomb() |> hamiltonian(hopping(I), orbitals = (:up,:down)) |> unitcell);

julia> sp[around = -1] |> orbitalstructure
OrbitalStructure:
  Orbital Type  : SVector{2, ComplexF64}
  Orbitals      : ((:up, :down), (:up, :down))
  Sublattices   : 2
  Dimensions    : 2
```
"""
orbitalstructure(h::Hamiltonian) = h.orbstruct

"""
    dims(lh::Union{Hamiltonian,AbstractLattice}) -> (E, L)

Return a tuple `(E, L)` of the embedding `E` and lattice dimensions `L` of `AbstractLattice`
or `Hamiltonian` `lh`
"""
dims(lat::AbstractLattice{E,L}) where {E,L} = E, L
dims(h::Hamiltonian) = dims(h.lattice)

"""
    sitepositions(lat::AbstractLattice; kw...)
    sitepositions(h::Hamiltonian; kw...)

Build a generator of the positions of sites in the lattice unitcell. Only sites specified
by `siteselector(kw...)` are selected, see `siteselector` for details.

"""
sitepositions(lat::AbstractLattice; kw...) = sitepositions(lat, siteselector(;kw...))
sitepositions(h::Hamiltonian; kw...) = sitepositions(h.lattice, siteselector(;kw...))

"""
    siteindices(lat::AbstractLattice; kw...)
    siteindices(h::Hamiltonian; kw...)
    siteindices(lat_or_h, s::SiteSelector)

Build a generator of the unique indices of sites in the lattice unitcell. Only sites
specified by `siteselector(kw...)` are selected, see `siteselector` for details.

"""
siteindices(lat::AbstractLattice; kw...) = siteindices(lat, siteselector(;kw...))
siteindices(h::Hamiltonian; kw...) = siteindices(h.lattice, siteselector(;kw...))
siteindices(h::Hamiltonian, s::SiteSelector) = siteindices(h.lattice, s)

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

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

latdim(h::Hamiltonian) = last(dims(h.lattice))

matrixtype(::Hamiltonian{LA,L,M,A}) where {LA,L,M,A} = A

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

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

# 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
blocktype(o::OrbitalStructure) = _blocktype(orbitaltype(o))

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

orbitaltype(h::Hamiltonian) = orbitaltype(orbitalstructure(h))
orbitaltype(o::OrbitalStructure) = o.orbtype
orbitaltype(::Type{M}) where {M<:Number} = M
orbitaltype(::Type{S}) where {N,T,S<:SMatrix{N,N,T}} = SVector{N,T}

# find SVector type that can hold all orbital amplitudes in any lattice sites
orbitaltype(orbs, type::Type{Tv}) where {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

coordination(ham) = nhoppings(ham) / nsites(ham)

function nhoppings(ham)
    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)
    count = 0
    for h in ham.harmonics
        iszero(h.dn) && (count += _nnzdiag(h.h))
    end
    return count
end

_nnz(h::AbstractSparseMatrix) = count(!iszero, nonzeros(h)) # Does not include stored zeros
_nnz(h::StridedMatrix) = count(!iszero, h)

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

_nnzdiag(s::Matrix) = count(!iszero, s[i,i] for i in 1:minimum(size(s)))

check_orbital_consistency(::Hamiltonian{<:Any,<:Any,<:Number}) = nothing

function check_orbital_consistency(h::Hamiltonian)
    slats = sublats(h.lattice)
    norbs = norbitals(h)
    M = maximum(norbs)
    all(norbs .== M) && return nothing
    for srow in slats, scol in slats
        nrow, ncol = norbs[srow], norbs[scol]
        nrow == ncol == M && continue
        rows = siterange(h.lattice, srow)
        cols = siterange(h.lattice, scol)
        for hh in h.harmonics
            hrows = rowvals(hh.h)
            nz = nonzeros(hh.h)
            for col in cols, ptr in nzrange(hh.h, col)
                row = hrows[ptr]
                row in rows || continue
                z = nz[ptr]
                is_consistent_orbs(z, nrow, ncol) || _throw_inconsistent_orbs()
            end
        end
    end
    return nothing
end

function is_consistent_orbs(z::SMatrix{N,N}, rows, cols) where {N}
    for j in cols + 1:N, i in 1:N
        iszero(z[i,j]) || return false
    end
    for j in 1:N-cols, i in rows+1:N
        iszero(z[i,j]) || return false
    end
    return true
end

@noinline _throw_inconsistent_orbs() = throw(ArgumentError("Internal error: orbital structure not correctly encoded into Hamiltonian harmonics"))

bravais(h::Hamiltonian) = bravais(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.(orbitals(h))

flatsize(h::Hamiltonian, n) = first(flatsize(h)) # h is always square

function flatsize(h::Hamiltonian)
    n = sum(sublatlengths(h.lattice) .* length.(orbitals(h)))
    return (n, n)
end

expand_supercell_mask(h::Hamiltonian{<:Superlattice}) =
    Hamiltonian(expand_supercell_mask(h.lattice), h.harmonics, orbitals(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)

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

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

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

Base.isequal(h1::Hamiltonian, h2::Hamiltonian) =
    isequal(h1.lattice, h2.lattice) && isequal(h1.harmonics, h2.harmonics) &&
    isequal(orbitalstructure(h1), orbitalstructure(h2))

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

Base.parent(h::Hamiltonian) = h

# Iterators #

function nonzero_indices(h::Hamiltonian, rowrange = 1:size(h, 1), colrange = 1:size(h, 2))
    rowrange´ = rclamp(rowrange, 1:size(h, 1))
    colrange´ = rclamp(colrange, 1:size(h, 2))
    gen = ((har.dn, rowvals(har.h)[ptr], col)
                for har in h.harmonics
                for col in colrange´
                for ptr in nzrange_inrows(har.h, col, rowrange´)
                if !iszero(nonzeros(har.h)[ptr]))
    return gen
end

function nonzero_indices(har::HamiltonianHarmonic, rowrange = 1:size(har, 1), colrange = 1:size(har, 2))
    rowrange´ = rclamp(rowrange, 1:size(har, 1))
    colrange´ = rclamp(colrange, 1:size(har, 2))
    gen = ((rowvals(har.h)[ptr], col)
                for col in colrange´
                for ptr in nzrange_inrows(har.h, col, rowrange´)
                if !iszero(nonzeros(har.h)[ptr]))
    return gen
end

function nzrange_inrows(h, col, rowrange)
    ptrs = nzrange(h, col)
    rows = rowvals(h)
    ptrmin = first(ptrs)
    ptrmax = last(ptrs)

    for p in ptrs
        rows[p] in rowrange && break
        ptrmin = p + 1
    end

    if ptrmin < ptrmax
        for p in ptrmax:-1:ptrmin
            ptrmax = p
            rows[p] in rowrange && break
        end
    end

    return ptrmin:ptrmax
end

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}, dn::SVector{L,Int}) where {L}
    get_or_push!(h.harmonics, dn, size(h))
    return h
end

function get_or_push!(harmonics::Vector{H}, dn::SVector{L,Int}, dims) where {L,M,A,H<:HamiltonianHarmonic{L,M,A}}
    for hh in harmonics
        hh.dn == dn && return hh
    end
    hh = HamiltonianHarmonic{L,M,A}(dn, dims...)
    push!(harmonics, hh)
    return hh
end

# Indexing #

Base.getindex(h::Hamiltonian{<:Any,L}) where {L} = h[zero(SVector{L,Int})]
Base.getindex(h::Hamiltonian, dn::NTuple{L}) where {L} = getindex(h, SVector{L,Int}(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.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

## Boolean masking
"""
    &(h1::Hamiltonian{<:Superlattice}, h2::Hamiltonian{<:Superlattice})

Construct a new `Hamiltonian{<:Superlattice}` using an `and` boolean mask, i.e. with a
supercell that contains cells that are both in the supercell of `h1` and `h2`

    &(s1::Superlattice, s2::Superlattice}

Equivalent of the above for `Superlattice`s

# See also
    `|`, `xor`
"""
(Base.:&)(s1::Hamiltonian{<:Superlattice}, s2::Hamiltonian{<:Superlattice}) =
    boolean_mask_hamiltonian(Base.:&, s1, s2)

"""
    |(h1::Hamiltonian{<:Superlattice}, h2::Hamiltonian{<:Superlattice})

Construct a new `Hamiltonian{<:Superlattice}` using an `or` boolean mask, i.e. with a
supercell that contains cells that are either in the supercell of `h1` or `h2`

    |(s1::Superlattice, s2::Superlattice}

Equivalent of the above for `Superlattice`s

# See also
    `&`, `xor`
"""
(Base.:|)(s1::Hamiltonian{<:Superlattice}, s2::Hamiltonian{<:Superlattice}) =
    boolean_mask_hamiltonian(Base.:|, s1, s2)

"""
    xor(h1::Hamiltonian{<:Superlattice}, h2::Hamiltonian{<:Superlattice})

Construct a new `Hamiltonian{<:Superlattice}` using a `xor` boolean mask, i.e. with a
supercell that contains cells that are either in the supercell of `h1` or `h2` but not in
both

    xor(s1::Superlattice, s2::Superlattice}

Equivalent of the above for `Superlattice`s

# See also
    `&`, `|`
"""
(Base.xor)(s1::Hamiltonian{<:Superlattice}, s2::Hamiltonian{<:Superlattice}) =
    boolean_mask_hamiltonian(Base.xor, s1, s2)

function boolean_mask_hamiltonian(f, s1::Hamiltonian{<:Superlattice}, s2::Hamiltonian{<:Superlattice})
    check_compatible_hsuper(s1, s2)
    return Hamiltonian(f(s1.lattice, s2.lattice), s1.harmonics, orbitals(s1))
end

function check_compatible_hsuper(s1, s2)
    compatible = isequal(s1.harmonics, s2.harmonics) && isequal(s1.orbstruct, s2.orbstruct)
    compatible || throw(ArgumentError("Hamiltonians are incompatible for boolean masking"))
    return nothing
end

## Hamiltonian algebra

function Base.literal_pow(::typeof(^), h::Hamiltonian, p::Val{P}) where {P}
    P > 0 || throw(ArgumentError("Only positive powers of Hamiltonians are supported"))
    P == 1 && return copy(h)
    hhs = similar(h.harmonics, 0)
    nh = length(h.harmonics)
    dims = size(h)
    hhiter = CartesianIndices(ntuple(i -> 1:nh, p))
    for Is in hhiter
        is = Tuple(Is)
        dn = sum(i -> h.harmonics[i].dn, is)
        hh = get_or_push!(hhs, dn, dims)
        hh.h .+= prod((i -> h.harmonics[i].h).(is))
    end
    return Hamiltonian(h.lattice, hhs, orbitalstructure(h))
end

function Base.:*(h1::Hamiltonian, h2::Hamiltonian) where {P}
    check_compatible_hamiltonians(h1, h2)
    hhs = similar(h1.harmonics, 0)
    dims = size(h1)
    for hh1 in h1.harmonics, hh2 in h2.harmonics
        dn = hh1.dn + hh2.dn
        hh = get_or_push!(hhs, dn, dims)
        mul!(hh.h, hh1.h, hh2.h, 1, 1)
    end
    return Hamiltonian(h1.lattice, hhs, orbitalstructure(h1))
end

Base.:*(h::Hamiltonian, p::Number) = p * h

function Base.:*(p::Number, h::Hamiltonian)
    hhs = copy.(h.harmonics)
    for hh in hhs
        hh.h .*= p
    end
    return Hamiltonian(h.lattice, hhs, orbitalstructure(h))
end

Base.:-(h1::Hamiltonian, h2) = h1 + (-h2)
Base.:-(h::Hamiltonian) = (-1) * h

function Base.:+(h1::Hamiltonian, h2::Hamiltonian)
    check_compatible_hamiltonians(h1, h2)
    hhs = copy.(h1.harmonics)
    dims = size(h1)
    for hh2 in h2.harmonics
        hh = get_or_push!(hhs, hh2.dn, dims)
        hh.h .+= hh2.h
    end
    return Hamiltonian(h1.lattice, hhs, orbitalstructure(h1))
end

Base.:+(id::UniformScaling, h::Hamiltonian) = h + id

function Base.:+(h::Hamiltonian, id::UniformScaling)
    hhs = copy.(h.harmonics)
    M = blocktype(h)
    mat = first(hhs).h
    shift_diagonal!(mat, blocktype(h), h, id)
    return Hamiltonian(h.lattice, hhs, orbitalstructure(h))
end

function shift_diagonal!(mat, ::Type{<:Number}, h, id::UniformScaling)
    mat .+= id
    return mat
end

function shift_diagonal!(mat, ::Type{S}, h, id::UniformScaling) where {N,S<:SMatrix{N,N}}
    for s in sublats(h.lattice)
        shift = id.λ * padprojector(S, orbitals(h)[s])
        for i in siterange(h.lattice, s)
            mat[i,i] += shift
        end
    end
    return mat
end

check_compatible_hamiltonians(h1, h2) =
    isequal(h1.lattice, h2.lattice) && isequal(orbitalstructure(h1), orbitalstructure(h2)) && size(h1) == size(h2) ||
        throw(ArgumentError("Cannot combine Hamiltonians with different lattices, dimensions or orbitals"))
#######################################################################
# 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

Base.eltype(b::IJVBuilder{L,M}) where {L,M} = M

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}
    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)
end

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

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

function applyterm!(builder::IJVBuilder{L}, term::OnsiteTerm) where {L}
    lat = builder.lat
    dn0 = zero(SVector{L,Int})
    ijv = builder[dn0]
    allpos = allsitepositions(lat)
    rsel = resolve(term.selector, lat)
    for s in sublats(rsel), i in siteindices(rsel, s)
        r = allpos[i]
        v = to_blocktype(term(r, r), eltype(builder), builder.orbs[s], builder.orbs[s])
        push!(ijv, (i, i, v))
    end
    return nothing
end

function applyterm!(builder::IJVBuilder{L}, term::HoppingTerm) where {L}
    lat = builder.lat
    rsel = resolve(term.selector, lat)
    L > 0 && checkinfinite(rsel)
    allpos = allsitepositions(lat)
    for (s2, s1) in sublats(rsel)  # Each is a Pair s2 => s1
        dns = dniter(rsel)
        for dn in dns
            keepgoing = false
            ijv = builder[dn]
            for j in source_candidates(rsel, s2)
                sitej = allpos[j]
                rsource = sitej - bravais(lat) * dn
                is = targets(builder, rsel.selector.range, rsource, s1)
                for i in is
                    # Make sure we don't stop searching until we reach minimum range
                    is_below_min_range((i, j), (dn, zero(dn)), rsel) && (keepgoing = true)
                    ((i, j), (dn, zero(dn))) in rsel || continue
                    keepgoing = true
                    rtarget = allsitepositions(lat)[i]
                    r, dr = _rdr(rsource, rtarget)
                    v = to_blocktype(term(r, dr), eltype(builder), builder.orbs[s1], builder.orbs[s2])
                    push!(ijv, (i, j, v))
                end
            end
            keepgoing && acceptcell!(dns, dn)
        end
    end
    return nothing
end

# For use in Hamiltonian building
to_blocktype(t::Number, ::Type{T}, t1::NTuple{1}, t2::NTuple{1}) where {T<:Number} = T(t)
to_blocktype(t::Number, ::Type{S}, t1::NTuple{1}, t2::NTuple{1}) where {S<:SMatrix} =
    padtotype(t, S)
to_blocktype(t::SMatrix{N1,N2}, ::Type{S}, t1::NTuple{N1}, t2::NTuple{N2}) where {N1,N2,S<:SMatrix} =
    padtotype(t, S)

to_blocktype(u::UniformScaling, ::Type{T}, t1::NTuple{1}, t2::NTuple{1}) where {T<:Number} = T(u.λ)
to_blocktype(u::UniformScaling, ::Type{S}, t1::NTuple{N1}, t2::NTuple{N2}) where {N1,N2,S<:SMatrix} =
    padtotype(SMatrix{N1,N2}(u), S)

# Fallback to catch mismatched or undesired block types
to_blocktype(t::Array, x...) = throw(ArgumentError("Array input in model, please use StaticArrays instead (e.g. SA[1 0; 0 1] instead of [1 0; 0 1])"))
to_blocktype(t, x...) = throw(DimensionMismatch("Dimension mismatch between model and Hamiltonian. Does the `orbitals` kwarg in your `hamiltonian` match your model?"))

# Although range can be (rmin, rmax) we return all targets within rmax.
# Those below rmin get filtered later by `in rsel`
function targets(builder, range, rsource, s1)
    rmax = maximum(range)
    !isfinite(rmax) && return targets(builder, missing, rsource, s1)
    if !isassigned(builder.kdtrees, s1)
        sitepos = sitepositions(builder.lat.unitcell, s1)
        (builder.kdtrees[s1] = KDTree(sitepos))
    end
    targetlist = inrange(builder.kdtrees[s1], rsource, rmax)
    targetlist .+= builder.lat.unitcell.offsets[s1]
    return targetlist
end

targets(builder, range::Missing, rsource, s1) = siterange(builder.lat, s1)

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

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

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

function unitcell(ham::Hamiltonian{LA,L}; modifiers = (), mincoordination = missing) where {E,L,T,L´,LA<:Superlattice{E,L,T,L´}}
    slat = ham.lattice
    sc = slat.supercell
    supercell_dn = r_to_dn(slat, sc.matrix, SVector{L´}(1:L´))
    pos = allsitepositions(slat)
    br = bravais(slat)
    modifiers´ = resolve.(ensuretuple(modifiers), Ref(slat))

    # Build a version of slat that has a filtered supercell mask according to mincoordination
    slat´ = filtered_superlat!(ham, mincoordination, supercell_dn, br, sc.matrix, pos)
    # store supersite indices newi
    mapping = OffsetArray{Int}(undef, sc.sites, sc.cells.indices...)
    mapping .= 0
    foreach_supersite((s, oldi, olddn, newi) -> mapping[oldi, Tuple(olddn)...] = newi, slat´)

    dim = nsites(slat´.supercell)
    B = blocktype(ham)
    S = typeof(SparseMatrixBuilder{B}(dim, dim))
    harmonic_builders = HamiltonianHarmonic{L´,B,S}[]

    foreach_supersite(slat´) 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
            for p in nzrange(oldh.h, source_i)
                target_i = rows[p]
                wrapped_dn, super_dn = wrap_super_dn(target_i, target_dn, supercell_dn, br, sc.matrix, pos)
                # check: wrapped_dn could exit bounding box along non-periodic direction
                checkbounds(Bool, mapping, target_i, Tuple(wrapped_dn)...) || continue
                newh = get_or_push!(harmonic_builders, super_dn, dim, newcol)
                newrow = mapping[target_i, Tuple(wrapped_dn)...]
                if !iszero(newrow)
                    val = applymodifiers(vals[p], slat, (source_i, target_i), (source_dn, target_dn), modifiers´...)
                    pushtocolumn!(newh.h, newrow, val)
                end
            end
        end
        foreach(h -> finalizecolumn!(h.h), harmonic_builders)
    end
    harmonics = [HamiltonianHarmonic(h.dn, sparse(h.h)) for h in harmonic_builders]
    unitlat = unitcell(slat´)
    return Hamiltonian(unitlat, harmonics, orbitals(ham))
end

filtered_superlat!(sham, ::Missing, args...) = sham.lattice
filtered_superlat!(sham, mc::Int, args...) =
    _filtered_superlat!(expand_supercell_mask(sham), mc, args...)

function _filtered_superlat!(sham::Hamiltonian{LA,L}, mincoord, args...) where {LA,L}
    slat = sham.lattice
    sc = slat.supercell
    mask = sc.mask
    if mincoord > 0
        delsites = NTuple{L+1,Int}[]
        while true
            foreach_supersite(sham.lattice) do _, source_i, source_dn, _
                nn = num_neighbors_supercell(sham.harmonics, source_i, source_dn, mask, args...)
                nn < mincoord && push!(delsites, (source_i, Tuple(source_dn)...))
            end
            foreach(p -> mask[p...] = false, delsites)
            isempty(delsites) && break
            resize!(delsites, 0)
        end
    end
    sc = Supercell(sc.matrix, sc.sites, sc.cells, mask)
    return Superlattice(slat.bravais, slat.unitcell, sc)
end

function num_neighbors_supercell(hhs, source_i, source_dn, mask, args...)
    nn = 0
    for hh in hhs
        ptrs = nzrange(hh.h, source_i)
        rows = rowvals(hh.h)
        nzel = nonzeros(hh.h)
        target_dn = source_dn + hh.dn
        for p in nzrange(hh.h, source_i)
            target_i = rows[p]
            wrapped_dn, _ = wrap_super_dn(target_i, target_dn, args...)
            isonsite = rows[p] == source_i && iszero(hh.dn)
            isincell = isinmask(mask, rows[p], Tuple(wrapped_dn)...)
            isnonzero = !iszero(nzel[p])
            nn += !isonsite && isincell && isnonzero
        end
    end
    return nn
end

function wrap_super_dn(target_i, target_dn, supercell_dn, br, smat, pos)
    r = pos[target_i] + br * target_dn
    super_dn = supercell_dn(r)
    wrapped_dn = target_dn - smat * super_dn
    return wrapped_dn, super_dn
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

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

function applymodifiers(val, lat, inds, dns, m::UniformModifier, ms...)
    selected = (inds, dns) in m.selector
    val´ = selected ? m.f(val) : val
    return applymodifiers(val´, lat, inds, dns, ms...)
end

function applymodifiers(val, lat, (row, col), (dnrow, dncol), m::OnsiteModifier, ms...)
    selected = ((row, col), (dnrow, dncol)) in m.selector
    if selected
        r = allsitepositions(lat)[col] + bravais(lat) * dncol
        val´ = selected ? m(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::HoppingModifier, ms...)
    selected = ((row, col), (dnrow, dncol)) in m.selector
    if selected
        br = bravais(lat)
        r, dr = _rdr(allsitepositions(lat)[col] + br * dncol, allsitepositions(lat)[row] + br * dnrow)
        val´ = selected ? m(val, r, dr) : val
    else
        val´ = val
    end
    return applymodifiers(val´, lat, (row, col), (dnrow, dncol), ms...)
end

#######################################################################
# wrap
#######################################################################
"""
    wrap(h::Hamiltonian, axes; phases = missing)

Build a new Hamiltonian from `h` reducing its dimensions from `L` to `L - length(axes)` by
wrapping the specified Bravais `axes` into a loop. `axes` can be an integer ∈ 1:L or a tuple
of such integers. If `phases` are given (with `length(axes) == length(phases)`), the wrapped
hoppings at a cell distance `dn` along `axes` will be multiplied by a factor
`cis(-dot(phases, dn))`. This is useful, for example, to represent a flux Φ through a loop,
using a single `axes = 1` and `phases = 2π * Φ/Φ₀`.

    wrap(h::Hamiltonian; kw...)

Wrap all axes of `h`, yielding a compactified zero-dimensional Hamiltonian.

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

Curried form equivalent to `wrap(h, axes; kw...)`.

# Examples

```jldoctest
julia> LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = 1/√3)) |>
       unitcell((1,-1), (10, 10)) |> wrap(2)
Hamiltonian{<:Lattice} : 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
```
"""
wrap(h::Hamiltonian, axis::Int; kw...) = wrap(h, (axis,); kw...)

wrap(h::Hamiltonian{<:Lattice,L}; kw...) where {L} = wrap(h, ntuple(identity, Val(L)); kw...)

function wrap(h::Hamiltonian{<:Lattice,L}, axes::NTuple{N,Int}; phases = missing) where {L,N}
    all(axis -> 1 <= axis <= L, axes) && allunique(axes) || throw(ArgumentError("wrap axes should be unique and between 1 and the lattice dimension $L"))
    lattice´ = _wrap(h.lattice, axes)
    phases´ = (phases === missing) ? filltuple(0, Val(N)) : phases
    harmonics´ = _wrap(h.harmonics, axes, phases´, size(h))
    return Hamiltonian(lattice´, harmonics´, orbitals(h))
end

wrap(axes::Union{Integer,Tuple}; kw...) = h -> wrap(h, axes; kw...)

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

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

function _wrap(harmonics::Vector{HamiltonianHarmonic{L,M,A}}, axes::NTuple{N,Int}, phases::NTuple{N,Number}, sizeh) where {L,M,A,N}
    harmonics´ = HamiltonianHarmonic{L-N,M,A}[]
    for har in harmonics
        dn = har.dn
        dn´ = deletemultiple_nocheck(dn, axes)
        phase = -sum(phases .* dn[SVector(axes)])
        newh = get_or_push!(harmonics´, dn´, sizeh)
        # map!(+, newh, newh, factor * har.h) # TODO: activate after resolving #37375
        newh.h .+= cis(phase) .* har.h
    end
    return harmonics´
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(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
#######################################################################

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

Build the Bloch Hamiltonian matrix of `h`, for Bloch phases `ϕs = (ϕ₁, ϕ₂,...)` (or an
`SVector(ϕs...)`). In terms of Bloch wavevector `k`, `ϕs = k * bravais(h)`, it is 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, axis::Int)

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(h::Hamiltonian{<:Lattice}, ϕs::NTuple{L,Real}, dnfunc::Function)

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

    bloch(ph::ParametricHamiltonian, [pϕs, [axis]])

Build the Bloch matrix for `ph`. `pϕs = (ϕs, (;kw...))` or `pϕs = (ϕs..., (;kw...))`
specifies both Bloch phases `ϕs` and the parameters `kw` passed to `ph(; kw...)`. If there are
no `ϕs`, the syntax `pϕs = (;kw...)` is also allowed, which is in that case equivalent to
`bloch(ph(; kw...))`. Similarly, `bloch(ph)` is equivalent to `bloch(ph())`.

    h |> bloch(ϕs, ...)
    ph |> bloch(pϕs, ...)

Curried forms of `bloch`, equivalent to `bloch(h, ϕs, ...)` and `bloch(ph, pϕs, ...)`

# Notes

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

`bloch` also supports slices of `Hamiltonians` and `ParametricHamiltonians`, produced with
`h[rows, cols]`, see the indexing section of `hamiltonian` and  `parametric` for details.

# Examples

```jldoctest
julia> h = LatticePresets.honeycomb() |> hamiltonian(onsite(1) + hopping(2)) |> bloch((0, 0))
2×2 SparseMatrixCSC{Complex{Float64},Int64} with 4 stored entries:
  [1, 1]  =  1.0+0.0im
  [2, 1]  =  6.0+0.0im
  [1, 2]  =  6.0+0.0im
  [2, 2]  =  1.0+0.0im
```

# See also
    `bloch!`, `similarmatrix`
"""
bloch(ϕs, axis = 0) = h -> bloch(h, ϕs, axis)
bloch(h::Hamiltonian, args...) = bloch!(similarmatrix(h), h, args...)

"""
    bloch!(matrix, h::Hamiltonian, [ϕs, [axis]])

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,matrix_type)`, 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 with `Number` eltype for a sparse `h`
with `SMatrix` block eltype).

    bloch!(matrix, ph::ParametricHamiltonian, [pϕs, [axis]])

Same as above but with `pϕs = (ϕs, (;kw...))`, `pϕs = (ϕs..., (;kw...))` or `pϕs = (;kw...)`
(see `bloch` for details).

# Examples

```jldoctest
julia> h = LatticePresets.honeycomb() |> hamiltonian(hopping(2I), orbitals = (Val(2), Val(1)));

julia> bloch!(similarmatrix(h), h, (0, 0))
2×2 SparseMatrixCSC{SMatrix{2, 2, ComplexF64, 4}, Int64} with 2 stored entries:
                     ⋅                       [6.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im]
 [6.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im]                      ⋅

julia> bloch!(similarmatrix(h, flatten), h, (0, 0))
3×3 SparseMatrixCSC{ComplexF64, Int64} with 9 stored entries:
 0.0+0.0im  0.0+0.0im  6.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im
 6.0+0.0im  0.0+0.0im  0.0+0.0im

julia> ph = parametric(h, @hopping!((t; α, β = 0) -> α * t .+ β));

julia> bloch!(similarmatrix(ph, flatten), ph, (0, 0, (; α = 2)))
3×3 SparseMatrixCSC{ComplexF64, Int64} with 9 stored entries:
  0.0+0.0im  0.0+0.0im  12.0+0.0im
  0.0+0.0im  0.0+0.0im   0.0+0.0im
 12.0+0.0im  0.0+0.0im   0.0+0.0im
```

# See also
    `bloch`, `similarmatrix`
"""
bloch!(matrix, h::Hamiltonian, ϕs, axis = 0) = _bloch!(matrix, h, toSVector(ϕs), axis)
bloch!(matrix, h::Hamiltonian, ϕs::Tuple{SVector,NamedTuple}, args...) = bloch!(matrix, h, first(ϕs), args...)

function bloch!(matrix, h::Hamiltonian)
    maybe_flatten_copy!(parent(matrix), first(h.harmonics).h, orbitalstructure(h)) # faster copy!(dense, sparse) specialization
    return matrix
end

function _bloch!(matrix::AbstractMatrix, h::Hamiltonian{<:Lattice,L,M}, ϕs::SVector{L}, axis::Number) where {L,M}
    rawmatrix = parent(matrix)
    if iszero(axis)
        maybe_flatten_copy!(rawmatrix, first(h.harmonics).h, orbitalstructure(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::SVector{L}, dnfunc::Function) where {L,M}
    prefactor0 = dnfunc(zero(ϕs))
    rawmatrix = parent(matrix)
    if iszero(prefactor0)
        fill!(rawmatrix, zero(eltype(rawmatrix)))
    else
        maybe_flatten_copy!(rawmatrix, first(h.harmonics).h, orbitalstructure(h))
        rmul!(rawmatrix, prefactor0)
    end
    add_harmonics!(rawmatrix, h, ϕs, dnfunc)
    return matrix
end

_bloch!(matrix, h::Hamiltonian{<:Lattice,L}, ϕs::SVector{L´}, axis) where {L,L´} =
    L == L´ ? throw(ArgumentError("Unexpected `bloch!` signature")) :
              throw(DimensionMismatch("Mismatch between $L-dimensional Hamiltonian and $L´-dimensional Bloch phases"))

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)
        maybe_flatten_add!(zerobloch, hhmatrix, orbitalstructure(h), ephi)
    end
    return zerobloch
end

############################################################################################
######## maybe_flatten_copy! and maybe_flatten_add! call specialized methods #########################################
############################################################################################

maybe_flatten_copy!(dest, src, _) = copy!(dest, src)
maybe_flatten_copy!(dst::AbstractMatrix{<:Number}, src::SparseMatrixCSC{<:Number}, o) = fast_sparse_copy!(dst, src)
maybe_flatten_copy!(dst::StridedMatrix{<:Number}, src::SparseMatrixCSC{<:Number}, o) = fast_sparse_copy!(dst, src)
maybe_flatten_copy!(dst::StridedMatrix{<:SMatrix{N,N}}, src::SparseMatrixCSC{<:SMatrix{N,N}}, o) where {N} = fast_sparse_copy!(dst, src)
maybe_flatten_copy!(dst::AbstractMatrix{<:Number}, src::SparseMatrixCSC{<:SMatrix}, o) = flatten_sparse_copy!(dst, src, o)
maybe_flatten_copy!(dst::StridedMatrix{<:Number}, src::StridedMatrix{<:SMatrix}, o) = flatten_dense_copy!(dst, src, o)

maybe_flatten_add!(dest, src, _, α) = _plain_muladd!(dest, src, α)
maybe_flatten_add!(dst::AbstractMatrix{<:Number}, src::SparseMatrixCSC{<:Number}, o, α = 1) = fast_sparse_muladd!(dst, src, α)
maybe_flatten_add!(dst::AbstractMatrix{<:SMatrix{N,N}}, src::SparseMatrixCSC{<:SMatrix{N,N}}, o, α = I) where {N} = fast_sparse_muladd!(dst, src, α)
maybe_flatten_add!(dst::AbstractMatrix{<:Number}, src::SparseMatrixCSC{<:SMatrix}, o, α = I) = flatten_sparse_muladd!(dst, src, o, α)
maybe_flatten_add!(dst::StridedMatrix{<:Number}, src::StridedMatrix{<:SMatrix}, o, α = I) = flatten_dense_muladd!(dst, src, o, α)

max_stars_repo_name: pablosanjose/Quantica

max_issues_repo_name: pablosanjose/Quantica

max_forks_repo_name: pablosanjose/Quantica


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

mutable struct TerminateNotificationProfile <: SwaggerModel
    notBeforeTimeout::Any # spec type: Union{ Nothing, String } # spec name: notBeforeTimeout
    enable::Any # spec type: Union{ Nothing, Bool } # spec name: enable

    function TerminateNotificationProfile(;notBeforeTimeout=nothing, enable=nothing)
        o = new()
        validate_property(TerminateNotificationProfile, Symbol("notBeforeTimeout"), notBeforeTimeout)
        setfield!(o, Symbol("notBeforeTimeout"), notBeforeTimeout)
        validate_property(TerminateNotificationProfile, Symbol("enable"), enable)
        setfield!(o, Symbol("enable"), enable)
        o
    end
end # type TerminateNotificationProfile

const _property_map_TerminateNotificationProfile = Dict{Symbol,Symbol}(Symbol("notBeforeTimeout")=>Symbol("notBeforeTimeout"), Symbol("enable")=>Symbol("enable"))
const _property_types_TerminateNotificationProfile = Dict{Symbol,String}(Symbol("notBeforeTimeout")=>"String", Symbol("enable")=>"Bool")
Base.propertynames(::Type{ TerminateNotificationProfile }) = collect(keys(_property_map_TerminateNotificationProfile))
Swagger.property_type(::Type{ TerminateNotificationProfile }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_TerminateNotificationProfile[name]))}
Swagger.field_name(::Type{ TerminateNotificationProfile }, property_name::Symbol) =  _property_map_TerminateNotificationProfile[property_name]

function check_required(o::TerminateNotificationProfile)
    true
end

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

ext: jl


text:
using HorizonSideRobots
include("../basics.jl")

md"""
            ДАНО: Где-то на неограниченном со всех сторон поле и без внутренних перегородок имеется единственный маркер.
                            Робот - в произвольной клетке поля.

            РЕЗУЛЬТАТ: Робот - в клетке с тем маркером.
"""

"""Двигает робота на {x} шагов или останавливается на маркере"""
function move_for_x_or_find_marker!(r, side, x)
    for _ in 1:x
        move!(r, side)
        if ismarker(r)
            return false
        end
    end
    return true
end

"""Двигает робота по спирали пока он не найдет маркер"""
function radial_search!(r)
    x = 1
    count = 0
    side = Nord

    while move_for_x_or_find_marker!(r, side, x)
        count += 1
        side = next_side(side)
        if count % 2 == 0
            x += 1
        end
    end
end

 """Решение задачи 9"""
 function solve_9!(r)
    radial_search!(r)

    show(r)
end


text:
using .DifferentialEquations

include("../models/diff_eq_monte_carlo.jl")

function predict_output_type_zero(x::absdiffeqmodel) where {absdiffeqmodel <: MonteCarloDiffEqModel}
    return zero(Float64)
end


text:
# plink related commands
using ABG

#allele_maf_stats  = ABG.allele_maf_stats
#bed_2_map_n_ped   = ABG.bed_2_map_n_ped
#bed_snp_subset    = ABG.bed_snp_subset
#fr2ped            = ABG.fr2ped
#hwe_stats         = ABG.hwe_stats
#merge_beds        = ABG.merge_beds
#miss_allele_stats = ABG.miss_allele_stats
#ped_n_map_to_bed  = ABG.ped_n_map_to_bed
#plink_2_vcf       = ABG.plink_2_vcf
#plink_filter_id   = ABG.plink_filter_id
#plink_filter_snp  = ABG.plink_filter_snp
#vcf_2_plink       = ABG.vcf_2_plink
#plink_keep_id     = ABG.plink_keep_id
#plink_012         = ABG.plink_012
#force_reference_allele = ABG.force_reference_allele

# styled printing
title    = ABG.title
subtitle = ABG.subtitle
item     = ABG.item
message  = ABG.message
warning  = ABG.warning
done     = ABG.done


text:
"""
This module is shared between the server and notebook processes.

It contains several data types definitions which should be understood on both
end by the serialization machinery.
"""
module PlutoShared

"""
Data structure which reprensets a single element of a User Interface.

Accidentally this was made to match the representation of React elements (the
data structure produced by React.createElement(...) or JSX syntax).
"""
struct UI
  type::String
  props::Dict{Symbol,Any}
end

UI(type::String; kwargs...) =
  UI(type, Dict(kwargs...))

end


text:
# General functions for an estimator
# Does not depend on dimentionality of support

export NeuralDensityEstimator, component_weights, condition!

"""
    autostatic_array(a)

function to automatically convert `a` into an appropriatly sized StaticArray
"""
autostatic_array(a) = _autostatic_array(Size(size(a)), a);
_autostatic_array(S, a) = S(a)

immutable NeuralDensityEstimator{N}
    sess::Session

    #Network nodes
    optimizer::Tensor
    conditioner::Tensor
    t::Tensor
    pdf::Tensor
end

for (outer, inner) in [(:(Distributions.pdf), :_pdf), (:(Distributions.loglikelihood), :_loglikelihood)]
    @eval begin
        # Entry
        function $outer(est::NeuralDensityEstimator, ts::AbstractArray)
        #    @show (1, typeof(ts))
            $outer(est, autostatic_array(ts))
        end

        # Single observation
        function $outer(est::NeuralDensityEstimator{1}, t::Real)
        #    @show (2, typeof(ts))
            $inner(est, SMatrix{1,1}(t)) |> first
        end

        function $outer(est::NeuralDensityEstimator{1}, ts::StaticVector{1})
        #    @show (3, typeof(ts))
            $inner(est, reshape(ts[1], (1,1)))
        end

        function $outer{N}(est::NeuralDensityEstimator{N}, ts::StaticVector{N})
        #    @show (4, typeof(ts))
            $inner(est, reshape(ts,Val(2)))
        end

        #M is number of observations
        function $outer{M}(est::NeuralDensityEstimator{1}, ts::StaticVector{M})
        #    @show (5, typeof(ts))
            $inner(est, reshape(ts, (1,M)))
        end
    end
end

function _pdf{N}(est::NeuralDensityEstimator{N}, ts::AbstractMatrix)
    size(ts,1) == N || error("wrong dimensions, expected $N, was given $(size(ts,1))")
    gr = est.sess.graph
    run(est.sess, est.pdf, Dict(est.t=>ts)) |> vec
end

function _loglikelihood{N}(est::NeuralDensityEstimator{N}, ts::AbstractMatrix)
    size(ts,1) == N || error("wrong dimensions, expected $N, was given $(size(ts,1))")
    gr = est.sess.graph
    run(est.sess, gr["loglikelihood"], Dict(est.t=>ts))
end

"""
    component_weights(est::NeuralDensityEstimator)

Extracts the final layer weights from the estimator.
These are similar to the mixture weights in a GMM
"""
function component_weights(est::NeuralDensityEstimator)
    node_names =  collect(keys(est.sess.graph))
    final_layer = node_names[ismatch.(r"^W_\d+_squared$",node_names)] |> sort |> last
    vec(run(est.sess, est.sess.graph[final_layer]))
end

ignore(args...) = nothing

function StatsBase.fit!(estimator::NeuralDensityEstimator{1}, observations::AbstractVector; kwargs...)
    fit!(estimator, reshape(observations, (1,length(observations))) ; kwargs...)
end

function StatsBase.fit!{N}(estimator::NeuralDensityEstimator{N}, observations::AbstractMatrix;
    callback=ignore, callback_vars=["ysmin", "ysmax", "loglikelihood", "working_loss"],
    epochs = 20_000)

    size(observations,1)==N || throw(ArgumentError("Estimator is for $(N) dimensional domain, but observations are $(size(observations,1)) dimensional"))

    gr = estimator.sess.graph
    callback_nodes = [gr[var] for var in callback_vars]

    for ii in 1:epochs
        outs = run(estimator.sess,
            [callback_nodes..., estimator.optimizer],
            Dict(estimator.t=>observations))

        if ii % 100 == 1
            callback(ii, Dict(zip(callback_vars, outs[1:end-1])))
        end

    end
    estimator
end

"""
    condition(est::NeuralDensityEstimator tol = 1e-15, max_epochs=2_000)

"Conditions" the neural density estimate so the support extrema are mapped to 1. and 2.
This improves training by adjusting the area the network has the learn over

"""
function condition!(est::NeuralDensityEstimator; tol = 1e-15, max_epochs=2_000)
    gr = est.sess.graph
    for ii in 1:max_epochs
        _, ysmin, ysmax, condition_loss = run(est.sess, [est.conditioner, gr["ysmin"],gr["ysmax"], gr["condition_loss"]])
        ii % 50 == 1 && @show (ii, ysmin, ysmax, condition_loss)
        if condition_loss[1] < 1e-15
            break
        end
    end
end


text:
import Base.Sort, Base.Sort.Algorithm
import Base:sort!

struct NOOPSortAlg <: Base.Sort.Algorithm end
const NOOPSort = NOOPSortAlg()

sort!(x, ::Integer, ::Integer, ::NOOPSortAlg, ::Base.Sort.Ordering) = x

struct BFS{T<:Base.Sort.Algorithm} <: TraversalAlgorithm
    sort_alg::T
end

BFS() = BFS(NOOPSort)

"""
    traverse_graph!(g, s, alg, state, neighborfn=outneighbors)
    traverse_graph!(g, ss, alg, state, neighborfn=outneighbors)

    Traverse a graph `g` starting at vertex `s` / vertices `ss` using algorithm `alg`, maintaining state in [`AbstractTraversalState`](@ref) `state`. Next vertices to be visited are determined by `neighborfn` (default `outneighbors`). Return `true` if traversal finished; `false` if one of the visit functions caused an early termination.

"""
function traverse_graph!(
    g::AbstractGraph{U},
    ss::AbstractVector,
    alg::BFS,
    state::AbstractTraversalState,
    neighborfn::Function=outneighbors
    ) where U<:Integer

    n = nv(g)
    visited = falses(n)
    cur_level = Vector{U}()
    sizehint!(cur_level, n)
    next_level = Vector{U}()
    sizehint!(next_level, n)
    @inbounds for s in ss
        us = U(s)
        visited[us] = true
        push!(cur_level, us)
        initfn!(state, us) || return false
    end
    while !isempty(cur_level)
        @inbounds for v in cur_level
            previsitfn!(state, v) || return false
            @inbounds @simd for i in neighborfn(g, v)
                visitfn!(state, v, i) || return false
                if !visited[i]
                    newvisitfn!(state, v, i) || return false
                    push!(next_level, i)
                    visited[i] = true
                end
            end
            postvisitfn!(state, v) || return false
        end
        postlevelfn!(state) || return false
        empty!(cur_level)
        cur_level, next_level = next_level, cur_level
        sort!(cur_level, alg=alg.sort_alg)
    end
    return true
end

mutable struct DistanceState{T<:Integer} <: AbstractTraversalState
    distances::Vector{T}
    n_level::T
end

@inline function initfn!(s::DistanceState{T}, u) where T
    s.distances[u] = zero(T)
    return true
end
@inline function newvisitfn!(s::DistanceState, u, v)
    s.distances[v] = s.n_level
    return true
end
@inline function postlevelfn!(s::DistanceState{T}) where T
    s.n_level += one(T)
    return true
end

"""
    distances(g, s, alg=BFS())
    distances(g, ss, alg=BFS())

Return a vector filled with the geodesic distances of vertices in  `g` from vertex `s` / unique vertices `ss`
using BFS traversal algorithm `alg`.  For vertices in disconnected components the default distance is `typemax(T)`.
"""
function distances(g::AbstractGraph{T}, s, alg::BFS=BFS()) where T
    d = fill(typemax(T), nv(g))
    state = DistanceState(d, one(T))
    traverse_graph!(g, s, alg, state)
    return state.distances
end


text:
using Statistics
using .DCHSBM: sumOfExteriorDegreesIntensityFunction, sampleEdges
using .Kronecker: hyperkron_graph, kron_params

## DCHSBM

function old_dchsbm(size::Real, density=1/40; kmax::Int64=3)
    density_modifier = density*40
    n = round(Integer, 4*(size/density_modifier)^(1/kmax)) # modified to change n from 20
    #begin copied from https://github.com/nveldt/HyperModularity.jl/blob/0478acb5e209550cbaf7572f4151c4569f64faf6/test/runtests.jl#L20
    Z = rand(1:5, n)
    Z = Int64.(Z)
    ϑ = dropdims(ones(1,n) + rand(1,n), dims = 1)
    μ = mean(ϑ)

    ω(x, α) = density_modifier*(x[1]+1)^(-α[x[2]]*x[2]) # modified to add density_modifier

    α = repeat([2.0], kmax)

    Ω = sumOfExteriorDegreesIntensityFunction(ω, kmax)
    sampleEdges(Z, ϑ, Ω;α=α, kmax=kmax, kmin = 1) # modified to skip computing degrees
    #end copied from https://github.com/nveldt/HyperModularity.jl/blob/0478acb5e209550cbaf7572f4151c4569f64faf6/test/runtests.jl#L20
end
function new_dchsbm(size::Real, density=1/40; kmax::Int64=3)
    #modeled after old_dchsbm
    density_modifier = density*40
    n = round(Integer, 4*(size/density_modifier)^(1/kmax))

    Z = rand(1:5, n)
    ϑ = rand(n) .+ 1

    Ω = FBD.inverse_power_intensity_function(8)

    sampler = DCHSBM_sampler(Ω, Z, ϑ, kmax)

    rand(sampler, round(Integer, size/kmax))
end

## Kronecker

# parameters copied from https://www.cs.purdue.edu/homes/dgleich/codes/hyperkron/
const KRON_PARAMS = kron_params(0.99, 0.2, 0.3, 0.05)
# With a little math, density can be turned into an optional argument with a default value
# dependent on size
function old_kronecker(size)
    power = max(1, round(Integer, log(2*size)/log(sum(KRON_PARAMS))))

    A,hedges = hyperkron_graph(KRON_PARAMS, power)
    hedges
end
function old_kronecker(size, density)
    power = max(1, round(Integer, log(2*size/density)/log(8)))
    density_modifier = density^(1/power)/mean(KRON_PARAMS)

    A,hedges = hyperkron_graph(min.(.99, density_modifier*KRON_PARAMS), power)
    hedges
end
function new_kronecker(size)
    power = max(1, round(Integer, log(2*size)/log(sum(KRON_PARAMS))))

    sampler = Kronecker_sampler(KRON_PARAMS, power, space=min(size÷10,30000))

    rand(sampler, round(Integer, size/3))
end
function new_kronecker(size, density)
    power = max(1, round(Integer, log(2*size/density)/log(8)))

    sampler = Kronecker_sampler(KRON_PARAMS, power, space=min(size÷10,30000))

    rand(sampler, round(Integer, size/3))
end

## HyperPA

function old_hyperpa(size::Integer=45000; name = "DAWN")
    output_directory = joinpath(@__DIR__, "..", "data", "hyper_pa")
    output_file = joinpath(output_directory, "$name.txt")
    mkpath(output_directory)
    rm(output_file, force=true)

    nodes::Integer = size ÷ 222 + 18
    nodes ≥ 200 && @warn "expected runtime on the order of $(round(Integer, (nodes/200)^2)) minute(s)\n"

    cd(HyperPA.source()) do
        run(`python3 hyper_preferential_attachment.py --name=$name --file_name=$name --num_nodes=$nodes --simplex_per_node_directory='simplex per node' --size_distribution_directory='size distribution' --output_directory=$output_directory`)
        isfile(output_file) ? [parse.(Int, split(line)) for line in readlines(output_file)] : Vector{Int}[]
    end
end
function new_hyperpa(size::Integer=45000)
    example(hyper_pa, size)
end

test_functions = [(old_dchsbm, new_dchsbm), (old_kronecker, new_kronecker), (old_hyperpa, new_hyperpa)]


text:
struct Boat_timing
    T_in::Float64
    T_out::Float64
    T_1::Float64
    T_2::Float64
    T_t::Float64
end

# T_1 = 2 * (0.7 + T_in + T_out)
Boat_timing(T_active, T_passive) =
    Boat_timing(0, 0, 2T_active, 2T_passive, T_active + T_passive)

Boat_timing() = Boat_timing(0.7, 0.5)

struct Boat_angle_ranges
    γ_OA::Tuple{Float64,Float64}
    θ_t::Tuple{Float64,Float64}
    θ_k::Tuple{Float64,Float64}
end

Boat_angle_ranges() = Boat_angle_ranges((30.0, 120.0), (22.5, 22.5), (90.0, 90.0))

function Base.deg2rad(angles::Boat_angle_ranges)
    Boat_angle_ranges(deg2rad.(angles.γ_OA), deg2rad.(angles.θ_t), deg2rad.(angles.θ_k))
end

struct Boat_angles
    γ_OA::Float64
    θ_t::Float64
    θ_k::Float64
    γ′_OA::Float64
    θ′_t::Float64
    θ′_k::Float64
    γ′′_OA::Float64
    θ′′_t::Float64
    θ′′_k::Float64
end

function Boat_angles(angle_ranges::Boat_angle_ranges, bt::Boat_timing, t)
    aranges = angle_ranges.γ_OA, angle_ranges.θ_t, angle_ranges.θ_k
    γoa, θt, θk = map(aranges) do r
        _angle(t, bt, r[1], r[2])
    end
    γ′oa, θ′t, θ′k = map(aranges) do r
        _angle′(t, bt, r[1], r[2])
    end
    γ′′oa, θ′′t, θ′′k = map(aranges) do r
        _angle′′(t, bt, r[1], r[2])
    end
    Boat_angles(γoa, θt, θk, γ′oa, θ′t, θ′k, γ′′oa, θ′′t, θ′′k)
end

Boat_angles() = Boat_angles(deg2rad(Boat_angle_ranges()), Boat_timing(), 0.0)

t_periodic(t, bt::Boat_timing) = mod(t, bt.T_t)

function η(t, bt::Boat_timing)
    # to - czas w danym okresie
    to = t_periodic(t, bt)

    if to < bt.T_in
        to / bt.T_in
    elseif to < 0.5bt.T_1 - bt.T_out
        one(t)
    elseif to < 0.5bt.T_1
        1 - (to - 0.5bt.T_1 + bt.T_out) / bt.T_out
    else
        zero(t)
    end
end

function kpars(to, bt::Boat_timing)
    #KPARS Helper to get phase parameters.
    k1 = to <= 0.5bt.T_1
    (k1, !k1) .* 1
end

function _angle(t, bt::Boat_timing, vmin, vmax)
    to = t_periodic(t, bt)
    k1, k2 = kpars(to, bt)
    T_1, T_2 = bt.T_1, bt.T_2
    0.5(vmax + vmin) -
    0.5(vmax - vmin) * (k1 * cospi(2to / T_1) - k2 * cospi(2(to - 0.5T_1) / T_2))
end

function _angle′(t, bt::Boat_timing, vmin, vmax)
    to = t_periodic(t, bt)
    k1, k2 = kpars(to, bt)
    T_1, T_2 = bt.T_1, bt.T_2
    -π *
    ((k2 * sinpi(2(T_1 / 2 - to) / T_2)) / T_2 + (k1 * sinpi(2to / T_1)) / T_1) *
    (vmin - vmax)
end

function _angle′′(t, bt::Boat_timing, vmin, vmax)
    to = t_periodic(t, bt)
    k1, k2 = kpars(to, bt)
    T_1, T_2 = bt.T_1, bt.T_2
    2π^2 *
    (k2 * cospi(2(T_1 / 2 - to) / T_2) / T_2^2 - k1 * cospi(2to / T_1) / T_1^2) *
    (vmin - vmax)
end


text:
module Octrees

using ClusterTrees
# using StaticArrays

struct Data{T}
    sector::Int
    values::Vector{T}
end

"""
An octree behaves like a pointer based tree, specifically designed for the
storage of data residing in Euclidian space. The routing and insertion methods
maintain a Hilbert space filling curve order.

This makes the Octree ideal as the basis for ordering algorithms aimed at
preserving data locality (in the context of e.g. parallellism).
"""
struct Octree{T}
    pbtree::ClusterTrees.PointerBasedTrees.PointerBasedTree{ClusterTrees.PointerBasedTrees.Node{Data{T}}}
    function Octree(value::Vector{T}) where {T}
        N = ClusterTrees.PointerBasedTrees.Node{Data{T}}
        new{T}(
            ClusterTrees.PointerBasedTrees.PointerBasedTree(
                N[N(Data(0,value),0,0,0,0)], 1))
    end
end

Base.length(itr::ClusterTrees.DepthFirstIterator{T}) where {T<:Octree} = length(itr.tree.pbtree.nodes)

# Tree API
ClusterTrees.data(tree::Octree, node) = data(tree.pbtree, node)
ClusterTrees.children(tree::Octree, node) = children(tree.pbtree, node)
ClusterTrees.root(tree::Octree) = root(tree.pbtree)
# ClusterTrees.insert!(tree::Octree, data; next, parent, prev) = ClusterTrees.insert!(tree.pbtree, data, next=next, parent=parent, prev=prev)

import ClusterTrees.ChildIterator
Base.insert!(chds::ChildIterator{<:Octree}, item, pos) = insert!(tree.pbtree, item, pos)

# PointerBasedTrees API
ClusterTrees.PointerBasedTrees.getnode(tree::Octree, node)  = ClusterTrees.PointerBasedTrees.getnode(tree.pbtree, node)

const hilbert_states = [
    [1, 2, 3, 2, 4, 5, 3, 5],
    [2, 6, 0, 7, 8, 8, 0, 7],
    [0, 9,10, 9, 1, 1,11,11],
    [6, 0, 6,11, 9, 0, 9, 8],
    [11,11, 0, 7, 5, 9, 0, 7],
    [4, 4, 8, 8, 0, 6,10, 6],
    [5, 7, 5, 3, 1, 1,11,11],
    [6, 1, 6,10, 9, 4, 9,10],
    [10, 3, 1, 1,10, 3, 5, 9],
    [4, 4, 8, 8, 2, 7, 2, 3],
    [7, 2,11, 2, 7, 5, 8, 5],
    [10, 3, 2, 6,10, 3, 4, 4]]

const hilbert_positions = [
    [0,1,3,2,7,6,4,5],
    [0,7,1,6,3,4,2,5],
    [0,3,7,4,1,2,6,5],
    [2,3,1,0,5,4,6,7],
    [4,3,5,2,7,0,6,1],
    [6,5,1,2,7,4,0,3],
    [4,7,3,0,5,6,2,1],
    [6,7,5,4,1,0,2,3],
    [2,5,3,4,1,6,0,7],
    [2,1,5,6,3,0,4,7],
    [4,5,7,6,3,2,0,1],
    [6,1,7,0,5,2,4,3]]

"""
Return the center of the child box and its sector based on the relative
positioning of the parent center `ct` and the target point `pt`. Descending
a level will half the *half size* `hs`.
"""
function sector_center_size(pt, ct, hs)
    hs = hs / 2
    bl = pt .> ct
    ct = ifelse.(bl, ct.+hs, ct.-hs)
    sc = sum(b ? 2^(i-1) : 0 for (i,b) in enumerate(bl))
    return sc, ct, hs
end

struct Router{T,P}
    smallest_box_size::T
    target_point::P
end

import ClusterTrees: start, next, done

function route!(tree::Octree, state, destination)

    point = destination.target_point
    smallest_box_size = destination.smallest_box_size

    node_idx, center, size, sfc_state = state
    size <= smallest_box_size && return state
    target_sector, target_center, target_size = sector_center_size(point, center, size)
    target_pos = hilbert_positions[sfc_state][target_sector+1] + 1
    target_sfc_state = hilbert_states[sfc_state][target_sector+1] + 1

    chds = children(tree, node_idx)
    pos = start(chds)
    while !done(chds, pos)
        child, newpos = next(chds, pos)
        child_sector = ClusterTrees.data(tree, child).sector
        child_pos = hilbert_positions[sfc_state][child_sector+1]+1
        target_pos < child_pos && break
        if child_sector == target_sector
            return child, target_center, target_size, target_sfc_state
        end
        pos = newpos
    end

    data = Data(target_sector, Int[])
    child = insert!(chds, data, pos)

    return child, target_center, target_size, target_sfc_state
end

end # module


text:
# # MSA and structures
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__notebooks/11_PDB.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__notebooks/11_PDB.ipynb)
#
# For this example, we use the Multiple Sequence Alignment (MSA) of the
# *Mu DNA-binding domain*.

using MIToS.Pfam
using MIToS.MSA

const pfam_file = downloadpfam("PF02316")
const msa = read(pfam_file,
				 Stockholm,
				 generatemapping=true,
				 useidcoordinates=true)

# We use the `getseq2pdb` function to look into the Pfam annotations for PDBs:

const seq2pdb = getseq2pdb(msa)

# For this example, we use the crystallographic structure 4FCY...

selected = [ (seq, pdb, chain) for (seq, pdbs) in seq2pdb
				for (pdb, chain) in pdbs if pdb == "4FCY" ]

# ...and we take the first one

seq_id, pdb_id, chain = selected[1]

# We download and read the PDB file

using MIToS.PDB
pdb_file = downloadpdb(pdb_id, format=PDBFile)

#-

const pdb_res = read(pdb_file, PDBFile)

# Each PDBResidue contains the information in the ATOM and HETATM PDB lines:

first_residue = pdb_res[1]

# You can access the information by accessing the field names:

first_residue.id.name

# #### Exercise 1
# 4
#
# How many missing residues are in the chain A

# How many chains are in the PDB?
#
# **Hint:** You can use `Set` or `unique`

#src chains = Set([ res.id.chain for res in pdb_res ])
#src length(chains)

# There are 5 chains: A, B, C, D, E

# ## Distances and contacts
#
# MIToS PDB module has functions to measure distance between residues and
# identify contacts:
#

res_i = pdb_res[1]
res_j = pdb_res[4]
distance(res_i, res_j)

# `contact(i, j, threshold)` is faster than `distance(i, j) < threshold`

contact(res_i, res_j, 8.0)

# `distance` and `contact` can take a `criteria` keyword argument with one of
# the following values: `Heavy`, `All`, `CA`, `CB` (it uses `CA` for `GLY`).

distance(res_i, res_j, criteria="CB")

# ### Exercise 2
# 4
#
# How many missing residues are in the chain A

# Write a function that returns the `Set` of residues in the first vector that
# are in contact with the residues in the second vector:

function get_contacts(residues_i, residues_j; threshold::Float64=8.0)
	result = Set{PDBResidue}()
##   ...your code here...
#src	for res_i in residues_i
#src		for res_j in residues_j
#src			if contact(res_i, res_j, threshold)
#src				push!(result, res_i)
#src			end
#src		end
#src	end
	result
end

#-

using Test
@test get_contacts(pdb_res[[1, end]], pdb_res[2:3], threshold=6.0) == Set([pdb_res[1]])

# ### Exercise 3
# 4
#
# How many missing residues are in the chain A

# Use the `get_contacts` function to get all the residues from the chain A
# (protein) that are in contact with the chains C, D and E (DNA).

## ...your solution...

#src get_contacts([res for res in pdb_res if res.id.chain == "A"],
#src 			  [res for res in pdb_res if res.id.chain in Set(["C", "D", "E"])])

# ## PDB Plots
#
# You can use Plots to get an idea of where your residues are:

using Plots
plotlyjs()

#-

plot(pdb_res)

# Or use Bio3DView to get an interactive view (but we need to save the residues
# in an uncompressed PDB file first):

write("selected_residues.pdb", pdb_res, PDBFile)

#-

using Bio3DView
viewfile("selected_residues.pdb", "pdb")

# ## SIFTS
#
# We can download the XML SIFTS file of the PDB using the SIFTS module of MIToS.
# It has a residue level mapping between databases and information about the
# structure.

using MIToS.SIFTS

const sifts_file = downloadsifts(pdb_id)

#-

const sifts_res = read(sifts_file, SIFTSXML)

# Similar to `PDBResidue`, a `SIFTSResidue` contains information about a
# single residue and their fields can be accessed using dots:

sifts_res[1]

# In this case, the first residue is missing in the structure:

sifts_res[1].missing

#-

sifts_res[end]

#-

sifts_res[end].missing

# #### Exercise 4
#
# How many missing residues are in the chain A?

## ...your code here...
# sum(res.missing for res in sifts_res if res.PDB.chain == "A")

max_stars_repo_name: diegozea/JuliaForBioinformatics

max_issues_repo_name: diegozea/JuliaForBioinformatics

max_forks_repo_name: diegozea/JuliaForBioinformatics


text:
@testset "test for voltage deviation" begin
    @testset "case_vstab: cigre mv" begin
        result = run_pandamodels_vstab(case_vstab)
        pm = _PdM.load_pm_from_json(case_vstab)
        params = _PdM.extract_params!(pm)

        @test string(result["termination_status"]) == "LOCALLY_SOLVED"
        @test string(result["dual_status"]) == "FEASIBLE_POINT"
        @test string(result["primal_status"]) == "FEASIBLE_POINT"

        bus = result["solution"]["bus"]
        for (idx, bus) in result["solution"]["bus"]
            if idx in keys(params[:setpoint_v])
                @test isapprox(bus["vm"], params[:setpoint_v][idx]["value"], atol=1e-1)
            end
        end

        @test isapprox(result["objective_lb"], -Inf)
        @test isapprox(result["objective"], 0.000453688; atol = 1e-8)
        @test result["solve_time"] > 0.0

    end

    @testset "case_qflex: cigre mv" begin
        result = run_pandamodels_qflex(case_qflex)
        pm = _PdM.load_pm_from_json(case_qflex)
        params = _PdM.extract_params!(pm)

        @test string(result["termination_status"]) == "LOCALLY_SOLVED"
        @test string(result["dual_status"]) == "FEASIBLE_POINT"
        @test string(result["primal_status"]) == "FEASIBLE_POINT"

        for (idx, br) in result["solution"]["branch"]
            if idx in keys(params[:setpoint_q])
                @test isapprox(br["qf"], params[:setpoint_q][idx]["value"], atol=1e-1)
            end
        end
        @test isapprox(result["objective_lb"], -Inf)
        @test result["solve_time"] > 0.0
    end

    @testset "case_multi_vstab: cigre mv" begin
        result = run_pandamodels_multi_vstab(case_multi_vstab)
        pm = _PdM.load_pm_from_json(case_multi_vstab)
        params = _PdM.extract_params!(pm)

        @test string(result["termination_status"]) == "LOCALLY_SOLVED"
        @test string(result["dual_status"]) == "FEASIBLE_POINT"
        @test string(result["primal_status"]) == "FEASIBLE_POINT"

        # bus = result["solution"]["bus"]
        # for (idx, bus) in result["solution"]["bus"]
        #     if idx in keys(params[:setpoint_v])
        #         @test isapprox(bus["vm"], params[:setpoint_v][idx]["value"], atol=1e-1)
        #     end
        # end

        @test result["solve_time"] > 0.0
    end
end