text:
#######################################################################
# Subspace
#######################################################################
struct Subspace{D,C,T,S<:AbstractMatrix{C},O<:OrbitalStructure}
    energy::T
    basis::S
    orbstruct::O
    basevert::SVector{D,T}
end
Subspace(h::Hamiltonian, args...) = Subspace(orbitalstructure(h), args...)
Subspace(::Missing, energy, basis, basevert...) =
    Subspace(OrbitalStructure(eltype(basis), size(basis, 1)), energy, basis, basevert...)
Subspace(o::OrbitalStructure, energy, basis, basevert...) =
    Subspace(energy, unflatten_orbitals_or_reinterpret(basis, o), o, basevert...)
Subspace(energy::T, basis, orbstruct) where {T} = Subspace(energy, basis, orbstruct, SVector{0,T}())
Subspace(energy::T, basis, orbstruct, basevert) where {T} = Subspace(energy, basis, orbstruct, SVector(T.(basevert)))

function Base.show(io::IO, s::Subspace{D,C,T}) where {D,C,T}
    i = get(io, :indent, "")
    print(io,
"$(i)Subspace{$D}: eigenenergy subspace on a $(D)D manifold
$i  Energy       : $(s.energy)
$i  Degeneracy   : $(degeneracy(s))
$i  Bloch/params : $(s.basevert)
$i  Basis eltype : $C")
end

"""
    degeneracy(s::Subspace)

Return the degeneracy of a given energy subspace. It is equal to `size(s.basis, 2)`.

# See also
    `spectrum`, `bandstructure`
"""
degeneracy(s::Subspace) = degeneracy(s.basis)
degeneracy(m::AbstractMatrix) = isempty(m) ? 1 : size(m, 2)  # To support sentinel empty projs

orbitalstructure(s::Subspace) = s.orbstruct

flatten(s::Subspace) =
    Subspace(s.energy, flatten(s.basis, s.orbstruct), s.orbstruct, s.basevert)

unflatten(s::Subspace, o::OrbitalStructure) = Subspace(s.energy, unflatten_orbitals(s.basis, o), o, s.basevert)

# destructuring
Base.iterate(s::Subspace) = s.energy, Val(:basis)
Base.iterate(s::Subspace, ::Val{:basis}) = s.basis, Val(:done)
Base.iterate(::Subspace, ::Val{:done}) = nothing
Base.IteratorSize(::Subspace) = Base.HasLength()
Base.first(s::Subspace) = s.energy
Base.last(s::Subspace) = s.basis
Base.length(s::Subspace) = 2

Base.copy(s::Subspace) = deepcopy(s)

#######################################################################
# Spectrum
#######################################################################
struct Spectrum{D,C,T,S<:AbstractMatrix{C},E<:AbstractVector{T},M<:Diagonalizer}
    energies::E
    states::S
    diag::M
    basevert::SVector{D,T}
    subs::Vector{UnitRange{Int}}
end

Spectrum(energies::AbstractVector{T}, states, diag) where {T} = Spectrum(energies, states, diag, zero(SVector{0,real(T)}))

function Spectrum(energies, states, diag, basevert::SVector{<:Any,T}) where {T}
    energies´ = maybereal(energies, T)
    subs = collect(approxruns(energies´))
    for rng in subs
        orthonormalize!(view(states, :, rng))
    end
    return Spectrum(energies´, states, diag, basevert, subs)
end

maybereal(energies, ::Type{T}) where {T<:Real} = T.(real.(energies))
maybereal(energies, ::Type{T}) where {T<:Complex} = T.(energies)

function Base.show(io::IO, s::Spectrum{D,C,T}) where {D,C,T}
    i = get(io, :indent, "")
    print(io,
"$(i)Spectrum{$D,$C,$T}: spectrum of a $(D)D Hamiltonian
$i  Orbital type : $C
$i  Energy type  : $T
$i  Energy range : $(extrema(real, s.energies))
$i  Bloch/params : $(Tuple(s.basevert))
$i  Eigenpairs   : $(length(s.energies))
$i  Subspaces    : $(length(s.subs))")
end

"""
    spectrum(h; method = LinearAlgebraPackage(), transform = missing)

Compute the spectrum of a 0D Hamiltonian `h` (or alternatively of the bounded unit cell of a
finite dimensional `h`) using one of the following `method`s

    method                    diagonalization function
    --------------------------------------------------------------
    LinearAlgebraPackage()     LinearAlgebra.eigen!
    ArpackPackage()            Arpack.eigs (must be `using Arpack`)

The option `transform = ε -> f(ε)` allows to transform eigenvalues by `f` in the returned
spectrum (useful for performing shifts or other postprocessing).

# Indexing

The eigenenergies `εv::Vector` and eigenstates `ψm::Matrix` in a `s::Spectrum` object can be
accessed via destructuring, `εv, ψm = sp`, or `εv, ψm = Tuple(sp)`, or `εv = first(sp) =
sp.energies, ψm = last(sp) = sp.states`. Any degenerate energies appear repeated in `εv`.
Alternatively, one can access one or more complete `sub::Subspace`s (eigenenergy together
with its eigenstates, including all degenerates) via the indexing syntax,

    s[1]                   : first `Subspace`
    s[2:4]                 : subspaces 2, 3 and 4
    s[[2,5,6]]             : subspaces 2, 5 and 6
    s[around = 0.2]        : single subspace with energy closest to 0.2
    s[around = (0.2, 10)]  : the ten subspaces with energies closest to 0.2

The eigenenergy `ε` and subspace basis `ψs` of a `sub::Subspace` can be obtained via
destructuring, `ε, ψs = sub`, or `ε = first(sub) = sub.energy, ψs = last(sub) = sub.basis`.
For performance reasons `ψs` is a `SubArray` view of the appropriate columns of `ψm`, not an
independent copy.

# See also
    `bandstructure`, `diagonalizer`
"""
function spectrum(h; method = LinearAlgebraPackage(), transform = missing)
    diag = diagonalizer(h; method = method)
    (ϵk, ψk) = diag((), NoUnflatten())
    s = Spectrum(ϵk, ψk, diag)
    transform === missing || transform!(transform, s)
    return s
end

"""
    transform!(f::Function, s::Spectrum)

Transform the energies of `s` by applying `f` to them in place.
"""
transform!(f, s::Spectrum) = (map!(f, s.energies, s.energies); s)

# destructuring
Base.first(s::Spectrum) = s.energies
Base.last(s::Spectrum) = s.states
Base.iterate(s::Spectrum) = first(s), Val(:states)
Base.iterate(s::Spectrum, ::Val{:states}) = last(s), Val(:done)
Base.iterate(::Spectrum, ::Val{:done}) = nothing
Base.Tuple(s::Spectrum) = (first(s), last(s))

Base.getindex(s::Spectrum, i::Int) = subspace(s, s.subs[i])
Base.getindex(s::Spectrum, is::Union{AbstractUnitRange,AbstractVector}) = getindex.(Ref(s), is)
Base.getindex(s::Spectrum; around) = get_around(s, around)

get_around(s::Spectrum, ε0::Number) = get_around(s, ε0, 1)
get_around(s::Spectrum, (ε0, n)::Tuple) = get_around(s, ε0, 1:n)

function get_around(s::Spectrum, ε0::Number, which)
    rngs = partialsort(s.subs, which, by = rng -> abs(s.energies[first(rng)] - ε0))
    return subspace(s, rngs)
end

subspace(s::Spectrum, rngs) = subspace.(Ref(s), rngs)

function subspace(s::Spectrum, rng::AbstractUnitRange)
    ε = mean(j -> s.energies[j], rng)
    ψs = view(s.states, :, rng)
    return Subspace(s.diag.orbstruct, ε, ψs)
end

nsubspaces(s::Spectrum) = length(s.subs)

basis_slice_type(s::Spectrum) = typeof(view(s.states, :, 1:0))
basis_block_type(s::Spectrum) = typeof(view(s.states, 1:0, 1:0))

######################################################################
# Band
######################################################################
struct Band{D,C,T,S<:AbstractMatrix{C},S´<:AbstractMatrix,D´}  # D´ is dimension of BaseMesh space plus one (energy)
    basemesh::CuboidMesh{D,T,D´}
    verts::Vector{SVector{D´,T}}
    vbases::Vector{S}
    vptrs::Array{UnitRange{Int},D}
    adjmat::SparseMatrixCSC{Bool,Int}
    simps::Vector{NTuple{D´,Int}}
    sbases::Vector{NTuple{D´,S´}}
    sptrs::Array{UnitRange{Int},D´}
end

Band(basemesh::CuboidMesh{D,T}, verts::Vector{SVector{D´,T´}}, args...) where {D,D´,T,T´} =
    Band(CuboidMesh{D,T´}(basemesh), verts, args...)

function Base.show(io::IO, b::Band{D}) where {D}
    i = get(io, :indent, "")
    print(io,
"$(i)Band{$D}: a band over a $(D)D parameter cuboid
$i  Vertices   : $(nvertices(b))
$i  Edges      : $(nedges(b))
$i  Simplices  : $(nsimplices(b))")
end

nvertices(m::Band) = length(m.verts)

nedges(m::Band) = div(nnz(m.adjmat), 2)

nsimplices(m::Band) = length(m.simps)

vertices(m::Band) = m.verts

edges(adjmat, src) = nzrange(adjmat, src)

edgedest(adjmat, edge) = rowvals(adjmat)[edge]

edgevertices(b::Band) =
    ((vsrc, b.verts[edgedest(b.adjmat, edge)]) for (i, vsrc) in enumerate(b.verts) for edge in edges(b.adjmat, i))

degeneracy(m::Band, i) = degeneracy(m.vbases[i])

transform!(f::Function, m::Band) = (map!(f, vertices(m), vertices(m)); m)

Base.getindex(b::Band, i::Int) = b.verts[i], b.vbases[i]

#######################################################################
# Bandstructure
#######################################################################
struct Bandstructure{D,C,T,B<:Band{D,C,T},M<:Diagonalizer}   # D is dimension of base mesh, D´ = D+1
    bands::Vector{B}  # band meshes (vertices + adjacencies)
    diag::M           # diagonalizer that can be used to add additional base-meshes for refinement
end

function Base.show(io::IO, bs::Bandstructure)
    i = get(io, :indent, "")
    ioindent = IOContext(io, :indent => string(i, "  "))
    print(io, i, summary(bs), "\n",
"$i  Bands         : $(length(bs.bands))
$i  Vertices      : $(nvertices(bs))
$i  Edges         : $(nedges(bs))
$i  Simplices     : $(nsimplices(bs))")
end

Base.summary(::Bandstructure{D}) where {D} =
    "Bandstructure{$D}: bands of a $(D)D Hamiltonian"

nvertices(bs::Bandstructure) = isempty(bands(bs)) ? 0 : sum(nvertices, bands(bs))

nedges(bs::Bandstructure) = isempty(bands(bs)) ? 0 : sum(nedges, bands(bs))

nsimplices(bs::Bandstructure) = isempty(bands(bs)) ? 0 : sum(nsimplices, bands(bs))

nbands(bs::Bandstructure) = length(bands(bs))

"""
    bands(bs::Bandstructure[, i])

Return a `bands::Vector{Band}` of all the bands in `bs`, or `bands[i]` if `i` is given.

# See also
    `bandstructure`
"""
bands(bs::Bandstructure) = bs.bands
bands(bs::Bandstructure, i) = bs.bands[i]

"""
    vertices(bs::Bandstructure, i)

Return the vertices `(k..., ϵ)` of the i-th band in `bs`, in the form of a
`Vector{SVector{L+1}}`, where `L` is the lattice dimension.
"""
vertices(bs::Bandstructure, i) = vertices(bands(bs, i))

"""
    transform!(f::Function, b::Bandstructure)

Transform the energies of all bands in `b` by applying `f` to them in place.

    transform!((fk, fε), b::Bandstructure)

Transform Bloch phases and energies of all bands in `b` by applying `fk` and `fε` to them in
place, respectively. If any of them is `missing`, it will be ignored.

"""
transform!(f, bs::Bandstructure) = transform!(sanitize_band_transform(f), bs)

function transform!((fk, fε)::Tuple{Function,Function}, bs::Bandstructure)
    for band in bands(bs)
        vs = vertices(band)
        for (i, v) in enumerate(vs)
            vs[i] = SVector((fk(SVector(Base.front(Tuple(v))))..., fε(last(v))))
        end
    end
    return bs
end

#######################################################################
# isometric and Brillouin zone points
#######################################################################
function isometric(h::Hamiltonian)
    r = qr(bravais(h)).R
    r = r * sign(r[1,1])
    ibr = inv(r')
    return ϕs -> ibr * ϕs
end

isometric(h::Hamiltonian{<:Any,L}, nodes) where {L} = _isometric(h, parsenode.(nodes, Val(L)))

_isometric(h, pts::Tuple) = _isometric(h, [pts...])

function _isometric(h, pts::Vector)
    br = bravais(h)
    pts´ = map(p -> br * p, pts)
    pathlength = pushfirst!(cumsum(norm.(diff(pts))), 0.0)
    isometric = piecewise_mapping(pathlength)
    return isometric
end

nodeindices(nodes::NTuple{N,Any}) where {N} = ntuple(i -> i-1, Val(N))

piecewise_mapping(nodes, ::Val{N}) where {N} = piecewise_mapping(parsenode.(nodes, Val(N)))

function piecewise_mapping(pts)
    N = length(pts) # could be a Tuple or a different container
    mapping = x -> begin
        x´ = clamp(only(x), 0, N-1)
        i = min(floor(Int, x´), N-2) + 1
        p = pts[i] + (x´ - i + 1) * (pts[i+1] - pts[i])
        return p
    end
    return mapping
end

parsenode(pt::SVector, ::Val{L}) where {L} = padright(pt, Val(L))
parsenode(pt::Tuple, val) = parsenode(SVector(float.(pt)), val)

function parsenode(node::Symbol, val)
    pt = get(BZpoints, node, missing)
    pt === missing && throw(ArgumentError("Unknown Brillouin zone point $pt, use one of $(keys(BZpoints))"))
    pt´ = parsenode(pt, val)
    return pt´
end

const BZpoints =
    ( Γ  = (0,)
    , X  = (pi,)
    , Y  = (0, pi)
    , Z  = (0, 0, pi)
    , K  = (2pi/3, -2pi/3)
    , K´ = (4pi/3, 2pi/3)
    , M  = (pi, 0)
    )

#######################################################################
# bandstructure API
#######################################################################
"""
    bandstructure(h::Hamiltonian; subticks = 13, kw...)

Compute `bandstructure(h, cuboid((-π,π)...; subticks = subticks); kw...)` using a base mesh
(of type `CuboidMesh`) over `h`'s full Brillouin zone with the specified `subticks` along
each [-π,π] reciprocal axis.

    bandstructure(h::Hamiltonian, nodes...; subticks = 13, kw...)

Create a linecut of a bandstructure of `h` along a polygonal line connecting two or more
`nodes`. Each node is either a `Tuple` or `SVector` of Bloch phases, or a symbolic name for
a Brillouin zone point (`:Γ`,`:K`, `:K´`, `:M`, `:X`, `:Y` or `:Z`). Each segment in the
polygon has the specified number of `subticks`. Different `subticks` per segments can be
specified with `subticks = (p1, p2...)`.

    bandstructure(h::Hamiltonian, mesh::CuboidMesh; mapping = missing, kw...)

Compute the bandstructure `bandstructure(h, mesh; kw...)` of Bloch Hamiltonian `bloch(h,
ϕ)`, with `ϕ = v` taken on each vertex `v` of the base `mesh` (or `ϕ = mapping(v...)` if a
`mapping` function is provided).

    bandstructure(ph::ParametricHamiltonian, ...; kw...)

Compute the bandstructure of a `ph`. Unless all parameters have default values, a `mapping`
is required between mesh vertices and Bloch/parameters for `ph`, see details on `mapping`
below.

    bandstructure(matrixf::Function, mesh::CuboidMesh; kw...)

Compute the bandstructure of the Hamiltonian matrix `m = matrixf(ϕ)`, with `ϕ` evaluated on
the vertices `v` of the `mesh`. Note that `ϕ` in `matrixf(ϕ)` is an unsplatted container.
Hence, i.e. `matrixf(x) = ...` or `matrixf(x, y) = ...` will not work. Use `matrixf((x,)) =
...`, `matrixf((x, y)) = ...` or matrixf(s::SVector) = ...` instead.

    h |> bandstructure([mesh,]; kw...)

Curried form of the above equivalent to `bandstructure(h[, mesh]; kw...)`.

# Options

The default options are

    (mapping = missing, method = LinearAlgebraPackage(), transform = missing, splitbands = true, showprogress = true)

`mapping`: when not `missing`, `mapping = v -> p` is a function that map base mesh vertices
`v` to Bloch phases and/or parameters `p`. The structure of `p` is whatever is accepted by
`bloch(h, p, ...)` (see `bloch`). For `h::Hamiltonian`, `p = ϕs::Union{Tuple,SVector}` are
Bloch phases. For `h::ParametricHamiltonian`, `p = (ϕs..., (; ps))` or `p = (ϕs, (; ps))`
combine Bloch phases `ϕs` and keyword parameters `ps` of `ph`. This allows to compute a
bandstructure along a cut in the Brillouin zone/parameter space of `ph`, see examples below.

`method`: it is chosen automatically if unspecified, and can be one of the following

    method                     diagonalization function
    --------------------------------------------------------------
    LinearAlgebraPackage()     LinearAlgebra.eigen!
    ArpackPackage()            Arpack.eigs (must be `using Arpack`)

Options passed to the `method` will be forwarded to the diagonalization function. For example,
`method = ArpackPackage(nev = 8, sigma = 1im)` will use `Arpack.eigs(matrix; nev = 8,
sigma = 1im)` to compute the bandstructure.

`transform`: the option `transform = ε -> fε(ε)` allows to transform eigenvalues by `fε` in
the returned bandstructure (useful for performing shifts or other postprocessing). We can
also do `transform -> (fφ, fε)` to transform also mesh vertices with fφ. Additionally,
`transform -> isometric` or `transform -> (isometric, fε)` will transform mesh vertices into
momenta, assuming they represent Bloch phases. This works both in full bandstructures and
linecuts.

`splitbands`: split all bands into disconnected subbands. See also `splitbands!`

`showprogress`: indicate whether progress bars are displayed during the calculation

# Indexing

The bands in a `bs::Bandstructure` object can be accessed with `bands`, while the indexing
syntax `bs[(φs...)]` gives access to one or more `sub::Subspace` objects, contructed by
linear interpolation of each band at base-mesh coordinates `φs`.

    bs[(φs...), 1]                  : first interpolated subspaces at base mesh coordinates `φs`, ordered by energy
    bs[(φs...), 1:3]                : interpolated subspaces 1 to 3 at base mesh coordinates `φs`, ordered by energy
    bs[(φs...)]                     : interpolated subspaces at base mesh coordinates `φs` in any band
    bs[(φs...), around = 0.2]       : the single interpolated subspaces at `φs` with energies closest to 0.2
    bs[(φs...), around = (0.2, 10)] : the ten interpolated subspaces at `φs` with energies closest to 0.2

The eigenenergy `ε` and subspace basis `ψs` of a `sub::Subspace` can themselves be obtained
via destructuring, `ε, ψs = sub`, or `ε = first(sub), ψs = last(sub)`.

# Examples
```jldoctest
julia> h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1)) |> unitcell(3);

julia> bandstructure(h; subticks = 25, method = LinearAlgebraPackage())
Bandstructure{2}: bands of a 2D Hamiltonian
  Bands         : 1
  Vertices      : 9917
  Edges         : 31044
  Simplices     : 20724

julia> bandstructure(h, :Γ, :X, :Y, :Γ; subticks = (10,15,10))
Bandstructure{1}: bands of a 1D Hamiltonian
  Bands         : 1
  Vertices      : 442
  Edges         : 468
  Simplices     : 468

julia> bandstructure(h, cuboid((0, 2π); subticks = 13); mapping = φ -> (φ, 0))
Bandstructure{1}: bands of a 1D Hamiltonian
  Bands         : 1
  Vertices      : 136
  Edges         : 144
  Simplices     : 144

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

julia> bandstructure(ph, cuboid((0, 2π); subticks = 13); mapping = φ -> (φ, 0, (; α = 2φ)))
Bandstructure{1}: bands of a 1D Hamiltonian
  Bands         : 1
  Vertices      : 132
  Edges         : 144
  Simplices     : 144
```

# See also
    `cuboid`, `diagonalizer`, `bloch`, `parametric`, `splitbands!`
"""
bandstructure(args...; kw...) = h -> bandstructure(h, args...; kw...)

function bandstructure(h::Hamiltonian{<:Any, L}; subticks = 13, kw...) where {L}
    L == 0 && throw(ArgumentError("Hamiltonian is 0D, use `spectrum` instead of `bandstructure`"))
    base = cuboid(filltuple((-π, π), Val(L))...; subticks = subticks)
    return bandstructure(h, base; kw...)
end

function bandstructure(h::Hamiltonian{<:Any,L}, node1, node2, nodes...; subticks = 13, transform = missing, kw...) where {L}
    allnodes = (node1, node2, nodes...)
    mapping´ = piecewise_mapping(allnodes, Val(L))
    base = cuboid(nodeindices(allnodes); subticks = subticks)
    transform´ = sanitize_band_transform(transform, h, allnodes)
    return bandstructure(h, base; mapping = mapping´, transform = transform´, kw...)
end

function bandstructure(h::Union{Hamiltonian,ParametricHamiltonian}, basemesh::CuboidMesh;
                       transform = missing, splitbands = true, showprogress = true, kw...)
    diag = diagonalizer(h; kw...)
    b = bandstructure(diag, basemesh, splitbands, showprogress)
    transform === missing || transform!(sanitize_band_transform(transform, h), b)
    return b
end

function bandstructure(matrixf::Function, basemesh::CuboidMesh;
                       mapping = missing, transform = missing, splitbands = true, showprogress = true)
    matrixf´ = wrapmapping(mapping, matrixf)
    dimh = size(samplematrix(matrixf´, basemesh), 1)
    diag = diagonalizer(matrixf´, dimh)
    b = bandstructure(diag, basemesh, splitbands, showprogress)
    transform === missing || transform!(transform, b)
    return b
end

@inline map_phiparams(mapping::Missing, basevertex) = sanitize_phiparams(basevertex)
@inline map_phiparams(mapping::Function, basevertex) = sanitize_phiparams(mapping(basevertex...))

wrapmapping(mapping::Missing, matrixf::Function) = matrixf
wrapmapping(mapping::Function, matrixf::Function) = basevertex -> matrixf(toSVector(mapping(basevertex...)))

sanitize_band_transform(::Missing, args...) = (identity, identity)
sanitize_band_transform(f::Function, args...) = (identity, f)
sanitize_band_transform(f::typeof(isometric), args...) = (isometric(args...), identity)
sanitize_band_transform((_,f)::Tuple{typeof(isometric),Function}, args...) = (isometric(args...), f)
sanitize_band_transform(fs::Tuple{Function,Function}, args...) = fs
sanitize_band_transform((_,f)::Tuple{Missing,Function}, args...) = (identity, f)
sanitize_band_transform((f,_)::Tuple{Function,Missing}, args...) = (f, identity)

samplematrix(matrixf, basemesh) = matrixf(Tuple(first(vertices(basemesh))))

#######################################################################
# bandstructure building
#######################################################################
function bandstructure(diag::Diagonalizer, basemesh::CuboidMesh, split, showprogress)
    # Step 1/2 - Diagonalise:
    spectra = bandstructure_diagonalize(diag, basemesh, showprogress)
    # Step 2/2 - Knit bands:
    bands = bandstructure_knit(basemesh, spectra, showprogress)
    bs = Bandstructure(bands, diag)
    split && (bs = splitbands(bs))
    return bs
end

## Diagonalize bands step
function bandstructure_diagonalize(diag, basemesh::CuboidMesh, showprogress)
    prog = Progress(length(basemesh), "Step 1/2 - Diagonalising: ")
    spectra = [build_spectrum(diag, vertex, showprogress, prog) for vertex in vertices(basemesh)]
    return spectra
end

function build_spectrum(diag, vertex, showprog, prog)
    (ϵs, ψs) = diag(Tuple(vertex), NoUnflatten())
    spectrum = Spectrum(ϵs, ψs, diag, vertex)
    showprog && ProgressMeter.next!(prog)
    return spectrum
end

## Knit bands step
function bandstructure_knit(basemesh::CuboidMesh{D}, spectra::AbstractArray{SP}, showprog) where {D,C,T,SP<:Spectrum{D,C,T}}
    simpitr = marchingsimplices(basemesh) # D+1 - dimensional iterator over simplices

    S       = basis_slice_type(first(spectra))
    verts   = SVector{D+1,T}[]
    vbases  = S[]
    vptrs   = Array{UnitRange{Int}}(undef, size(vertices(basemesh)))
    simps   = NTuple{D+1,Int}[]
    sbases  = NTuple{D+1,Matrix{C}}[]
    sptrs   = Array{UnitRange{Int}}(undef, size(simpitr))

    cbase = eachindex(basemesh)
    lbase = LinearIndices(cbase)

    # Collect vertices
    for csrc in cbase
        len = length(verts)
        push_verts!((verts, vbases), spectra[csrc])
        vptrs[csrc] = len+1:length(verts)
    end

    # Store subspace projections in vertex adjacency matrix
    prog = Progress(length(cbase), "Step 2/2 - Knitting bands: ")
    S´ = basis_block_type(first(spectra))
    I, J, V = Int[], Int[], S´[]
    for csrc in cbase
        for cdst in neighbors_forward(basemesh, csrc)
            push_adjs!((I, J, V), spectra[csrc], spectra[cdst], vptrs[csrc], vptrs[cdst])
        end
        showprog && ProgressMeter.next!(prog)
    end
    n = length(verts)
    adjprojs = sparse(I, J, V, n, n)
    adjmat = sparsealiasbool(adjprojs)

    # Build simplices from stable cycles around base simplices
    buffers = NTuple{D+1,Int}[], NTuple{D+1,Int}[]
    emptybases = filltuple(fill(zero(C), 0, 0), Val(D+1))  # sentinel bases for deg == 1 simplices
    for (csimp, vs) in zip(CartesianIndices(simpitr), simpitr)  # vs isa NTuple{D+1,CartesianIndex{D}}
        len = length(simps)
        ranges = getindex.(Ref(vptrs), vs)
        push_simps!((simps, sbases, verts), buffers, ranges, adjprojs, emptybases)
        sptrs[csimp] = len+1:length(simps)
    end

    bands = [Band(basemesh, verts, vbases, vptrs, adjmat, simps, sbases, sptrs)]
    return bands
end

function push_verts!((verts, vbases), spectrum)
    for rng in spectrum.subs
        vert = SVector(Tuple(spectrum.basevert)..., spectrum.energies[first(rng)])
        basis = view(spectrum.states, :, rng)
        push!(verts, vert)
        push!(vbases, basis)
    end
    return nothing
end

function push_adjs!((I, J, V), ssrc, sdst, rsrc, rdst)
    ψdst = sdst.states
    ψsrc = ssrc.states
    proj = ψdst' * ψsrc
    proj´ = copy(proj')
    for (is, rs) in enumerate(ssrc.subs)
        srcdim = length(rs)
        for (id, rd) in enumerate(sdst.subs)
            crange = CartesianIndices((rd, rs))
            crange´ = CartesianIndices((rs, rd))
            rank = rankproj(proj, crange)
            if !iszero(rank)
                srcdim -= rank
                srcdim < 0 && @warn("Unexpected band connectivity between $(ssrc.basevert) and $(sdst.basevert). Rank $rank in $(size(crange)) projector.")
                append!(I, (rdst[id], rsrc[is]))
                append!(J, (rsrc[is], rdst[id]))
                append!(V, (view(proj, crange), view(proj´, crange´)))
            end
        end
    end
    return nothing
end

# equivalent to r = round(Int, tr(m'm)) over crange, but if r > 0, must compute and count singular values
function rankproj(m, crange)
    r = round(Int, sum(c -> abs2(m[c]), crange))
    (iszero(r) || minimum(size(crange)) == 1) && return r
    s = svdvals(view(m, crange))
    r = count(s -> abs2(s) >= 0.5, s)
    return r
end

sparsealiasbool(s) = SparseMatrixCSC(s.m, s.n, s.colptr, s.rowval, fill(true, length(s.nzval)))

function push_simps!((simps, sbases, verts), buffers, ranges, adjprojs, emptybases)
    cycles = build_cycles!(buffers, ranges, adjprojs)
    for cycle in cycles
        done = false
        projs = getindex.(Ref(adjprojs), shiftleft(cycle), cycle)
        degs = degeneracy.(projs)
        if all(isequal(1), degs)  # optimization with sentinel emptybases
            push!(simps, cycle)
            push!(sbases, emptybases)
            continue
        end
        d, start = findmin(degs) # minimum degeneracy of source
        projs = shiftleft(projs, start - 1)
        bases = shiftright(accumulate(orthomul, projs; init = I(d)))
        while !done
            if rankbasis(first(bases)) < d
                basis_start = discard_zero_cols(first(bases))
                d = rankbasis(basis_start)
                d == 0 && break
                bases = shiftright(accumulate(orthomul, projs; init = basis_start))
            else
                done = true
            end
        end
        if done
            push!(simps, cycle)
            bases = shiftright(bases, start - 1)
            push!(sbases, bases)
        end
    end
    return nothing
end

# build cycles of elements in ranges that are connected by adjprojs
function build_cycles!((partials, partials´), ranges::NTuple{D´}, adjprojs) where {D´}
    resize!(partials, length(first(ranges)))
    for (n, i) in enumerate(first(ranges))
        partials[n] = padright((i,), Val(D´))
    end
    rows = rowvals(adjprojs)
    for n in 1:D´
        empty!(partials´)
        for partial in partials
            if n == D´
                seed = first(partial)
                nextrng = seed:seed
            else
                nextrng = ranges[n+1]
            end
            src = partial[n]
            for ptr in nzrange(adjprojs, src)
                dst = rows[ptr]
                dst in nextrng || continue
                partial´ = n == D´ ? partial : tuplesplice(partial, n+1, dst)
                push!(partials´, partial´)
            end
        end
        partials, partials´ = partials´, partials
    end
    return partials
end

function discard_zero_cols(mat)
    ncols = rankbasis(mat)
    mat´ = similar(mat, size(mat, 1), ncols)
    j = 0
    for col in eachcol(mat)
        if !iszero(col)
            j += 1
            mat´[:, j] .= col
        end
    end
    return mat´
end

# cheap-rank for number of non-zero columns
rankbasis(m) = count(!iszero, eachcol(m))

# The minimum projected norm2 per src state (m column) is 0.5/dim(src) by default.
# 0.5 due to adjacency, and 1/dim(src) because the projection eigenstate has a minimum
# 1/√din(src) component over some src basis vector
function orthomul(m1, m2)
    m = m2 * m1
    threshold = 0.5 / size(m, 2)
    return orthonormalize!(m, threshold)
end

#######################################################################
# splitbands
#######################################################################
"""
    splitbands(bs::Bandstructure)

Splits the bands in `bs` into disconnected subbands that share no vertices. See also `splitbands` option in `bandstructure`.

# See also
    `bandstructure`
"""
splitbands(b::Bandstructure) = Bandstructure(splitbands(b.bands), b.diag)

function splitbands(bands::Vector{<:Band})
    bands´ = similar(bands, 0)
    splitbands!.(Ref(bands´), bands)
    return bands´
end

function splitbands!(bands, band)
    bs, is, nbands = subgraphs(band.adjmat)
    nbands == 1 && return push!(bands, band)
    adjmats = split_adjmats(band.adjmat, bs, is, nbands)
    verts, vbases, vptrs = split_list(band.verts, band.vbases, band.vptrs, bs, is, nbands)
    simps, sbases, sptrs = split_list(band.simps, band.sbases, band.sptrs, bs, is, nbands)
    for n in 1:nbands
        newband = Band(band.basemesh, verts[n], vbases[n], vptrs[n], adjmats[n], simps[n], sbases[n], sptrs[n])
        nsimplices(newband) > 0 && push!(bands, newband)
    end
    return bands
end

struct BandIndex
    n::Int
end

Base.:*(_, b::BandIndex) = b
Base.:*(b::BandIndex, _) = b
Base.:+(b1::BandIndex, b2::BandIndex) = BandIndex(min(b1.n, b2.n))
Base.zero(b::Type{BandIndex}) = BandIndex(typemax(Int))

# for each row or column i in a, compute its subgraph bs[i] and the column/row index within that subgraph, is[i]
function subgraphs(a)
    bs = BandIndex.(1:size(a, 2))
    bs´ = copy(bs)
    total = sum(i -> i.n, bs)
    done = false
    while !done
        mul!(bs´, a, bs, true, true)
        bs, bs´ = bs´, bs
        total´ = sum(i -> i.n, bs)
        done = total == total´
        total = total´
    end
    dict = Dict(1 => 1)
    m = 1
    for b in bs
        if !haskey(dict, b.n)
            m += 1
            dict[b.n] = m
        end
    end
    is = similar(bs, Int)
    ngraphs = m
    counters = fill(0, ngraphs)
    for (i, b) in enumerate(bs)
        n = dict[b.n]
        bs[i] = BandIndex(n)
        counters[n] += 1
        is[i] = counters[n]
    end
    return bs, is, ngraphs
end

function split_adjmats(adjmat, bs, is, nbands)
    I0, J0, _ = findnz(adjmat)
    Is = [Int[] for _ in 1:nbands]
    Js = copy.(Is)
    for (i, j) in zip(I0, J0)
        n = bs[i].n
        n == bs[j].n || throw(ArgumentError("Unexpected error in subgraphs function"))
        push!(Is[n], is[i])
        push!(Js[n], is[j])
    end
    nverts = [count(isequal(BandIndex(n)), bs) for n in 1:nbands]
    adjmats = [sparse(Is[n], Js[n], true, nverts[n], nverts[n]) for n in 1:nbands]
    return adjmats
end

function split_list(vs, vbs, vptrs, bs, is, nbands)
    vs´ = [similar(vs, 0) for _ in 1:nbands]
    vbs´ = [similar(vbs, 0) for _ in 1:nbands]
    vptrs´ = [similar(vptrs) for _ in 1:nbands]
    counters = fill(0, nbands)
    counters´ = copy(counters)
    for (iptr, rng) in enumerate(vptrs)
        for i in rng
            n, v = get_vert_or_simp(vs[i], i, bs, is)
            push!(vs´[n], v)
            push!(vbs´[n], vbs[i])
        end
        counters´ = length.(vs´)
        for n in 1:nbands
            vptrs´[n][iptr] = counters[n]+1:counters´[n]
        end
        counters, counters´ = counters´, counters
    end
    return vs´, vbs´, vptrs´
end

# return band index and vertex/simplex
get_vert_or_simp(v::SVector, i, bs, is) = bs[i].n, v
get_vert_or_simp(s::NTuple, i, bs, is) = bs[first(s)].n, getindex.(Ref(is), s)

#######################################################################
# Bandstructure indexing
#######################################################################
Base.getindex(bs::Bandstructure, ϕs::Tuple; around = missing) = interpolate_bandstructure(bs, ϕs, around)

function interpolate_bandstructure(bs::Bandstructure, ϕs, around)
    subs = subspace_type(bs, ϕs)[]
    foreach(band -> interpolate_band!(subs, band, ϕs, bs.diag.orbstruct), bs.bands)
    return filter_around!(subs, around)
end

function subspace_type(bs::Bandstructure, ϕs)
    (isempty(bs.bands) || isempty(first(bs.bands).vbases)) && throw(ArgumentError("Cannot index into empty bandstructure"))
    band = first(bs.bands)
    ε0 = last(first(band.verts))
    b0 = first(band.vbases)
    m0 = Matrix{eltype(b0)}(undef, size(b0, 1), 0)
    return typeof(Subspace(bs.diag.orbstruct, ε0, m0, ϕs))
end

function interpolate_band!(subs, band::Band{D}, ϕs, orbstruct) where {D}
    D == length(ϕs) || throw(ArgumentError("Bandstructure needs a NTuple{$D} of base coordinates for interpolation"))
    s = find_basemesh_simplex(promote(ϕs...), band.basemesh)
    s === nothing && return subs
    (ϕs, inds, dϕinds) = s
    for i in band.sptrs[inds...]
        si = band.simps[i]
        ε0, εs = firsttail(last.(getindex.(Ref(band.verts), si)))
        b0, bs = firsttail(getindex.(Ref(band.vbases), si))
        p0, ps = firsttail(band.sbases[i])
        energy = ε0 * (1 - sum(dϕinds)) + sum(εs .* dϕinds)
        basis = isempty(p0) ? copy(b0) : b0 * p0
        lmul!(1 - sum(dϕinds), basis)
        if isempty(p0)  # deg == 1 sentinel optimization
            for j in 1:D
                basis .+= bs[j] .* (sentinel_phase(bs[j], b0) * dϕinds[j])
            end
        else
            foreach(i -> mul!(basis, bs[i], ps[i], dϕinds[i], true), 1:D)
        end
        # h is necessary here to perhaps unflatten basis
        push!(subs, Subspace(orbstruct, energy, basis, ϕs))
    end
    return subs
end

sentinel_phase(bj::AbstractVector{<:Real}, b0::AbstractVector{<:Real}) = 1
sentinel_phase(bj, b0) = cis(angle(dot(bj, b0)))

# The simplex sought has unitperm such that reverse(dϕs[unitperm]) is sorted
function find_basemesh_simplex(ϕs::NTuple{D}, basemesh) where {D}
    baseinds = find_basemesh_interval.(ϕs, basemesh.ticks)
    any(iszero, baseinds) && return nothing
    vertex0 = getindex.(basemesh.ticks, baseinds)
    vertex1 = getindex.(basemesh.ticks, baseinds .+ 1)
    dϕs = (ϕs .- vertex0) ./ (vertex1 .- vertex0)
    dϕperm = Tuple(sortperm(SVector(dϕs), rev = true))
    baseperms = unitperms(basemesh)
    i = findfirst(isequal(dϕperm), baseperms)
    i === nothing && return nothing
    sverts = SVector.(Tuple.(Base.tail(unitsimplex(basemesh, i))))
    smat = hcat(sverts...)
    dϕinds = inv(smat) * SVector(dϕs)
    inds = (baseinds..., i)
    return ϕs, inds, dϕinds
end

function find_basemesh_interval(ϕ, ticks)
    @inbounds for m = 1:length(ticks)-1
        ticks[m] <= ϕ < ticks[m+1] && return m
    end
    ϕ ≈ last(ticks) && return length(ticks) - 1
    return 0
end

filter_around!(ss, ::Missing) = sort!(ss; by = s -> s.energy, alg = Base.DEFAULT_UNSTABLE)
filter_around!(ss, ε0::Number) = filter_around!(ss, ε0, 1)

function filter_around!(ss, (ε0, n)::Tuple)
    0 <= n <= length(ss) || throw(ArgumentError("Cannot retrieve more than $(length(ss)) subspaces"))
    return filter_around!(ss, ε0, 1:n)
end

filter_around!(ss::Vector{<:Subspace}, ε0, which) = partialsort!(ss, which; by = s -> abs(s.energy - ε0))

#######################################################################
# Bandstructure minima, maxima, gap, gapedge
#######################################################################
"""
    minima(b::Bandstructure{1}; refinesteps = 0)

For a 1D bandstructure `b`, compute a vector of `Vector{Tuple{T,T}}`s (one per band),
containing pairs `(φ, ε)` of Bloch phase and energy where the band, as sampled, has a local
minimum. The minima will be further refined by a number `refinesteps` of bisections steps.
Only band vertices with one neighbors on each side will be considered as potential local
minimum.

# See also:
    `maxima`, `gapedge`, `gap`
"""
minima(b::Bandstructure{1,<:Any,T}; kw...) where {T} =
    Vector{Tuple{real(T),real(T)}}[band_extrema(band, minimum_criterion, b.diag; kw...) for band in b.bands]

minimum_criterion(ε, εs...) = all(>=(ε), εs)

"""
    maxima(b::Bandstructure{1}; refinesteps = 0)

For a 1D bandstructure `b`, compute a vector of `Vector{Tuple{T,T}}`s (one per band),
containing pairs `(φ, ε)` of Bloch phase and energy where the band, as sampled, has a local
maximum. The maxima will be further refined by a number `refinesteps` of bisections steps.
Only band vertices with one neighbors on each side will be considered as potential local
maximum.

# See also:
    `minima`, `gapedge`, `gap`
"""
maxima(b::Bandstructure{1,<:Any,T}; kw...) where {T} =
    Vector{Tuple{real(T),real(T)}}[band_extrema(band, maximum_criterion, b.diag; kw...) for band in b.bands]

maximum_criterion(ε, εs...) = all(<=(ε), εs)

function band_extrema(b::Band{1}, criterion, diag; refinesteps = 0)
    found, neighbors = findall_with_neighbors(criterion, b)
    found´ = refine_bisection.(found, neighbors, refinesteps, criterion, Ref(diag))
    unique!(x -> round.(chop.(x), digits = 8), found´)
    return found´
end

function findall_with_neighbors(criterion::Function, b::Band{D,<:Any,T´}) where {D,T´}
    vertices = b.verts
    T = real(T´)
    found = NTuple{2,T}[]
    neighs = Tuple{NTuple{2,T},NTuple{2,T}}[]
    for (i, vertex) in enumerate(vertices)
        ns = neighbors(b, i)
        if length(ns) == max_neighbors(D)
            (φ0, ε0) = real.(vertex)
            i1, i2 = ns
            (φ1, ε1) = real.(vertices[i1])
            (φ2, ε2) = real.(vertices[i2])
            φ1 < φ0 < φ2 || φ2 < φ0 < φ1 || continue
            if criterion(ε0, ε1, ε2)
                push!(found, (φ0,  ε0))
                push!(neighs, ((φ1, ε1), (φ2, ε2)))
            end
        end
    end
    return found, neighs
end

neighbors(b::Band, i::Int) = (rowvals(b.adjmat)[ptr] for ptr in nzrange(b.adjmat, i))

max_neighbors(D::Int) = 2^D # sum(m->binomial(D+1, m), 1:D), where D = 1 for 1D bands

function refine_bisection(found, neighs, steps, criterion, diag)
    φ0, ε0 = found
    (φ1, ε1), (φ2, ε2) = neighs
    if φ1 > φ2
        φ1, ε1, φ2, ε2 = φ2, ε2, φ1, ε1
    end
    φ1 < φ0 < φ2 || return φ0, ε0
    for _ in 1:steps
        a, b = SA[(φ1-φ0)^2 φ1-φ0; (φ2-φ0)^2 φ2-φ0] \ SA[ε1-ε0; ε2-ε0]
        dφ = -b / 2a  # corresponds to estimated dε = -b^2/4a
        if iszero(dφ) || isnan(dφ) || ifelse(dφ > 0, dφ/(φ2 - φ0), -dφ/(φ0 - φ1)) > 1
            break
        end
        εs, _ = diag(φ0 + dφ)
        ε0´ = select_closest(real.(εs), criterion, ε0, ε1, ε2)
        if dφ > 0
            φ1, φ0, φ2 = φ0, φ0 + dφ, φ2
            ε1, ε0, ε2 = ε0, ε0´, ε2
        else
            φ1, φ0, φ2 = φ1, φ0 + dφ, φ0
            ε1, ε0, ε2 = ε1, ε0´, ε0
        end
    end
    return chop(φ0), ε0
end

select_closest(εs, criterion, ε0, εi...) =
    argmin(ε -> ifelse(criterion(ε, ε0, εi...), abs(ε - ε0), Inf), εs)

"""
    gapedge(b::Bandstructure{1}, ε₀; refinesteps = 0)

For a 1D bandstructure `b`, compute two tuples, `(φ₊, ε₊)` and `(φ₋, ε₋)`, of band points
closest in energy to `ε₀`, from above and below, respectively. If ε₀ is inside a band or
outside the global bandwidth, `φ₊` and `φ₋` will be `missing`. See `minima` or `maxima` for
details about `refinesteps`.

    gapedge(b::Bandstructure{1}, ε₀, +; kw...)
    gapedge(b::Bandstructure{1}, ε₀, -; kw...)

Compute only `(φ₊, ε₊)` or `(φ₋, ε₋)`, respectively.

# See also:
    `gap`, `minima`, `maxima`
"""
gapedge(b::Bandstructure{1}, ε0; kw...) = gapedge(b, ε0, +; kw...), gapedge(b, ε0, -; kw...)

function gapedge(b::Bandstructure{1,<:Any,T}, ε0, ::typeof(+); kw...) where {T}
    isinband(b, ε0) && return (missing, zero(T))
    minbands = Iterators.flatten(filter!.(φε -> last(φε) > ε0, minima(b; kw...)))
    isempty(minbands) && return (missing, T(Inf))
    (φ₊, ε₊) = argmin(last, minbands)
    return (φ₊, ε₊)
end

function gapedge(b::Bandstructure{1,<:Any,T}, ε0, ::typeof(-); kw...) where {T}
    isinband(b, ε0) && return (missing, zero(T))
    maxbands = Iterators.flatten(filter!.(φε -> last(φε) < ε0, maxima(b; kw...)))
    isempty(maxbands) && return (missing, T(-Inf))
    (φ₋, ε₋) = argmax(last, maxbands)
    return (φ₋, ε₋)
end

"""
    isinband(b::Bandstructure, ε)
    isinband(b::Band, ε)

Returns true if `ε` is contained within a band, false otherwise.

# See also:
    `gap`, `minima`, `maxima`
"""
isinband(b::Bandstructure, ε) = any(band -> isinband(band, ε), b.bands)
isinband(b::Band, ε) = maximum(last, b.verts) > ε > minimum(last, b.verts)

"""
    gap(b::Bandstructure{1}, ε₀; refinesteps = 0)

Compute the gap if a 1D bandstructure `b` around ε₀, if any.

# See also:
    `gapedge`, `minima`, `maxima`
"""
function gap(b::Bandstructure{1}, ε0; kw...)
    (_, ε₊), (_, ε₋) = gapedge(b, ε0; kw...)
    return ε₊ - ε₋
end


text:
# compiler driver and main interface

# (::CompilerJob)
const compile_hook = Ref{Union{Nothing,Function}}(nothing)

"""
    compile(target::Symbol, cap::VersionNumber, f, tt, kernel=true;
            libraries=true, dynamic_parallelism=true,
            optimize=true, strip=false, strict=true, ...)

Compile a function `f` invoked with types `tt` for device capability `cap` to one of the
following formats as specified by the `target` argument: `:julia` for Julia IR, `:llvm` for
LLVM IR, `:ptx` for PTX assembly and `:cuda` for CUDA driver objects. If the `kernel` flag
is set, specialized code generation and optimization for kernel functions is enabled.

The following keyword arguments are supported:
- `libraries`: link the CUDAnative runtime and `libdevice` libraries (if required)
- `dynamic_parallelism`: resolve dynamic parallelism (if required)
- `optimize`: optimize the code (default: true)
- `strip`: strip non-functional metadata and debug information  (default: false)
- `strict`: perform code validation either as early or as late as possible

Other keyword arguments can be found in the documentation of [`cufunction`](@ref).
"""
compile(target::Symbol, cap::VersionNumber, @nospecialize(f::Core.Function),
        @nospecialize(tt), kernel::Bool=true; libraries::Bool=true,
        dynamic_parallelism::Bool=true, optimize::Bool=true,
        strip::Bool=false, strict::Bool=true, kwargs...) =
    compile(target, CompilerJob(f, tt, cap, kernel; kwargs...);
            libraries=libraries, dynamic_parallelism=dynamic_parallelism,
            optimize=optimize, strip=strip, strict=strict)

function compile(target::Symbol, job::CompilerJob;
                 libraries::Bool=true, dynamic_parallelism::Bool=true,
                 optimize::Bool=true, strip::Bool=false, strict::Bool=true)
    @debug "(Re)compiling function" job

    if compile_hook[] != nothing
        global globalUnique
        previous_globalUnique = globalUnique

        compile_hook[](job)

        globalUnique = previous_globalUnique
    end

    return codegen(target, job;
                   libraries=libraries, dynamic_parallelism=dynamic_parallelism,
                   optimize=optimize, strip=strip, strict=strict)
end

function codegen(target::Symbol, job::CompilerJob;
                 libraries::Bool=true, dynamic_parallelism::Bool=true, optimize::Bool=true,
                 strip::Bool=false,strict::Bool=true)
    ## Julia IR

    @timeit_debug to "validation" check_method(job)

    @timeit_debug to "Julia front-end" begin

        # get the method instance
        world = typemax(UInt)
        meth = which(job.f, job.tt)
        sig = Base.signature_type(job.f, job.tt)::Type
        (ti, env) = ccall(:jl_type_intersection_with_env, Any,
                          (Any, Any), sig, meth.sig)::Core.SimpleVector
        if VERSION >= v"1.2.0-DEV.320"
            meth = Base.func_for_method_checked(meth, ti, env)
        else
            meth = Base.func_for_method_checked(meth, ti)
        end
        method_instance = ccall(:jl_specializations_get_linfo, Ref{Core.MethodInstance},
                      (Any, Any, Any, UInt), meth, ti, env, world)

        for var in env
            if var isa TypeVar
                throw(KernelError(job, "method captures a typevar (you probably use an unbound type variable)"))
            end
        end
    end

    target == :julia && return method_instance

    ## LLVM IR

    defs(mod)  = filter(f -> !isdeclaration(f), collect(functions(mod)))
    decls(mod) = filter(f ->  isdeclaration(f) && intrinsic_id(f) == 0,
                        collect(functions(mod)))

    # always preload the runtime, and do so early; it cannot be part of any timing block
    # because it recurses into the compiler
    if libraries
        runtime = load_runtime(job.cap)
        runtime_fns = LLVM.name.(defs(runtime))
    end

    @timeit_debug to "LLVM middle-end" begin
        ir, kernel = @timeit_debug to "IR generation" irgen(job, method_instance, world)

        if libraries
            undefined_fns = LLVM.name.(decls(ir))
            if any(fn->startswith(fn, "__nv_"), undefined_fns)
                libdevice = load_libdevice(job.cap)
                @timeit_debug to "device library" link_libdevice!(job, ir, libdevice)
            end
        end

        if optimize
            kernel = @timeit_debug to "optimization" optimize!(job, ir, kernel)
        end

        if libraries
            undefined_fns = LLVM.name.(decls(ir))
            if any(fn -> fn in runtime_fns, undefined_fns)
                @timeit_debug to "runtime library" link_library!(job, ir, runtime)
            end
        end

        if ccall(:jl_is_debugbuild, Cint, ()) == 1
            @timeit_debug to "verification" verify(ir)
        end

        kernel_fn = LLVM.name(kernel)
    end

    # dynamic parallelism
    if dynamic_parallelism && haskey(functions(ir), "cudanativeCompileKernel")
        dyn_marker = functions(ir)["cudanativeCompileKernel"]

        cache = Dict{CompilerJob, String}(job => kernel_fn)

        # iterative compilation (non-recursive)
        changed = true
        while changed
            changed = false

            # find dynamic kernel invocations
            # TODO: recover this information earlier, from the Julia IR
            worklist = MultiDict{CompilerJob, LLVM.CallInst}()
            for use in uses(dyn_marker)
                # decode the call
                call = user(use)::LLVM.CallInst
                id = convert(Int, first(operands(call)))

                global delayed_cufunctions
                dyn_f, dyn_tt = delayed_cufunctions[id]
                dyn_job = CompilerJob(dyn_f, dyn_tt, job.cap, #=kernel=# true)
                push!(worklist, dyn_job => call)
            end

            # compile and link
            for dyn_job in keys(worklist)
                # cached compilation
                dyn_kernel_fn = get!(cache, dyn_job) do
                    dyn_ir, dyn_kernel = codegen(:llvm, dyn_job;
                                                 optimize=optimize, strip=strip,
                                                 dynamic_parallelism=false, strict=false)
                    dyn_kernel_fn = LLVM.name(dyn_kernel)
                    link!(ir, dyn_ir)
                    changed = true
                    dyn_kernel_fn
                end
                dyn_kernel = functions(ir)[dyn_kernel_fn]

                # insert a pointer to the function everywhere the kernel is used
                T_ptr = convert(LLVMType, Ptr{Cvoid})
                for call in worklist[dyn_job]
                    Builder(JuliaContext()) do builder
                        position!(builder, call)
                        fptr = ptrtoint!(builder, dyn_kernel, T_ptr)
                        replace_uses!(call, fptr)
                    end
                    unsafe_delete!(LLVM.parent(call), call)
                end
            end
        end

        # all dynamic launches should have been resolved
        @compiler_assert isempty(uses(dyn_marker)) job
        unsafe_delete!(ir, dyn_marker)
    end

    if strict
        # NOTE: keep in sync with non-strict check below
        @timeit_debug to "validation" begin
            check_invocation(job, kernel)
            check_ir(job, ir)
        end
    end

    if strip
        @timeit_debug to "strip debug info" strip_debuginfo!(ir)
    end

    target == :llvm && return ir, kernel

    ## PTX machine code

    @timeit_debug to "LLVM back-end" begin
        @timeit_debug to "preparation" prepare_execution!(job, ir)

        asm = @timeit_debug to "machine-code generation" mcgen(job, ir, kernel)
    end

    undefined_fns = LLVM.name.(decls(ir))
    target == :ptx && return asm, kernel_fn, undefined_fns

    error("Unknown compilation target $target")
end

ext: jl


text:
#!/usr/bin/env julia

using Getopt

include("rsa.jl")

function main()
    infile = stdin
    outfile = stdout
    pubkey = open("/tmp/rsa.pub")

    for (opt, arg) in getopt(ARGS, "i:o:k:", ["infile=", "outfile=", "pubkey="])
        if opt == "-i" || opt == "--infile"
            infile = open(arg, "r")
        elseif opt == "-o" || opt == "--outfile"
            outfile = open(arg, "w")
        elseif opt == "-k" || opt == "--pubkey"
            pubkey = open(arg, "r")
        end
    end

    # Public key: (n, e)
    n = parse(BigInt, readline(pubkey))
    e = parse(BigInt, readline(pubkey))

    block = (floor(UInt64, log2(n)) - 1) ÷ 8
    while !eof(infile)
        bytes = prepend!(read(infile, block - 1), [0xFF])
        m = parse(BigInt, "0x" * bytes2hex(bytes))
        c = encrypt(m, e, n)
        println(outfile, c)
    end

    close(infile)
    close(outfile)
    close(pubkey)
end

main()

max_stars_repo_name: darrelllong/RSA-in-Julia

max_issues_repo_name: darrelllong/RSA-in-Julia

max_forks_repo_name: darrelllong/RSA-in-Julia


text:
module TestMaps

using Test
using ThreadPools

include("util.jl")

@testset "maps" begin

    @testset "tmap" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        primary = Threads.nthreads() == 1
        fn! = (x) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data
        end
        @test tmap(fn!, objs) == collect(1:N)
        @test primary
        @inferred tmap(fn!, objs)
    end

    @testset "bmap" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        fn! = (x) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data
        end
        @test bmap(fn!, objs) == collect(1:N)
        @inferred bmap(fn!, objs)
    end

    @testset "qmap" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        primary = Threads.nthreads() == 1
        fn! = (x) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data
        end
        @test qmap(fn!, objs) == collect(1:N)
        @test primary
        @inferred qmap(fn!, objs)
    end

    @testset "qbmap" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        fn! = (x) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data
        end
        @test qbmap(fn!, objs) == collect(1:N)
        @inferred qbmap(fn!, objs)
    end

    @testset "logtmap" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        primary = Threads.nthreads() == 1
        fn! = (x) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data
        end
        p, r = logtmap(fn!, objs)
        @test r == collect(1:N)
        @test length(p.recs) == N*2
        @test primary
        @inferred logtmap(fn!, objs)
    end

    @testset "logbmap" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        fn! = (x) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data
        end
        p, r = logbmap(fn!, objs)
        @test r == collect(1:N)
        @test length(p.recs) == N*2
        @inferred logbmap(fn!, objs)
    end

    @testset "logqmap" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        primary = Threads.nthreads() == 1
        fn! = (x) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data
        end
        p, r = logqmap(fn!, objs)
        @test r == collect(1:N)
        @test length(p.recs) == N*2
        @test primary
        @inferred logqmap(fn!, objs)
    end

    @testset "logqbmap" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        fn! = (x) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data
        end
        p, r = logqbmap(fn!, objs)
        @test r == collect(1:N)
        @test length(p.recs) == N*2
        @inferred logqbmap(fn!, objs)
    end

end

end # module


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module Fibonacci

using Nemo, IterTools

export ModuleFibonacci
export I000045, F000045, L000045, V000045, R000045, is000045, L039834

"""

* I000045, F000045, L000045, V000045, R000045, is000045
"""
const ModuleFibonacci = ""

struct FiboIterate
    length
end

function Base.iterate(I::FiboIterate)
    I.length == 0 && return nothing
    state = (fmpz(0), (fmpz(0), fmpz(1), fmpz(1)))
end

function Base.iterate(I::FiboIterate, (a, b, c))
    c >= I.length && return nothing
    (b, (b, a + b, c + 1))
end

Base.length(f::FiboIterate) = f.length
Base.eltype(f::FiboIterate) = fmpz

"""

Iterate over the first ``n`` Fibonacci numbers.
"""
I000045(n) = FiboIterate(n)

"""

Iterate over the Fibonacci numbers which do not exceed ``n``.
"""
F000045(n) = IterTools.takewhile(k -> k <= n, FiboIterate(n + 1))

"""

Return the first ``n`` Fibonacci numbers in an array.
"""
L000045(n) = Base.collect(FiboIterate(n))

"""

Return the ``n``-th Fibonacci number.
"""
function V000045(n)
    F = fmpz[1 1; 1 0]
    Fn = F^n
    Fn[2, 1]
end

"""

Return the ``n``-th Fibonacci number, explicite formula by Paul Hankin.
"""
function HankinFibonacci(n)
    n < 2 && return n
    a = <<(fmpz(4), (n - 1) * (n + 2))
    b = <<(fmpz(4), 2 * (n - 1))
    c = <<(fmpz(2), n - 1)
    div(a, b - c - 1) & (c - 1)
end

"""

Fibonacci function for real values, returns a Float64.
"""
function R000045(x::Float64)
    (Base.MathConstants.golden^x -
     cos(x * Base.MathConstants.pi) * Base.MathConstants.golden^(-x)) / sqrt(5)
end

"""

Query if ``n`` is a Fibonacci number, returns a Bool.
"""
function is000045(n)
    d = 0
    for f ∈ FiboIterate(n + 2)
        d = n - f
        d <= 0 && break
    end
    d == 0
end

"""

Return the signed Fibonacci numbers, dubbed negaFibonacci numbers.
"""
function L039834(n)
    L = Array{fmpz,1}(undef, n)
    x, y = 1, 1
    for n in 1:n
        L[n] = x
        x, y = y, x - y
    end
    L
end

#START-TEST-########################################################

using Test

function test()
    @testset "Fibonacci" begin
        @test isa(V000045(30), fmpz)
    end
end

function demo()
    println("V000045")
    for n ∈ 0:12
        print(V000045(n), ", ")
    end
    println()

    println("I000045")
    for f ∈ I000045(20)
        print(f, ", ")
    end
    println()

    println("F000045")
    for f ∈ F000045(20)
        print(f, ", ")
    end
    println()

    println("L000045")
    print(L000045(17))
    println()

    println("R000045")
    println([R000045(Float64(x / 2 + 0.5)) for x ∈ 0:9])
    println(typeof(R000045(2 + 0.5)))

    println("is000045")
    for n ∈ 0:150
        is000045(n) && print(n, " ")
    end
    println()
end

"""

for n ∈ 1:1000 V000045(n) end
    0.004861 seconds (8.49 k allocations: 320.141 KiB)
for fib ∈ I000045(10000) end
    0.006211 seconds (40.00 k allocations: 937.547 KiB)
L000045(10000)
    0.006533 seconds (49.49 k allocations: 1.137 MiB)
"""
function perf()
    GC.gc()
    @time (for n ∈ 1:1000
        V000045(n)
    end)
    @time (for fib ∈ I000045(10000)
    end)
    @time L000045(10000)
    # println("----")
    # @time (for n ∈ 1:1000 HankinFibonacci(n) end)
end

function main()
    test()
    demo()
    perf()
end

#main()
L039834(12) |> println

end # module

#=

Test Summary: | Pass  Total
Fibonacci     |    1      1

0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144,
0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181,

fmpz[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987]

[0.568864, 1.0, 0.920442, 1.0, 1.48931, 2.0, 2.40975, 3.0, 3.89906, 5.0]

  0.323110 seconds (13.20 k allocations: 206.234 KiB)
  0.007040 seconds (40.00 k allocations: 937.547 KiB)
  0.006942 seconds (49.49 k allocations: 1.137 MiB)

=#

### EOF ############################################################

#= COMMENTS ########################################################
# An alternative implementation mimicking Python's elegant way with yield:

using ResumableFunctions

@resumable function fibonnaci(n::Int) :: Int
  a, b = 0, 1
  for i ∈ 1:n-1
    @yield a
    a, b = b, a + b
  end
  a
end

for fib ∈ fibonnaci(10)
    println(fib)
end

=#

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
module CGP

using Logging
using PaddedViews
using Distributions

Logging.configure(level=INFO)

include("config.jl")
include("chromosome.jl")
include("distance.jl")
include("mutation.jl")
include("crossover.jl")
include("chromosomes/cgp.jl")
include("chromosomes/pcgp.jl")
include("logging.jl")
include("EAs/oneplus.jl")
include("EAs/cgpneat.jl")
include("EAs/ga.jl")
include("EAs/cmaes.jl")

EAs = [oneplus, GA]
chromosomes = [CGPChromo, PCGPChromo]
mutations = [:gene_mutate, :mixed_node_mutate, :mixed_subtree_mutate, :adaptive_node_mutate,
             :adaptive_subtree_mutate]
crossovers = [:single_point_crossover, :random_node_crossover, :aligned_node_crossover,
              :proportional_crossover, :output_graph_crossover, :subgraph_crossover]
distances = [:genetic_distance, :positional_distance, :constant_functional_distance,
             :random_functional_distance, :active_distance]

end


text:
using AbstractPlotting
using Meshing: SignedDistanceField, HyperRectangle, Vec, HomogenousMesh, MarchingTetrahedra

# algebraic dual boundary of the elliptopes
f = x -> (x[1]^2*x[2]^2 + x[1]^2*x[3]^2 + x[2]^2*x[3]^2 - 2*prod(x))

sdf = SignedDistanceField(f, HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2)),0.01)

m = HomogenousMesh(sdf, MarchingTetrahedra())

using MeshCat
vis = Visualizer()
open(vis)

setobject!(vis, m,
           MeshPhongMaterial(color=RGBA{Float32}(1, 0, 0, 0.5)))

# Makie
# mesh_scene = Scene()
# mesh!(mesh_scene, m, color=:blue)
# display(mesh_scene)


text:
"""
    Code for a top hat mask shape for use with CCF
Author: Eric Ford
Created: August 2020
"""

"""   TopHatCCFMask
The standard tophat mask with one parameter, it's half width as a velocity in m/s.
Mask weights are stored separately in a line list.
"""
struct TopHatCCFMask <: AbstractCCFMaskShape
    half_width::Float64
end

""" `TopHatCCFMask( full_width )`` """
function TopHatCCFMask(w::Real=default_v_width)
    @assert 0 < w < default_v_max
    TopHatCCFMask(w/2)
end

""" `TopHatCCFMask( inst  )` """
function TopHatCCFMask(inst::AbstractInstrument; scale_factor::Real = 1)
    w = scale_factor * default_ccf_mask_v_width(inst)
    TopHatCCFMask(w/2)
end

λ_min(m::TopHatCCFMask,λ::Real) = λ/calc_doppler_factor(m.half_width)
λ_max(m::TopHatCCFMask,λ::Real) = λ*calc_doppler_factor(m.half_width)

#global num_error_msgs = 0

function integrate(m::TopHatCCFMask, v_lo::Real,v_hi::Real)
#=    if ! (-m.half_width*1.01 <= v_lo <= m.half_width*1.01) || !(-m.half_width*1.01 <= v_hi <= m.half_width*1.01)
        global num_error_msgs += 1
        if num_error_msgs < 100
        println("# v_lo = ", v_lo, "  v_hi = ", v_hi, "   half_width = ", m.half_width)
#        @assert -m.half_width <= v_hi <= m.half_width
#        @assert -m.half_width <= v_lo <= m.half_width
        end
    end
    =#
    #quadgk(m, v_lo, v_hi)[1]
    0.5*(v_hi-v_lo)/m.half_width
    #0.5*(v_hi-v_lo)/m.half_width * (1+m.half_width/speed_of_light_mps)/(1+(m.half_width/speed_of_light_mps)^2/2)
end

""" Functor for returning constant for any Δv <= width.  """
function (m::TopHatCCFMask)(Δv::Real)
    # TODO: Change once implement generic CCF mask shapes (currently, I don't understand why the normalization differs for the two)
    #return abs2(Δv)<=abs2(m.half_width) ? 0.5/m.half_width : zero(Δv)  # Old for project_mask that assumes tophat
    return abs(Δv)<=abs(m.half_width) ? 0.5/m.half_width : zero(Δv)  # Old for project_mask that assumes tophat
    #return abs2(Δv)<=abs2(m.half_width) ? 1.0 : zero(Δv)    # New version for project_mask that works with a probability density
end

function mask_with_increased_fwhm(m::TopHatCCFMask, Δfwhm::Real )
    width = 2*m.half_width+Δfwhm
    return TopHatCCFMask(width)
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
# ==================================
# FOM
#
# system type: continuous LTI system
# state dimension: 1006
# input dimension: 1
# ==================================
using ReachabilityBenchmarks, MathematicalSystems, LazySets, MAT

function fom_model()
    file = matopen(@relpath "fom.mat")

    # system matrix
    A = float(read(file, "A"))  # the matrix has Int entries

    # input matrix
    B = read(file, "B")

    # state domain
    X = Universe(1006)

    # input domain
    U = BallInf([0.0], 1.0)

    # continuous LTI system
    S = @system(x' = Ax + Bu, x ∈ X, u ∈ U)

    return S
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
# Entry for script
function main_FVGQ_1_kalman(args=ARGS)
    d = parse_commandline_FVGQ_1_kalman(args)
    return estimate_FVGQ_1_kalman((; d...)) # to named tuple
end

function estimate_FVGQ_1_kalman(d)
    # Or move these into main package when loading?
    Turing.setadbackend(:zygote)
    HMCExamples.set_BLAS_threads()

    # load data relative to the current path
    data_path = joinpath(pkgdir(HMCExamples), d.data_path)
    z = collect(Matrix(DataFrame(CSV.File(data_path)))')

    # Create the perturbation and the turing models
    m = PerturbationModel(HMCExamples.FVGQ20)
    p_d = (β = d.beta, h = d.h, κ = d.kappa, χ = d.chi, γR = d.gamma_R, γΠ = d.gamma_Pi, Πbar = d.Pi_bar, ρd = d.rho_d, ρφ = d.rho_psi, ρg = d.rho_g, g_bar = d.g_bar, σ_A = d.sigma_A, σ_d = d.sigma_d, σ_φ = d.sigma_psi, σ_μ = d.sigma_mu, σ_m = d.sigma_m, σ_g = d.sigma_g, Λμ = d.Lambda_mu, ΛA = d.Lambda_A)
    p_f = (ϑ = 1.0, δ = d.delta, ε = d.epsilon, ϕ = d.phi, γ2 = d.gamma2, Ω_ii = d.Omega_ii, α = d.alpha, γy = d.gamma_y, θp = d.theta_p)
    c = SolverCache(m, Val(1), p_d)
    #create H prior
    Hx = zeros(6, m.n_x)
    Hy = zeros(6, m.n_y)
    Hy[1, 19] = 1 # Π
    Hy[2, 20] = 1 # R
    Hy[3, 18] = 1 # dw, trend is μz
    Hy[4, 18] = 1 # dy, trend is μz
    Hy[6, 24] = 1 # μ-1
    params = (β = Gamma_tr(d.beta_prior[1],d.beta_prior[2]),
    h = Beta_tr(d.h_prior[1], d.h_prior[2]),
    κ = (d.kappa_prior[1], d.kappa_prior[2], d.kappa_prior[3], d.kappa_prior[4]),
    γΠ = (d.gamma_Pi_prior[1], d.gamma_Pi_prior[2], d.gamma_Pi_prior[3], d.gamma_Pi_prior[4]),
    χ = Beta_tr(d.chi_prior[1], d.chi_prior[2]),
    γR = Beta_tr(d.gamma_R_prior[1], d.gamma_R_prior[2]),
    Πbar = Gamma_tr(d.Pi_bar_prior[1], d.Pi_bar_prior[2]),
    ρd = Beta_tr(d.rho_d_prior[1], d.rho_d_prior[2]),
    ρφ = Beta_tr(d.rho_psi_prior[1], d.rho_psi_prior[2]),
    ρg = Beta_tr(d.rho_g_prior[1], d.rho_g_prior[2]),
    g_bar = Beta_tr(d.g_bar_prior[1], d.g_bar_prior[2]),
    σ_A = InvGamma_tr(d.sigma_A_prior[1], d.sigma_A_prior[2]),
    σ_d = InvGamma_tr(d.sigma_d_prior[1], d.sigma_d_prior[2]),
    σ_φ = InvGamma_tr(d.sigma_psi_prior[1], d.sigma_psi_prior[2]),
    σ_μ = InvGamma_tr(d.sigma_mu_prior[1], d.sigma_mu_prior[2]),
    σ_m = InvGamma_tr(d.sigma_m_prior[1], d.sigma_m_prior[2]),
    σ_g = InvGamma_tr(d.sigma_g_prior[1], d.sigma_g_prior[2]),
    Λμ = Gamma_tr(d.Lambda_mu_prior[1], d.Lambda_mu_prior[2]),
    ΛA = Gamma_tr(d.Lambda_A_prior[1], d.Lambda_A_prior[2]),
    Hx = Hx,
    Hy = Hy)

    turing_model = FVGQ20_kalman(
        z, m, p_f, params, c, PerturbationSolverSettings(;ϵ_BK = d.epsilon_BK, print_level = d.print_level)
    )

    # Sampler
    name = "FQGV-kalman-s$(d.num_samples)-seed$(d.seed)"
    include_vars = ["β_draw", "h", "κ", "χ", "γR", "γΠ", "Πbar_draw", "ρd", "ρφ", "ρg", "g_bar", "σ_A", "σ_d", "σ_φ", "σ_μ", "σ_m", "σ_g", "Λμ", "ΛA"]  # variables to log
    callback = TensorBoardCallback(d.results_path; name, include=include_vars)
    num_adapts = convert(Int64, floor(d.num_samples * d.adapts_burnin_prop))

    Random.seed!(d.seed)
    @info "Generating $(d.num_samples) samples with $(num_adapts) adapts across $(d.num_chains) chains"

    chain = sample(
        turing_model,
        NUTS(num_adapts, d.target_acceptance_rate; d.max_depth),
        MCMCThreads(),
        d.num_samples,
        d.num_chains;
        init_params=[p_d...],
        progress=true,
        save_state=true,
        callback,
    )

    # Store parameters in log directory
    parameter_save_path = joinpath(callback.logger.logdir, "parameters.json")

    @info "Storing Parameters at $(parameter_save_path) "
    open(parameter_save_path, "w") do f
        write(f, JSON.json(d))
    end

    # Calculate and save results into the logdir
    calculate_experiment_results(chain, callback.logger, include_vars)
end

function parse_commandline_FVGQ_1_kalman(args)
    s = ArgParseSettings(; fromfile_prefix_chars=['@'])

    # See the appropriate _defaults.txt file for the default vvalues.
    @add_arg_table! s begin
        "--data_path"
        help = "relative path to data from the root of the package"
        arg_type = String
        "--beta"
        help = "Initialization of parameters"
        arg_type = Float64
        "--h"
        help = "Initialization of parameters"
        arg_type = Float64
        "--kappa"
        help = "Initialization of parameters"
        arg_type = Float64
        "--chi"
        help = "Initialization of parameters"
        arg_type = Float64
        "--gamma_R"
        help = "Initialization of parameters"
        arg_type = Float64
        "--gamma_Pi"
        help = "Initialization of parameters"
        arg_type = Float64
        "--Pi_bar"
        help = "Initialization of parameters"
        arg_type = Float64
        "--rho_d"
        help = "Initialization of parameters"
        arg_type = Float64
        "--rho_psi"
        help = "Initialization of parameters"
        arg_type = Float64
        "--rho_g"
        help = "Initialization of parameters"
        arg_type = Float64
        "--g_bar"
        help = "Initialization of parameters"
        arg_type = Float64
        "--sigma_A"
        help = "Initialization of parameters"
        arg_type = Float64
        "--sigma_d"
        help = "Initialization of parameters"
        arg_type = Float64
        "--sigma_psi"
        help = "Initialization of parameters"
        arg_type = Float64
        "--sigma_mu"
        help = "Initialization of parameters"
        arg_type = Float64
        "--sigma_m"
        help = "Initialization of parameters"
        arg_type = Float64
        "--sigma_g"
        help = "Initialization of parameters"
        arg_type = Float64
        "--Lambda_mu"
        help = "Initialization of parameters"
        arg_type = Float64
        "--Lambda_A"
        help = "Initialization of parameters"
        arg_type = Float64
        "--delta"
        help = "Value of fixed parameters"
        arg_type = Float64
        "--epsilon"
        help = "Value of fixed parameters"
        arg_type = Float64
        "--phi"
        help = "Value of fixed parameters"
        arg_type = Float64
        "--gamma2"
        help = "Value of fixed parameters"
        arg_type = Float64
        "--Omega_ii"
        help = "Value of fixed parameters"
        arg_type = Float64
        "--alpha"
        help = "Value of fixed parameters"
        arg_type = Float64
        "--gamma_y"
        help = "Value of fixed parameters"
        arg_type = Float64
        "--theta_p"
        help = "Value of fixed parameters"
        arg_type = Float64

        "--kappa_prior"
        help = "Value of fixed parameters"
        arg_type = Vector{Float64}
        "--gamma_Pi_prior"
        help = "Value of fixed parameters"
        arg_type = Vector{Float64}

        "--beta_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--h_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--chi_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--gamma_R_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--Pi_bar_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--rho_d_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--rho_psi_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--rho_g_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--g_bar_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}

        "--sigma_A_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--sigma_d_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--sigma_psi_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--sigma_mu_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--sigma_m_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--sigma_g_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}

        "--Lambda_mu_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--Lambda_A_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}

        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--num_samples"
        help = "samples to draw in chain"
        arg_type = Int64
        "--num_chains"
        help = "number of chains to sample"
        arg_type = Int64
        "--adapts_burnin_prop"
        help = "Proportion of Adaptations burned in"
        arg_type = Float64
        "--target_acceptance_rate"
        help = "Target acceptance rate for dual averaging for NUTS"
        arg_type = Float64
        "--max_depth"
        help = "Maximum depth for NUTS"
        arg_type = Int64
        "--seed"
        help = "Random number seed"
        arg_type = Int64
        "--results_path"
        arg_type = String
        help = "Location to store results and logs"
        "--print_level"
        arg_type = Int64
        help = "Print level for output during sampling"
        "--epsilon_BK"
        arg_type = Float64
        help = "Threshold for Checking Blanchard-Khan condition"
        "--use_solution_cache"
        arg_type = Bool
        help = "Use solution cache in perturbation solutions"

    end

    args_with_default = vcat("@$(pkgdir(HMCExamples))/src/FVGQ_1_kalman_defaults.txt", args)
    return parse_args(args_with_default, s; as_symbols=true)

end

max_stars_repo_name: HighDimensionalEconLab/HMCExamples

max_issues_repo_name: HighDimensionalEconLab/HMCExamples

max_forks_repo_name: HighDimensionalEconLab/HMCExamples


text:
using Snopt
using DelimitedFiles
using PyPlot
import ForwardDiff

# set up boundary constraint wrapper function
function boundary_wrapper(x, params)
    # include relevant globals
    params.boundary_center
    params.boundary_radius

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # get and return boundary distances
    return ff.circle_boundary(boundary_center, boundary_radius, turbine_x, turbine_y)
end

# set up spacing constraint wrapper function
function spacing_wrapper(x, params)
    # include relevant globals
    params.rotor_diameter

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # get and return spacing distances
    return 2.0*rotor_diameter[1] .- ff.turbine_spacing(turbine_x,turbine_y)
end

# set up objective wrapper function
@everywhere function aep_wrapper(x, params)

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # calculate AEP
    AEP = obj_scale*ff.calculate_aep(turbine_x, turbine_y, params.turbine_z, params.rotor_diameter,
    params.hub_height, params.turbine_yaw, params.ct_models, params.generator_efficiency, params.cut_in_speed,
    params.cut_out_speed, params.rated_speed, params.rated_power, params.windresource, params.power_models, params.model_set,
                rotor_sample_points_y=params.rotor_points_y,rotor_sample_points_z=params.rotor_points_z)

    # return the objective as an array
    return [AEP]
end

# set up optimization problem wrapper function
function wind_farm_opt(x)

    # calculate spacing constraint value and jacobian
    spacing_con = spacing_wrapper(x)
    ds_dx = ForwardDiff.jacobian(spacing_wrapper, x)

    # calculate boundary constraint and jacobian
    boundary_con = boundary_wrapper(x)
    db_dx = ForwardDiff.jacobian(boundary_wrapper, x)

    # combine constaint values and jacobians into overall constaint value and jacobian arrays
    c = [spacing_con; boundary_con]
    dcdx = [ds_dx; db_dx]

    # calculate the objective function and jacobian (negative sign in order to maximize AEP)
    AEP = -aep_wrapper(x)[1]
    dAEP_dx = -ForwardDiff.jacobian(aep_wrapper,x)

    # set fail flag to false
    fail = false

    # return objective, constraint, and jacobian values
    return AEP, c, dAEP_dx, dcdx, fail
end

# set globals for use in wrapper functions
struct params_struct2{MS, AF, F, I, ACTM, WR, APM}
    model_set::MS
    rotor_points_y::AF
    rotor_points_z::AF
    turbine_z::AF
    rotor_diameter::AF
    boundary_center::AF
    boundary_radius::F
    obj_scale::I
    hub_height::AF
    turbine_yaw::AF
    ct_models::ACTM
    generator_efficiency::AF
    cut_in_speed::AF
    cut_out_speed::AF
    rated_speed::AF
    rated_power::AF
    windresource::WR
    power_models::APM
end

# import model set with wind farm and related details
include("../model_sets/model_set_6.jl")

# scale objective to be between 0 and 1
obj_scale = 1E-11

# set wind farm boundary parameters
boundary_center = [0.0,0.0]
boundary_radius = 300.0

# initialize struct for opt params
params = params_struct2(model_set, rotor_points_y, rotor_points_z, turbine_z,
    rotor_diameter, boundary_center, boundary_radius, obj_scale, hub_height, turbine_yaw,
    ct_models, generator_efficiency, cut_in_speed, cut_out_speed, rated_speed, rated_power,
    windresource, power_models)

# initialize design variable array
x = [copy(turbine_x);copy(turbine_y)]
xinit = deepcopy(x)
# report initial objective value
println("starting objective value: ", aep_wrapper(x, params)[1])

# add initial turbine location to plot
for i = 1:length(turbine_x)
    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[1]/2.0, fill=false,color="C0"))
end

# set general lower and upper bounds for design variables
lb = zeros(length(x)) .- boundary_radius
ub = zeros(length(x)) .+ boundary_radius

# set up options for SNOPT
options = Dict{String, Any}()
options["Derivative option"] = 1
options["Verify level"] = 3
options["Major optimality tolerance"] = 1e-5
options["Major iteration limit"] = 1e6
options["Summary file"] = "./snopt-opt2-summary.out"
options["Print file"] = "./snopt-opt2-print.out"

# generate wrapper function surrogates
spacing_wrapper(x) = spacing_wrapper(x, params)
aep_wrapper(x) = aep_wrapper(x, params)
boundary_wrapper(x) = boundary_wrapper(x, params)
obj_func(x) = wind_farm_opt(x)

# run and time optimization
t1 = time()
xopt, fopt, info = snopt(obj_func, x, lb, ub, options)
t2 = time()
clkt = t2-t1

# print optimization results
println("Finished in : ", clkt, " (s)")
println("info: ", info)
println("end objective value: ", aep_wrapper(xopt))

# extract final turbine locations
turbine_x = copy(xopt[1:nturbines])
turbine_y = copy(xopt[nturbines+1:end])

# add final turbine locations to plot
for i = 1:length(turbine_x)
    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[1]/2.0, fill=false,color="C1", linestyle="--"))
end

# add wind farm boundary to plot
plt.gcf().gca().add_artist(plt.Circle((boundary_center[1],boundary_center[2]), boundary_radius, fill=false,color="C2"))

# set up and show plot
axis("square")
xlim(-boundary_radius-200,boundary_radius+200)
ylim(-boundary_radius-200,boundary_radius+200)
plt.show()


text:
using Sockets

users = eachline("./smtp-users.txt")

function smtp_timeout()
    sleep(1.5)
    return nothing
end

function connect_smtp(target, port)
    try
        c = Channel(2)
        @async put!(c, connect(target, port))
        @async put!(c, smtp_timeout())
        socket = take!(c)
        if socket !== nothing
            # Read until carriage return
            connection_banner = String(readuntil(socket, 0x0d))
            # SMTP Connected
            if occursin("220", connection_banner)
                return socket
            end
            return nothing
        end
        return nothing
    catch e
        return nothing
    end
end

function vrfy(socket, name)
    write(socket, "VRFY $name \r\n")
    response = String(readuntil(socket, 0x0d))
    return occursin("252", response) && occursin(name, response)
end

function scan_smtp_target(target, port)
    connection = connect_smtp(target, port)
    if connection !== nothing
        print("Users found on $target SMTP server on port $port:\n")
        for user in users
            try
                if vrfy(connection, user)
                    println(user)
                end
            catch e
                if isa(e, Base.IOError)
                    # Connection reset by peer, retry...
                    if occursin("ECONNRESET", e.msg)
                        connection = connect_smtp(target, port)
                        if vrfy(connection, user)
                            println(user)
                        end
                    end
                end
            end
        end
        close(connection)
        return true
    end
    return false
end

function smtp_scan(target, openPorts)
    got_smtp_server = false
    if !got_smtp_server && openPorts[25] !== nothing
        got_smtp_server = scan_smtp_target(target, 25)
    end
    if !got_smtp_server && openPorts[465] !== nothing
        got_smtp_server = scan_smtp_target(target, 465)
    end
    if !got_smtp_server && openPorts[587] !== nothing
        got_smtp_server = scan_smtp_target(target, 587)
    end
    if !got_smtp_server
        print("\nNo SMTP server, skipping VRFY scan...\n")
    end
end

max_stars_repo_name: thecrypters/nis

max_issues_repo_name: thecrypters/nis

max_forks_repo_name: thecrypters/julia-host-scanner


text:
"""
        summary(decisionMat, weights, fns, methods)

Apply more methods for a given decision problem. The methods accept standart number of arguments.

# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n candidate (or strategy) and m criteria
 - `weights::Array{Float64, 1}`: m-vector of weights that sum up to 1.0. If the sum of weights is not 1.0, it is automatically normalized.
 - `fns::Array{Function, 1}`: m-vector of function that are either minimize or maximize.
 - `methods::Array{Symbol, 1}`: Array of symbols. The elements can be :topsis, :electre, :cocoso, :copras, :moora, :vikor, :grey, :aras, :saw, :wpm, :waspas, :edas, :marcos, :mabac, :mairca, :copras, :critic

# Description
    This method outputs a summarized output using more than MCDM methods in a comparable way.

    # Output
- `::DataFrame`: A DataFrame object, methods in columns, and alternatives in rows. Green check symbol indicates the selected alternative as the best by the corresponding method.

# Examples
```julia-repl
julia> df = DataFrame(
:age        => [6.0, 4, 12],
:size       => [140.0, 90, 140],
:price      => [150000.0, 100000, 75000],
:distance   => [950.0, 1500, 550],
:population => [1500.0, 2000, 1100]);

julia> methods = [:topsis, :electre, :vikor, :moora, :cocoso, :wpm, :waspas]

julia> w  = [0.036, 0.192, 0.326, 0.326, 0.12];

julia> fns = [maximum, minimum, maximum, maximum, maximum];

julia> result = summary(df, w, fns, methods)
3×7 DataFrame
 Row │ topsis  electre  cocoso  moora   vikor   wpm     waspas
     │ String  String   String  String  String  String  String
─────┼─────────────────────────────────────────────────────────
   1 │                           ✅      ✅
   2 │  ✅      ✅       ✅                      ✅      ✅
   3 │

```
"""
function summary(
    decisionMat::DataFrame,
    weights::Array{Float64,1},
    fns::Array{Function,1},
    methods::Array{Symbol,1})

    nmethods = length(methods)
    nalternatives, _ = size(decisionMat)

    resultdf = DataFrame()

    check = " ✅ "

    if :topsis in methods
        topresult = topsis(decisionMat, weights, fns)
        resultdf[:,:topsis] = map(x -> if topresult.bestIndex == x check else " " end, 1:nalternatives)
    end

    if :electre in methods
        # electre best index is a tuple
        result = electre(decisionMat, weights, fns)
        resultdf[:,:electre] = map(x -> if x in result.bestIndex check else " " end, 1:nalternatives)
    end

    if :cocoso in methods
        result = cocoso(decisionMat, weights, fns)
        resultdf[:,:cocoso] = map(x -> if result.bestIndex == x check else " " end, 1:nalternatives)
    end

    if :copras in methods
        result = copras(decisionMat, weights, fns)
        resultdf[:,:copras] = map(x -> if result.bestIndex == x check else " " end, 1:nalternatives)
    end

    if :moora in methods
        result = moora(decisionMat, weights, fns)
        resultdf[:,:moora] = map(x -> if result.bestIndex == x check else " " end, 1:nalternatives)
    end

    if :vikor in methods
        result = vikor(decisionMat, weights, fns)
        resultdf[:,:vikor] = map(x -> if result.bestIndex == x check else " " end, 1:nalternatives)
    end

    if :grey in methods
        result = grey(decisionMat, weights, fns)
        resultdf[:,:grey] = map(x -> if result.bestIndex == x check else " " end, 1:nalternatives)
    end

    if :aras in methods
        result = aras(decisionMat, weights, fns)
        resultdf[:,:aras] = map(x -> if result.bestIndex == x check else " " end, 1:nalternatives)
    end

    if :saw in methods
        result = saw(decisionMat, weights, fns)
        resultdf[:,:saw] = map(x -> if result.bestIndex == x check else " " end, 1:nalternatives)
    end

    if :wpm in methods
        result = wpm(decisionMat, weights, fns)
        resultdf[:,:wpm] = map(x -> if result.bestIndex == x check else " " end, 1:nalternatives)
    end

    if :waspas in methods
        result = waspas(decisionMat, weights, fns)
        resultdf[:,:waspas] = map(x -> if result.bestIndex == x check else " " end, 1:nalternatives)
    end

    if :edas in methods
        result = edas(decisionMat, weights, fns)
        resultdf[:,:edas] = map(x -> if result.bestIndex == x check else " " end, 1:nalternatives)
    end

    if :marcos in methods
        result = marcos(decisionMat, weights, fns)
        resultdf[:,:marcos] = map(x -> if result.bestIndex == x check else " " end, 1:nalternatives)
    end

    if :mabac in methods
        result = mabac(decisionMat, weights, fns)
        resultdf[:,:mabac] = map(x -> if result.bestIndex == x check else " " end, 1:nalternatives)
    end

    if :mairca in methods
        result = mairca(decisionMat, weights, fns)
        resultdf[:,:mairca] = map(x -> if result.bestIndex == x check else " " end, 1:nalternatives)
    end

    if :copras in methods
        result = copras(decisionMat, weights, fns)
        resultdf[:,:copras] = map(x -> if result.bestIndex == x check else " " end, 1:nalternatives)
    end

    if :critic in methods
        result = critic(decisionMat, fns)
        resultdf[:,:critic] = map(x -> if result.bestIndex == x check else " " end, 1:nalternatives)
    end

    return resultdf

end

"""
        summary(setting, methods)

Apply more methods for a given decision problem. The methods accept standart number of arguments.

# Arguments:
 - `setting::MCDMSetting`: MCDMSetting object.
 - `methods::Array{Symbol, 1}`: Array of symbols. The elements can be :topsis, :electre, :cocoso, :copras, :moora, :vikor, :grey, :aras, :saw, :wpm, :waspas, :edas, :marcos, :mabac, :mairca, :copras, :critic

# Description
    This method outputs a summarized output using more than MCDM methods in a comparable way.
"""
function summary(setting::MCDMSetting, methods::Array{Symbol,1})
    summary(
        setting.df,
        setting.weights,
        setting.fns,
        methods
    )
end

max_stars_repo_name: drvinceknight/JMcDM

max_issues_repo_name: drvinceknight/JMcDM

max_forks_repo_name: drvinceknight/JMcDM


text:
import .StatsFuns
using .StatsFuns: xlogx, xlogy, logistic, logit, log1psq, log1pexp,
    logsumexp, logaddexp, logsubexp
using Base.Broadcast: broadcasted

@adjoint function xlogx(x::Real)
    result, dx = ∇xlogx(x)
    back(δ) = (dx * δ,)
    return result, back
end
@adjoint function broadcasted(::typeof(xlogx), x::Numeric)
    result, dx = ∇xlogx(x)
    back(δ) = (nothing, unbroadcast(x, δ .* dx))
    return result, back
end
function ∇xlogx(x::Numeric)
    logx = log.(x)
    xlogx = x .* logx
    result = @. ifelse(iszero(x), zero(xlogx), xlogx)
    dx = @. one(x) + logx
    return result, dx
end

@adjoint function logistic(x::Real)
    y = logistic(x)
    return y, Δ->(Δ * y * (1 - y),)
end

@adjoint logit(x::Real) = logit(x), Δ->(Δ / (x * (1 - x)),)

@adjoint log1psq(x::Real) = log1psq(x), Δ->(Δ * 2x / (1 + abs2(x)),)

@adjoint function log1pexp(x::Real)
    dx = ∂log1pexp(x)
    return log1pexp(x), δ -> (δ * dx,)
end
@adjoint function broadcasted(::typeof(log1pexp), x::Numeric)
    dx = ∂log1pexp.(x)
    return log1pexp.(x), δ -> (nothing, unbroadcast(x, δ .* dx))
end
∂log1pexp(x::Real)    = x < 18.0 ? logistic(x) : x < 33.3 ? one(x) - exp(-x) : oftype(exp(x), 1)
∂log1pexp(x::Float32) = x < 9f0  ? logistic(x) : x < 16f0 ? one(x) - exp(-x) : oftype(exp(x), 1)

@adjoint function logsumexp(X::AbstractArray{<:Real}; dims=:)
    lse = logsumexp(X; dims=dims)
    return lse, Δ -> (Δ .* exp.(X .- lse),)
end

@adjoint function xlogy(x::Real, y::Real)
    result, dx, dy = ∇xlogy(x, y)
    back(δ) = (δ * dx, δ * dy)
    return result, back
end
@adjoint function broadcasted(::typeof(xlogy), x::Numeric, y::Numeric)
    result, dx, dy = ∇xlogy(x, y)
    back(δ) = (nothing, unbroadcast(x, δ .* dx), unbroadcast(y, δ .* dy))
    return result, back
end
function ∇xlogy(x::Numeric, y::Numeric)
    dx = logy = log.(y)
    dy = x ./ y
    xlogy = x .* logy
    result = @. ifelse(iszero(x) & !isnan(y), zero(xlogy), xlogy)
    return result, dx, dy
end

@adjoint function logaddexp(x::Real, y::Real)
    result, dx, dy = ∇logaddexp(x, y)
    back(δ) = (δ * dx, δ * dy)
    return result, back
end
@adjoint function broadcasted(::typeof(logaddexp), x::Numeric, y::Numeric)
    result, dx, dy = ∇logaddexp(x, y)
    back(δ) = (nothing, unbroadcast(x, δ .* dx), unbroadcast(y, δ .* dy))
    return result, back
end
function ∇logaddexp(x::Numeric, y::Numeric)
    result = logaddexp.(x, y)
    t = @. exp(-abs(x - y))
    dx, dy = select(x .≥ y, inv.(one.(t) .+ t), t ./ (one.(t) .+ t))
    return result, dx, dy
end

@adjoint function logsubexp(x::Real, y::Real)
    result, dx, dy = ∇logsubexp(x, y)
    back(δ) = (δ * dx, δ * dy)
    return result, back
end
@adjoint function broadcasted(::typeof(logsubexp), x::Numeric, y::Numeric)
    result, dx, dy = ∇logsubexp(x, y)
    back(δ) = (nothing, unbroadcast(x, δ .* dx), unbroadcast(y, δ .* dy))
    return result, back
end
function ∇logsubexp(x::Numeric, y::Numeric)
    result = logsubexp.(x, y)
    t = @. -inv(expm1(-abs(x - y)))
    dx, dy = select(x .≥ y, t, one.(t) .- t)
    return result, dx, dy
end

"""
	select(cond, x, y)

The call `a, b = select(cond, x, y)` constructs two arrays `a, b`, where
`a[i], b[i] = x[i], y[i]` if `cond[i]` is `true`, and `a[i], b[i] = y[i], x[i]`
 if `cond[i]` is `false`.
"""
select(cond, x, y) = ifelse.(cond, x, y), ifelse.(cond, y, x)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
__precompile__()
module RandomObjects

export random

@inline random(::Type{String}; l::Integer=rand(2:10)) = randstring(l)
@inline random(T, n::Int) = T[random(T) for _ in 1:n]
random(T, dims::Dims) = reshape(random(T, prod(dims)), dims)
random(T, dims::Int...) = random(T, dims)

@inline random{T}(::Type{Vector{T}}; l::Integer=rand(0:20)) = random(T,l)
@inline random{T <: Number}(::Type{T}) = rand(T)
random{K,V}(D::Type{Dict{K,V}}) = D(random(Vector{Tuple{K, V}}))

function random{T, n}(::Type{Array{T, n}})
    dims = tuple(rand(1:10, n)...) :: NTuple{n, Int}
    random(T, dims)
end

function random_from_fields_impl{T}(::Type{T})
    fieldtypes = T.types
    random_types = [:(random($S)) for S in fieldtypes]
    Expr(:call, T, random_types...)
end

function random_tuple_impl{T <: Tuple}(::Type{T})
    fieldtypes = T.parameters
    random_types = [:(random($S)) for S in fieldtypes]
    Expr(:tuple, random_types...)
end

@generated random{T}(::Type{T}) = random_from_fields_impl(T)
@generated random{T <: Tuple}(::Type{T}) = random_tuple_impl(T)

end # module


text:
module FirstPkg

using Distributions

export pdf_normal

"""
   pdf_normal(x, mu, sigma)

Return pdf of a normal distribution with mean mu and standard deviation sigma at point x.
"""
function pdf_normal(x::Number, mu::Number, sigma::Number)
  dist = Normal(mu,sigma)
  pdf(dist,x)
end

end # module


text:
@polly function kernel_gesummv(alpha, beta, A, B, tmp, x, y)
    n = size(tmp,1)
    for i = 1:n
        tmp[i] = zero(eltype(tmp))
        y[i] = zero(eltype(tmp))
        for j = 1:n
            tmp[i] = A[i,j] * x[j] + tmp[i]
            y[i] = B[i,j] * x[j] + y[i]
        end
        y[i] = alpha * tmp[i] + beta * y[i]
    end
end

let
    n = 1300

    A = zeros(Float32, n, n)
    B = zeros(Float32, n, n)
    tmp = zeros(Float32, n)
    x = zeros(Float32, n)
    y = zeros(Float32, n)

    alpha = Float32(1.5)
    beta = Float32(1.2)
    for i = 1:n
        x[i] = (i % n) / n
        for j = 1:n
            A[i,j] = ((i*j+1) % n) / n
            B[i,j] = ((i*j+2) % n) / n
        end
    end

    SUITE["gesummv"] = @benchmarkable kernel_gesummv($alpha, $beta, A, B, tmp, x, y) setup = (A = copy($A); B = copy($B); tmp = copy($tmp); x = copy($x); y = copy($y))
end


text:
__precompile__()

module Intan

using HDF5, DistributedArrays, SpikeSorting, Gtk.ShortNames, Cairo, MAT,JLD, DSP

using SharedArrays, Libdl, Distributed, Dates

import SpikeSorting.Vec2, SpikeSorting.RubberBand, SpikeSorting.rb_set, SpikeSorting.rb_draw, SpikeSorting.identity_matrix, SpikeSorting.select_color

if VERSION > v"0.7-"
    const Void = Nothing
    const is_linux() = Sys.islinux()
    const is_apple() = Sys.isapple()
    const is_windows() = Sys.iswindows()
    const is_unix() = Sys.isunix()
    const setproperty! = set_gtk_property!
    const getproperty = get_gtk_property
    const linspace(x,y,z) = range(x,stop=y,length=z)
    const indmin = argmin
    const find = findall
    const indmax = argmax
    Base.round(x::Number, digits) = round(x; digits=digits)
end

export makeRHD, RHD2164, RHD2132,  makegui, Debug, Intan_GUI

#const MyFilter = DSP.Filters.DF2TFilter{DSP.Filters.SecondOrderSections{Float64,Float64},Array{Float64,2}}
const MyFilter = DSP.DF2TFilter{SecondOrderSections{:z, Float64, Float64}, Matrix{Float64}}

include("types.jl")
include("constants.jl")
include("gtk_helpers.jl")
include("gui/single_channel_canvas.jl")
include("gui/cluster_treeview.jl")
include("gui/multi_canvas_selection.jl")
include("gui/soft_scope.jl")
include("gui/event_display.jl")
include("gui/spectrogram.jl")
include("gui/multi_channel_display.jl")
include("gui/reference_channel.jl")
include("gui/parameter_table.jl")
include("gui/filters.jl")
include("gui/thresholds.jl")
include("gui/save_pref.jl")
include("gui/clock.jl")
include("gui/digital_out_gui.jl")

#functions needed for evaluation board from Intan
include("ICs/Opal_Kelly_FPGA/opalkelly_gui.jl")
include("ICs/Opal_Kelly_FPGA/c_wrapper.jl")
include("ICs/Opal_Kelly_FPGA/dacs.jl")
include("ICs/Opal_Kelly_FPGA/ttls.jl")
include("ICs/Opal_Kelly_FPGA/intan_expansion_board.jl")
include("ICs/Opal_Kelly_FPGA/cable_delay.jl")
include("ICs/Opal_Kelly_FPGA/headstage_cmd.jl")

include("rhd2000evalboard.jl")
include("rhd2000registers.jl")
include("gui.jl")
include("digital_output.jl")
include("plotting.jl")
include("offline.jl")
include("task_list.jl")
include("save_load.jl")
include("benchmark.jl")

#Spike Sorting
include("sorting/windows.jl")
include("sorting/templates.jl")

#Graphics
include("graphics_help/cairo.jl")

end