text:
abstract type AbstractLattice{E,L,T} end

#######################################################################
# Sublat (sublattice)
#######################################################################
struct Sublat{E,T}
    sites::Vector{SVector{E,T}}
    name::NameType
end

Base.empty(s::Sublat) = Sublat(empty(s.sites), s.name)

Base.copy(s::Sublat) = Sublat(copy(s.sites), s.name)

Base.show(io::IO, s::Sublat{E,T}) where {E,T} = print(io,
"Sublat{$E,$T} : sublattice of $T-typed sites in $(E)D space
  Sites    : $(length(s.sites))
  Name     : $(displayname(s))")

displayname(s::Sublat) = s.name == nametype(:_) ? "pending" : string(":", s.name)
# displayorbitals(s::Sublat) = string("(", join(string.(":", s.orbitals), ", "), ")")
nsites(s::Sublat) = length(s.sites)

# External API #

"""
    sublat(sites...; name::$(NameType), orbitals = (:noname,))
    sublat(sites::Vector{<:SVector}; name::$(NameType), orbitals = (:noname,))

Create a `Sublat{E,T,D}` that adds a sublattice, of name `name`, with sites at positions
`sites` in `E` dimensional space, each of which hosts `D` different orbitals, with orbital
names specified by `orbitals`. Sites can be entered as tuples or `SVectors`.

# Examples
```jldoctest
julia> sublat((0.0, 0), (1, 1), (1, -1), name = :A, orbitals = (:upspin, :downspin))
Sublat{2,Float64,2} : sublattice in 2D space with 2 orbitals per site
  Sites    : 3
  Name     : :A
  Orbitals : (:upspin, :downspin)
```
"""
sublat(sites::Vector{<:SVector}; name = :_, kw...) =
    Sublat(sites, nametype(name))
sublat(vs::Union{Tuple,AbstractVector{<:Number}}...; kw...) = sublat(toSVectors(vs...); kw...)

# transform!(s::S, f::F) where {S <: Sublat,F <: Function} = (s.sites .= f.(s.sites); s)

#######################################################################
# Bravais
#######################################################################
struct Bravais{E,L,T,EL}
    matrix::SMatrix{E,L,T,EL}
    semibounded::SVector{L,Bool}
end

Bravais{E,T}() where {E,T} = Bravais(SMatrix{E,0,T,0}(), SVector{0,Bool}())

displayvectors(br::Bravais) = displayvectors(br.matrix)

Base.show(io::IO, b::Bravais{E,L,T}) where {E,L,T} = print(io,
"Bravais{$E,$L,$T} : set of $L Bravais vectors in $(E)D space.
  Vectors     : $(displayvectors(b))
  Matrix      : $(b.matrix),
  Semibounded : $(issemibounded(b) ? display_as_tuple(findall(b.semibounded)) : "none")")

# External API #

"""
    bravais(vecs...; semibounded = false)
    bravais(matrix; semibounded = false)

Create a `Bravais{E,L}` that adds `L` Bravais vectors `vecs` in `E` dimensional space,
alternatively given as the columns of matrix `mat`. For higher instantiation efficiency
enter `vecs` as `Tuple`s or `SVector`s and `mat` as `SMatrix`.

To create semibounded lattices along some or all Bravais vectors, use `semibounded`. A
`semibounded = true` makes all axes semibounded. A `semibounded = (axes::Int...)` indicates
the indices of axes to be made semibounded. A tuple of Booleans, one per Bravais axes,
`semibounded = issemi::NTuple{L,Bool}` is also allowed. Note that semibounded lattices
always extend toward positive multiples of Bravais vectors. To invert the direction, invert
the vectors.

We can scale a `b::Bravais` simply by multiplying it with a factor `a`, like `a * b`.

    bravais(lat::Lattice)
    bravais(h::Hamiltonian)

Obtain the Bravais matrix of lattice `lat` or Hamiltonian `h`

# Examples
```jldoctest
julia> bravais((1.0, 2), (3, 4))
Bravais{2,2,Float64} : set of 2 Bravais vectors in 2D space.
  Vectors : ((1.0, 2.0), (3.0, 4.0))
  Matrix  : [1.0 3.0; 2.0 4.0]
```

# See also:
    semibounded
"""
bravais(vs::Union{Tuple,AbstractVector,AbstractMatrix}...; semibounded = false, kw...) =
    (s = toSMatrix(vs...); Bravais(s, sanitize_semibounded(semibounded, s)))

bravais(lat::AbstractLattice) = lat.bravais.matrix

sanitize_semibounded(sb::Bool, ::SMatrix{E,L}) where {E,L} =
    SVector{L,Bool}(filltuple(sb, Val(L))) # need to specify type becaus of L=0 case
sanitize_semibounded(sb::Int, s::SMatrix{E,L}) where {E,L} =
    sanitize_semibounded((sb,), s)
sanitize_semibounded(sb::NTuple{L,Bool}, ::SMatrix{E,L}) where {E,L} = SVector{L,Bool}(sb)
sanitize_semibounded(sb, ::SMatrix{E,L}) where {E,L} =
    SVector{L,Bool}(ntuple(i -> i in sb, Val(L)))

transform(b::Bravais{E,0}, f::F) where {E,F<:Function} = b

function transform(b::Bravais{E,L,T}, f::F) where {E,L,T,F<:Function}
    svecs = let z = zero(SVector{E,T})
        ntuple(i -> f(b.matrix[:, i]) - f(z), Val(L))
    end
    matrix = hcat(svecs...)
    return Bravais(matrix, b.semibounded)
end

Base.:*(factor::Number, b::Bravais) = Bravais(factor * b.matrix, b.semibounded)
Base.:*(b::Bravais, factor::Number) = Bravais(b.matrix * factor, b.semibounded)

issemibounded(b::Bravais) = !iszero(b.semibounded)

#######################################################################
# Unitcell
#######################################################################
struct Unitcell{E,T,N}
    sites::Vector{SVector{E,T}}
    names::NTuple{N,NameType}
    offsets::Vector{Int}        # Linear site number offsets for each sublat
end                             # so that diff(offset) == sublatsites

Unitcell(sublats::Sublat...; kw...) = Unitcell(promote(sublats...); kw...)
Unitcell(sublats::NTuple{N,Sublat{E,T}};
    dim = Val(E), type = float(T), names = [s.name for s in sublats], kw...) where {N,E,T} =
    _Unitcell(sublats, dim, type, names)

# Dynamic dispatch
_Unitcell(sublats, dim::Integer, type, names) = _Unitcell(sublats, Val(dim), type, names)

function _Unitcell(sublats::NTuple{N,Sublat}, dim::Val{E}, type::Type{T}, names) where {N,E,T}
    sites = SVector{E,T}[]
    offsets = [0]  # length(offsets) == length(sublats) + 1
    for s in eachindex(sublats)
        for site in sublats[s].sites
            push!(sites, padright(site, Val(E)))
        end
        push!(offsets, length(sites))
    end
    names´ = uniquenames(sanitize_names(names, Val(N)))
    return Unitcell(sites, names´, offsets)
end

sanitize_names(name::Union{NameType,Int}, ::Val{N}) where {N} = ntuple(_ -> NameType(name), Val(N))
sanitize_names(names::AbstractVector, ::Val{N}) where {N} = ntuple(i -> NameType(names[i]), Val(N))
sanitize_names(names::NTuple{N,Union{NameType,Int}}, ::Val{N}) where {N} = NameType.(names)

function uniquenames(names::NTuple{N,NameType}) where {N}
    namesvec = [names...]
    allnames = NameType[:_]
    for i in eachindex(names)
        namesvec[i] in allnames && (namesvec[i] = uniquename(allnames, namesvec[i], i))
        push!(allnames, namesvec[i])
    end
    names´ = ntuple(i -> namesvec[i], Val(N))
    return names´
end

function uniquename(allnames, name, i)
    newname = nametype(Char(64+i)) # Lexicographic, starting from Char(65) = 'A'
    return newname in allnames ? uniquename(allnames, name, i + 1) : newname
end

sites(u::Unitcell) = u.sites
sites(u::Unitcell, s::Int) = view(u.sites, siterange(u, s))

siteindex(u::Unitcell, sublat, idx) = idx + u.offsets[sublat]

siterange(u::Unitcell, sublat) = (1+u.offsets[sublat]):u.offsets[sublat+1]

nsites(u::Unitcell) = length(u.sites)
nsites(u::Unitcell, sublat) = sublatsites(u)[sublat]

offsets(u::Unitcell) = u.offsets

sublat(u::Unitcell, siteidx) = Int(findlast(o -> o < siteidx, u.offsets))

sublatsites(u::Unitcell) = diff(u.offsets)

nsublats(u::Unitcell) = length(u.names)

sublats(u::Unitcell) = 1:nsublats(u)

transform!(u::Unitcell, f::Function) = (u.sites .= f.(u.sites); u)

Base.copy(u::Unitcell) = Unitcell(copy(u.sites), u.names, copy(u.offsets))

#######################################################################
# Lattice
#######################################################################
struct Lattice{E,L,T<:AbstractFloat,B<:Bravais{E,L,T},U<:Unitcell{E,T}} <: AbstractLattice{E,L,T}
    bravais::B
    unitcell::U
end

displaynames(l::AbstractLattice) = display_as_tuple(l.unitcell.names, ":")

function Base.show(io::IO, lat::Lattice)
    i = get(io, :indent, "")
    hassemi = issemibounded(lat.bravais)
    print(io, i, summary(lat), "\n",
"$i  Bravais vectors : $(displayvectors(lat.bravais.matrix; digits = 6))", hassemi ? "
$i    Semibounded   : $(display_as_tuple(findall(lat.bravais.semibounded)))" : "", "
$i  Sublattices     : $(nsublats(lat))
$i    Names         : $(displaynames(lat))
$i    Sites         : $(display_as_tuple(sublatsites(lat))) --> $(nsites(lat)) total per unit cell")
end

Base.summary(::Lattice{E,L,T}) where {E,L,T} =
    "Lattice{$E,$L,$T} : $(L)D lattice in $(E)D space"

# External API #

"""
    lattice([bravais::Bravais,] sublats::Sublat...; dim::Val{E}, type::T, names)

Create a `Lattice{E,L,T}` with Bravais matrix `bravais` and sublattices `sublats`
converted to a common  `E`-dimensional embedding space and type `T`. To override the
embedding  dimension `E`, use keyword `dim = Val(E)`. Similarly, override type `T` with
`type = T`.

The keywords `names` can be used to rename `sublats`. Given names can be replaced to ensure
that all sublattice names are unique.

See also `LatticePresets` for built-in lattices.

# Examples
```jldoctest
julia> lattice(bravais((1, 0)), sublat((0, 0)), sublat((0, Float32(1))); dim = Val(3))
Lattice{3,1,Float32} : 1D lattice in 3D space
  Bravais vectors : ((1.0, 0.0, 0.0),)
  Sublattices     : 2
    Names         : (:A, :B)
    Sites         : (1, 1) --> 2 total per unit cell

julia> LatticePresets.honeycomb(semibounded = 1, names = (:C, :D))
Lattice{2,2,Float64} : 2D lattice in 2D space
  Bravais vectors : ((0.5, 0.866025), (-0.5, 0.866025))
    Semibounded   : (1)
  Sublattices     : 2
    Names         : (:C, :D)
    Sites         : (1, 1) --> 2 total per unit cell
```

# See also:
    LatticePresets, bravais, sublat, supercell, intracell
"""
lattice(s::Sublat, ss::Sublat...; kw...) = _lattice(Unitcell(s, ss...; kw...))
_lattice(u::Unitcell{E,T}) where {E,T} = Lattice(Bravais{E,T}(), u)
lattice(br::Bravais, s::Sublat, ss::Sublat...; kw...) = lattice(br, Unitcell(s, ss...; kw...))

function lattice(bravais::Bravais{E2,L2}, unitcell::Unitcell{E,T}) where {E2,L2,E,T}
    L = min(E,L2) # L should not exceed E
    Lattice(convert(Bravais{E,L,T}, bravais), unitcell)
end

issemibounded(lat::Lattice) = issemibounded(lat.bravais)

"""
    dims(lat::Lattice{E,L}) -> (E, L)

Return a tuple `(E, L)` of the embedding and lattice dimensions of `lat`
"""
dims(lat::Lattice{E,L}) where {E,L} = E, L

#######################################################################
# Supercell
#######################################################################
struct Supercell{L,L´,M<:Union{Missing,OffsetArray{Bool}},S<:SMatrix{L,L´}} # L´ is supercell dim
    matrix::S
    sites::UnitRange{Int}
    cells::CartesianIndices{L,NTuple{L,UnitRange{Int}}}
    mask::M
    semibounded::SVector{L´,Bool}
end

# Supercell{L}(ns::Integer) where {L} =
#     Supercell(SMatrix{L,L,Int,L*L}(I), 1:ns, CartesianIndices(ntuple(_->0:0, Val(L))),
#               missing, filltuple(false, Val(L)))

dim(::Supercell{L,L´}) where {L,L´} = L´

nsites(s::Supercell{L,L´,<:OffsetArray}) where {L,L´} = sum(s.mask)
nsites(s::Supercell{L,L´,Missing}) where {L,L´} = length(s.sites) * length(s.cells)

Base.CartesianIndices(s::Supercell) = s.cells

function Base.show(io::IO, s::Supercell{L,L´}) where {L,L´}
    hassemi = issemibounded(s)
    i = get(io, :indent, "")
    print(io, i,
"Supercell{$L,$(L´)} for $(L´)D superlattice of the base $(L)D lattice
$i  Supervectors  : $(displayvectors(s.matrix))
$i  Supersites    : $(nsites(s))", hassemi ? "
$i  Semibounded   : $(display_as_tuple(findall(s.semibounded)))" : "")
end

ismasked(s::Supercell{L,L´,<:OffsetArray})  where {L,L´} = true
ismasked(s::Supercell{L,L´,Missing})  where {L,L´} = false

isinmask(s::Supercell{L,L´,<:OffsetArray}, site, dn) where {L,L´} = s.mask[site, Tuple(dn)...]
isinmask(s::Supercell{L,L´,Missing}, site, dn) where {L,L´} = true
isinmask(s::Supercell{L,L´,<:OffsetArray}, site) where {L,L´} = s.mask[site]
isinmask(s::Supercell{L,L´,Missing}, site) where {L,L´} = true

issemibounded(sc::Supercell) = !iszero(sc.semibounded)

Base.copy(s::Supercell{<:Any,<:Any,Missing}) =
    Supercell(copy(s.matrix), s.sites, s.cells, s.mask, s.semibounded)

Base.copy(s::Supercell{<:Any,<:Any,<:AbstractArray}) =
    Supercell(copy(s.matrix), s.sites, s.cells, copy(s.mask), s.semibounded)

#######################################################################
# Superlattice
#######################################################################
struct Superlattice{E,L,T<:AbstractFloat,L´,S<:Supercell{L,L´},B<:Bravais{E,L,T},U<:Unitcell{E,T}} <: AbstractLattice{E,L,T}
    bravais::B
    unitcell::U
    supercell::S
end

function Base.show(io::IO, lat::Superlattice)
    i = get(io, :indent, "")
    hassemi = issemibounded(lat.bravais)
    ioindent = IOContext(io, :indent => string(i, "  "))
    print(io, i, summary(lat), "\n",
"$i  Bravais vectors : $(displayvectors(lat.bravais.matrix; digits = 6))", hassemi ? "
$i    Semibounded   : $(display_as_tuple(findall(lat.bravais.semibounded)))" : "", "
$i  Sublattices     : $(nsublats(lat))
$i    Names         : $(displaynames(lat))
$i    Sites         : $(display_as_tuple(sublatsites(lat))) --> $(nsites(lat)) total per unit cell\n")
    print(ioindent, lat.supercell)
end

Base.summary(::Superlattice{E,L,T,L´}) where {E,L,T,L´} =
    "Superlattice{$E,$L,$T,$L´} : $(L)D lattice in $(E)D space, filling a $(L´)D supercell"

# apply f to trues in mask. Arguments are s = sublat, oldi = old site, dn, newi = new site
function foreach_supersite(f::F, lat::Superlattice) where {F<:Function}
    newi = 0
    for s in 1:nsublats(lat), oldi in siterange(lat, s)
        for dn in CartesianIndices(lat.supercell)
            if isinmask(lat.supercell, oldi, dn)
                newi += 1
                f(s, oldi, toSVector(Int, Tuple(dn)), newi)
            end
        end
    end
    return nothing
end

issemibounded(lat::Superlattice) = issemibounded(lat.supercell)

#######################################################################
# AbstractLattice interface
#######################################################################

Base.copy(lat::Lattice) = Lattice(lat.bravais, copy(lat.unitcell))
Base.copy(lat::Superlattice) = Superlattice(lat.bravais, copy(lat.unitcell), copy(lat.supercell))

numbertype(::AbstractLattice{E,L,T}) where {E,L,T} = T
positiontype(::AbstractLattice{E,L,T}) where {E,L,T} = SVector{E,T}
dntype(::AbstractLattice{E,L}) where {E,L} = SVector{L,Int}

sublat(lat::AbstractLattice, siteidx) = sublat(lat.unitcell, siteidx)
sublats(lat::AbstractLattice) = sublats(lat.unitcell)

siterange(lat::AbstractLattice, sublat) = siterange(lat.unitcell, sublat)

siteindex(lat::AbstractLattice, sublat, idx) = siteindex(lat.unitcell, sublat, idx)

offsets(lat::AbstractLattice) = offsets(lat.unitcell)

sublatsites(lat::AbstractLattice) = sublatsites(lat.unitcell)

nsites(lat::AbstractLattice) = nsites(lat.unitcell)
nsites(lat::AbstractLattice, sublat) = nsites(lat.unitcell, sublat)

nsublats(lat::AbstractLattice) = nsublats(lat.unitcell)

issuperlattice(lat::Lattice) = false
issuperlattice(lat::Superlattice) = true

ismasked(lat::Lattice) = false
ismasked(lat::Superlattice) = ismasked(lat.supercell)

maskranges(lat::Superlattice) = (1:nsites(lat), lat.supercell.cells.indices...)
maskranges(lat::Lattice) = (1:nsites(lat),)

# External API #

"""
    sites(lat[, sublat::Int])

Extract the positions of all sites in a lattice, or in a specific sublattice
"""
sites(lat::AbstractLattice) = sites(lat.unitcell)
sites(lat::AbstractLattice, s) = sites(lat.unitcell, s)

"""
    transform!(lat::Lattice, f::Function)

Transform the site positions of `lat` by applying `f` to them in place.
"""
function transform!(lat::Lattice, f::Function)
    transform!(lat.unitcell, f)
    bravais´ = transform(lat.bravais, f)
    return Lattice(bravais´, lat.unitcell)
end

"""
    combine(lats::Lattice...)

If all `lats` have compatible Bravais vectors, combine them into a single lattice.
Sublattice names are renamed to be unique if necessary.
"""
function combine(lats::Lattice...)
    is_bravais_compatible(lats...) || throw(ArgumentError("Lattices must share all Bravais vectors, $(bravais.(lats))"))
    bravais´ = first(lats).bravais
    unitcell´ = combine((l -> l.unitcell).(lats)...)
    return Lattice(bravais´, unitcell´)
end

function combine(us::Vararg{Unitcell,N}) where {N}
    sites = vcat(ntuple(i -> us[i].sites, Val(N))...)
    names = uniquenames(tuplejoin(ntuple(i -> us[i].names, Val(N))...))
    offsets = combined_offsets(us...)
    return Unitcell(sites, names, offsets)
end

is_bravais_compatible() = true
is_bravais_compatible(lat::Lattice, lats::Lattice...) = all(l -> isequal(lat.bravais, l.bravais), lats)

function Base.isequal(b1::Bravais{E,L}, b2::Bravais{E,L}) where {E,L}
    vs1 = ntuple(i -> b1.matrix[:, i], Val(L))
    vs2 = ntuple(i -> b2.matrix[:, i], Val(L))
    for v2 in vs2
        found = false
        for v1 in vs1
            (isapprox(v1, v2) || isapprox(v1, -v2)) && (found = true; break)
        end
        !found && return false
    end
    return true
end

function combined_offsets(us::Unitcell...)
    nsubs = sum(nsublats, us)
    offsets = Vector{Int}(undef, nsubs + 1)
    lastoffset = 0
    idx = 1
    for u in us
        for idx´ in 1:(length(u.offsets) - 1)
            offsets[idx] = lastoffset + u.offsets[idx´]
            idx += 1
        end
        lastoffset = u.offsets[end]
    end
    offsets[end] = lastoffset + last(us).offsets[end]
    return offsets
end

#######################################################################
# supercell
#######################################################################
"""
    supercell(lat::AbstractLattice{E,L}, v::NTuple{L,Integer}...; region = missing)
    supercell(lat::AbstractLattice{E,L}, sc::SMatrix{L,L´,Int}; region = missing)

Generates a `Superlattice` from an `L`-dimensional lattice `lat` with Bravais vectors
`br´= br * sc`, where `sc::SMatrix{L,L´,Int}` is the integer supercell matrix with the `L´`
vectors `v`s as columns. If no `v` are given, the superlattice will be bounded. If `lat` has
semibounded axes, these cannot be mixed with any other axes (they can only be removed, kept
intact or scaled by a factor, see below).

Only sites at position `r` such that `region(r) == true` will be included in the supercell.
If `region` is missing, a Bravais unit cell perpendicular to the `v` axes will be selected
for the `L-L´` non-periodic directions.

    supercell(lattice::AbstractLattice{E,L}, factor::Integer; region = missing)

Calls `supercell` with a uniformly scaled `sc = SMatrix{L,L}(factor * I)`

    supercell(lattice::AbstractLattice, factors::Integer...; region = missing)

Calls `supercell` with different scaling along each Bravais vector (diagonal supercell
with factors along the diagonal)

    lat |> supercell(v...; kw...)

Functional syntax, equivalent to `supercell(lat, v...; kw...)

    supercell(h::Hamiltonian, v...; kw...)

Promotes the `Lattice` of `h` to a `Superlattice` without changing the Hamiltonian itself,
which always refers to the unitcell of the lattice.

# Examples
```jldoctest
julia> supercell(LatticePresets.honeycomb(), region = RegionPresets.circle(300))
Superlattice{2,2,Float64,0} : 2D lattice in 2D space, filling a 0D supercell
  Bravais vectors : ((0.5, 0.866025), (-0.5, 0.866025))
  Sublattices     : 2
    Names         : (:A, :B)
    Orbitals      : ((:noname,), (:noname,))
    Sites         : (1, 1) --> 2 total per unit cell
  Supercell{2,0} for 0D superlattice of the base 2D lattice
    Supervectors  : ()
    Supersites    : 652966

julia> supercell(LatticePresets.triangular(), (1,1), (1, -1))
Superlattice{2,2,Float64,2} : 2D lattice in 2D space, filling a 2D supercell
  Bravais vectors : ((0.5, 0.866025), (-0.5, 0.866025))
  Sublattices     : 1
    Names         : (:A)
    Orbitals      : ((:noname,),)
    Sites         : (1) --> 1 total per unit cell
  Supercell{2,2} for 2D superlattice of the base 2D lattice
    Supervectors  : ((1, 1), (1, -1))
    Supersites    : 2

julia> LatticePresets.square() |> supercell(3)
Superlattice{2,2,Float64,2} : 2D lattice in 2D space, filling a 2D supercell
  Bravais vectors : ((1.0, 0.0), (0.0, 1.0))
  Sublattices     : 1
    Names         : (:A)
    Orbitals      : ((:noname,),)
    Sites         : (1) --> 1 total per unit cell
  Supercell{2,2} for 2D superlattice of the base 2D lattice
    Supervectors  : ((3, 0), (0, 3))
    Supersites    : 9
```

# See also:
    unitcell
"""
supercell(v::Union{SMatrix,Tuple,SVector,Integer}...; kw...) = lat -> supercell(lat, v...; kw...)

supercell(lat::AbstractLattice{E,L}; kw...) where {E,L} =
    supercell(lat, SMatrix{L,0,Int}(); kw...)
supercell(lat::AbstractLattice{E,L}, factors::Vararg{<:Integer,L}; kw...) where {E,L} =
    _supercell(lat, factors...)
supercell(lat::AbstractLattice{E,L}, factors::Vararg{<:Integer,L´}; kw...) where {E,L,L´} =
    throw(ArgumentError("Provide either a single scaling factor or one for each of the $L lattice dimensions"))
supercell(lat::AbstractLattice{E,L}, factor::Integer; kw...) where {E,L} =
    _supercell(lat, ntuple(_ -> factor, Val(L))...)
supercell(lat::AbstractLattice, vecs::NTuple{L,Integer}...; region = missing) where {L} =
    _supercell(lat, toSMatrix(Int, vecs...), region)
supercell(lat::AbstractLattice, s::SMatrix; region = missing) = _supercell(lat, s, region)

function _supercell(lat::AbstractLattice{E,L}, factors::Vararg{Integer,L}) where {E,L}
    scmatrix = SMatrix{L,L,Int}(Diagonal(SVector(factors)))
    sites = 1:nsites(lat)
    cells = CartesianIndices((i -> 0 : i - 1).(factors))
    mask = missing
    semibounded = supercell_semibounded(lat, scmatrix)
    supercell = Supercell(scmatrix, sites, cells, mask, semibounded)
    return Superlattice(lat.bravais, lat.unitcell, supercell)
end

function _supercell(lat::AbstractLattice{E,L}, scmatrix::SMatrix{L,L´,Int}, region) where {E,L,L´}
    semibounded = supercell_semibounded(lat, scmatrix)
    brmatrix = lat.bravais.matrix
    regionfunc = region === missing ? ribbonfunc(brmatrix, scmatrix) : region
    in_supercell_func = is_perp_dir(scmatrix)
    cells = supercell_cells(lat, regionfunc, in_supercell_func)
    ns = nsites(lat)
    mask = OffsetArray(BitArray(undef, ns, size(cells)...), 1:ns, cells.indices...)
    @inbounds for dn in cells
        dntup = Tuple(dn)
        dnvec = toSVector(Int, dntup)
        in_supercell = in_supercell_func(dnvec)
        in_supercell || (mask[:, dntup...] .= false; continue)
        r0 = brmatrix * dnvec
        for (i, site) in enumerate(lat.unitcell.sites)
            r = site + r0
            mask[i, dntup...] = in_supercell && regionfunc(r)
        end
    end
    supercell = Supercell(scmatrix, 1:ns, cells, all(mask) ? missing : mask, semibounded)
    return Superlattice(lat.bravais, lat.unitcell, supercell)
end

# This is true whenever old ndist is perpendicular to new lattice
is_perp_dir(supercell) = let invs = pinvmultiple(supercell); dn -> iszero(new_dn(dn, invs)); end

new_dn(oldndist, (pinvs, n)) = fld.(pinvs * oldndist, n)
new_dn(oldndist, ::Tuple{<:SMatrix{0,0},Int}) = SVector{0,Int}()

function ribbonfunc(bravais::SMatrix{E,L,T}, supercell::SMatrix{L,L´}) where {E,L,T,L´}
    L <= L´ && return r -> true
    # real-space supercell axes + all space
    s = hcat(bravais * supercell, SMatrix{E,E,T}(I))
    q = qr(s).Q
    # last vecs in Q are orthogonal to axes
    os = ntuple(i -> SVector(q[:,i+L´]), Val(E-L´))
    # range of projected bravais, including zero
    ranges = (o -> extrema(vcat(SVector(zero(T)), bravais' * o)) .- sqrt(eps(T))).(os)
    # projector * r gives the projection of r on orthogonal vectors
    projector = hcat(os...)'
    # ribbon defined by r's with projection within ranges for all orthogonal vectors
    regionfunc(r) = all(first.(ranges) .<= Tuple(projector * r) .< last.(ranges))
    return regionfunc
end

function supercell_cells(lat::Lattice{E,L}, regionfunc, in_supercell_func) where {E,L}
    bravais = lat.bravais.matrix
    iter = BoxIterator(zero(SVector{L,Int}))
    foundfirst = false
    counter = 0
    for dn in iter
        found = false
        counter += 1; counter == TOOMANYITERS &&
            throw(ArgumentError("`region` seems unbounded (after $TOOMANYITERS iterations)"))
        in_supercell = in_supercell_func(toSVector(Int, dn))
        r0 = bravais * toSVector(Int, dn)
        for site in lat.unitcell.sites
            r = r0 + site
            found = in_supercell && regionfunc(r)
            if found || !foundfirst
                acceptcell!(iter, dn)
                foundfirst = found
                break
            end
        end
    end
    c = CartesianIndices(iter)
    return c
end

supercell_semibounded(lat::Lattice, s::SMatrix) = supercell_semibounded(lat.bravais, s)
function supercell_semibounded(b::Bravais, s::SMatrix{L,L´,Int}) where {L,L´}
    lsb = b.semibounded
    err = ArgumentError("Cannot mix semibounded axes in a supercell")
    ssb = ntuple(Val(L´)) do i
        vec_semi = s[lsb, i]
        vec_other = s[(!).(lsb), i]
        if iszero(vec_other)
            count(!iszero, vec_semi) != 1 && throw(err)
            return true
        else
            iszero(vec_semi) || throw(err)
            return false
        end
    end
    return SVector{L´,Bool}(ssb)
end

#######################################################################
# unitcell
#######################################################################
"""
    unitcell(lat::Lattice{E,L}, v::NTuple{L,Integer}...; region = missing)
    unitcell(lat::Lattice{E,L}, uc::SMatrix{L,L´,Int}; region = missing)

Generates a `Lattice` from an `L`-dimensional lattice `lat` and a larger unit cell, such
that its Bravais vectors are `br´= br * uc`. Here `uc::SMatrix{L,L´,Int}` is the integer
unitcell matrix, with the `L´` vectors `v`s as columns. If no `v` are given, the new lattice
will be bounded.

Only sites at position `r` such that `region(r) == true` will be included in the new
unitcell. If `region` is missing, a Bravais unitcell perpendicular to the `v` axes will be
selected for the `L-L´` non-periodic directions.

    unitcell(lattice::Lattice{E,L}, factor::Integer; region = missing)

Calls `unitcell` with a uniformly scaled `uc = SMatrix{L,L}(factor * I)`

    unitcell(lattice::Lattice{E,L}, factors::Integer...; region = missing)

Calls `unitcell` with different scaling along each Bravais vector (diagonal supercell
with factors along the diagonal)

    unitcell(slat::Superlattice)

Convert Superlattice `slat` into a lattice with its unit cell matching `slat`'s supercell.

    unitcell(h::Hamiltonian, v...; modifiers = (), kw...)

Transforms the `Lattice` of `h` to have a larger unitcell, while expanding the Hamiltonian
accordingly. The modifiers (a tuple of `ElementModifier`s, either `onsite!` or `hopping!`)
will be applied to onsite and hoppings as the hamiltonian is expanded. See `onsite!` and
`hopping!` for details

Note: for performance reasons, in sparse hamiltonians only the stored onsites and hoppings
will be transformed by `ElementModifier`s, so you might want to add zero onsites or hoppings
when building `h` to have a modifier applied to them later. Note also that additional
onsites and hoppings may be stored when calling `optimize!` or `bloch`/`bloch!` on `h` for
the first time.

    lat_or_h |> unitcell(v...; kw...)

Functional syntax, equivalent to `unitcell(lat_or_h, v...; kw...)

# Examples
```jldoctest
julia> unitcell(LatticePresets.honeycomb(), region = RegionPresets.circle(300))
Lattice{2,0,Float64} : 0D lattice in 2D space
  Bravais vectors : ()
  Sublattices     : 2
    Names         : (:A, :B)
    Orbitals      : ((:noname,), (:noname,))
    Sites         : (326483, 326483) --> 652966 total per unit cell

julia> unitcell(LatticePresets.triangular(), (1,1), (1, -1))
Lattice{2,2,Float64} : 2D lattice in 2D space
  Bravais vectors : ((0.0, 1.732051), (1.0, 0.0))
  Sublattices     : 1
    Names         : (:A)
    Orbitals      : ((:noname,),)
    Sites         : (2) --> 2 total per unit cell

julia> LatticePresets.square() |> unitcell(3)
Lattice{2,2,Float64} : 2D lattice in 2D space
  Bravais vectors : ((3.0, 0.0), (0.0, 3.0))
  Sublattices     : 1
    Names         : (:A)
    Orbitals      : ((:noname,),)
    Sites         : (9) --> 9 total per unit cell

julia> supercell(LatticePresets.square(), 3) |> unitcell
Lattice{2,2,Float64} : 2D lattice in 2D space
  Bravais vectors : ((3.0, 0.0), (0.0, 3.0))
  Sublattices     : 1
    Names         : (:A)
    Orbitals      : ((:noname,),)
    Sites         : (9) --> 9 total per unit cell
```

# See also:
    supercell
"""
unitcell(v::Union{SMatrix,Tuple,SVector,Integer}...; kw...) = lat -> unitcell(lat, v...; kw...)
unitcell(lat::Lattice, args...; kw...) = unitcell(supercell(lat, args...; kw...))

function unitcell(lat::Superlattice)
    newoffsets = supercell_offsets(lat)
    newsites = supercell_sites(lat)
    unitcell = Unitcell(newsites, lat.unitcell.names, newoffsets)
    bravais = Bravais(lat.bravais.matrix * lat.supercell.matrix, lat.supercell.semibounded)
    return Lattice(bravais, unitcell)
end

function supercell_offsets(lat::Superlattice)
    sitecounts = zeros(Int, nsublats(lat) + 1)
    foreach_supersite((s, oldi, dn, newi) -> sitecounts[s + 1] += 1, lat)
    newoffsets = cumsum!(sitecounts, sitecounts)
    return newoffsets
end

function supercell_sites(lat::Superlattice)
    newsites = similar(lat.unitcell.sites, nsites(lat.supercell))
    oldsites = lat.unitcell.sites
    bravais = lat.bravais.matrix
    foreach_supersite((s, oldi, dn, newi) -> newsites[newi] = bravais * dn + oldsites[oldi], lat)
    return newsites
end

# #######################################################################
# # semibounded
# #######################################################################
# """
#     semibounded(lat::AbstractLattice, axes)

# Create an AbstractLattice like `lat` but with a specific set of semibound axes. These can be
# a boolean (i.e. all or none are semibounded), a tuple of booleans (one per axis) or a
# collection of axis indices to make semibounded

#     lat |> semibounded(axes)

# Function form equivalent to semibounded(lat, axes)

# # Examples
# ```jldoctest
# julia> semibounded(LatticePresets.honeycomb(), 2)
# Lattice{2,2,Float64} : 2D lattice in 2D space
#   Bravais vectors : ((0.5, 0.866025), (-0.5, 0.866025))
#     Semibounded   : (2)
#   Sublattices     : 2
#     Names         : (:A, :B)
#     Sites         : (1, 1) --> 2 total per unit cell

# julia> LatticePresets.honeycomb() |> supercell(4) |> semibounded((true, false))
# Superlattice{2,2,Float64,2} : 2D lattice in 2D space, filling a 2D supercell
#   Bravais vectors : ((0.5, 0.866025), (-0.5, 0.866025))
#   Sublattices     : 2
#     Names         : (:A, :B)
#     Sites         : (1, 1) --> 2 total per unit cell
#   Supercell{2,2} for 2D superlattice of the base 2D lattice
#     Supervectors  : ((4, 0), (0, 4))
#     Supersites    : 32
#     Semibounded   : (1)
# ```
# """
# semibounded(axes) = lat -> semibounded(lat, axes)

# function semibounded(lat::Lattice, axes = true)
#     matrix = lat.bravais.matrix
#     sb = sanitize_semibounded(axes, matrix)
#     br = Bravais(matrix, sb)
#     return Lattice(br, lat.unitcell)
# end

# function semibounded(lat::Superlattice, axes = true)
#     sc = lat.supercell
#     matrix = lat.bravais.matrix * sc.matrix
#     sb = sanitize_semibounded(axes, matrix)
#     scell = Supercell(sc.matrix, sc.sites, sc.cells, sc.mask, sb)
#     return Superlattice(lat.bravais, lat.unitcell, scell)
# end

max_issues_repo_name: pablosanjose/QBox

max_forks_repo_name: pablosanjose/QBox


text:
include("cg.jl")
using CG;
include("camera.jl")
include("sceneobjects.jl")
include("raytracing.jl")
include("lighting.jl");
include("shaders.jl")
include("render.jl")
include("objects.jl")
using Objects;
using Camera;
using Raytracing;
using Shaders;
using Lighting;

cam = PinholeCamera(Vec4f(0, 0, 0, 1), Vec4f(0, 0, -1, 1), Vec4f(0, 1, 0, 1));
sphere1 = Sphere(Vec4f(-0.5, 0.5, 0, 1), 0.25f0);
sphere2 = Sphere(Vec4f(-0.5, -0.5, 0, 1), 0.5f0);
aabb1 = AABB(Vec4f(0.5, -0.5, 0, 1), 0.25f0, 0.25f0, 0.25f0);
aabb2 = AABB(Vec4f(0.5, 0.5, 0, 1), 0.5f0, 0.5f0, 0.5f0);

light1 = PointLight(Vec4f(0.9, -0.5, 0.3, 1));
light2 = PointLight(Vec4f(-1, 0, 5, 1));

scene = Scene(SceneObject[sphere1, sphere2, aabb1, aabb2]);
lights = SceneLights([light1, light2]);
tracerays(scene, lights, cam, blinnPhongShader);
show();

ext: jl

max_stars_repo_name: fpoppinga/cg

max_issues_repo_name: fpoppinga/cg

max_forks_repo_name: fpoppinga/cg


text:
module Kagome

using LatticeTools

NN_BOND_TYPES = [
    ([ 0, 0], "A", [ 0, 0], "B", 1),
    ([ 0, 0], "B", [ 0, 0], "C", 1),
    ([ 0, 0], "C", [ 0, 0], "A", 1),

    ([ 1, 1], "A", [ 1, 0], "B",-1),
    ([ 1, 0], "B", [ 0, 0], "C",-1),
    ([ 0, 0], "C", [ 1, 1], "A",-1),
]

NNN_BOND_TYPES = [
    ([ 0, 0], "A", [ 1, 0], "B", 1), # ◁
    ([ 1, 0], "B", [ 0,-1], "C", 1),
    ([ 0,-1], "C", [ 0, 0], "A", 1),
    ([ 0, 0], "C", [ 1, 0], "A",-1),
    ([ 1, 0], "A", [ 0,-1], "B",-1), # ▷
    ([ 0,-1], "B", [ 0, 0], "C",-1),
]

export make_kagome_lattice
export make_kagome_obc

function make_kagome_lattice(size_matrix ::AbstractMatrix{<:Integer}; compute_symmetry::Bool=false)
    latticevectors = [1 -0.5; 0 0.5*sqrt(3.0)];
    unitcell = makeunitcell(latticevectors, SiteType=String)
    addsite!(unitcell, "A", carte2fract(unitcell, [0.5, 0.0]))
    addsite!(unitcell, "B", carte2fract(unitcell, [0.25, 0.25*sqrt(3.0)]))
    addsite!(unitcell, "C", carte2fract(unitcell, [0.5+0.25, 0.25*sqrt(3.0)]))

    lattice = make_lattice(unitcell, size_matrix)
    hypercube = lattice.hypercube
    supercell = lattice.supercell
    tsym = FiniteTranslationSymmetry(lattice)
    psym = little_symmetry(tsym, project(PointSymmetryDatabase.get(25), [1 0 0; 0 1 0]))
    tsymbed = embed(lattice, tsym)
    psymbed = embed(lattice, psym)
    ssymbed = SymmorphicSymmetry(tsymbed, psymbed)

    nnbonds = []
    nnnbonds = []

    for r in lattice.bravais_coordinates
        for (rowvec, rowsite, colvec, colsite, bondsign) in NN_BOND_TYPES
            R_row, r_row = hypercube.wrap(r .+ rowvec)
            R_col, r_col = hypercube.wrap(r .+ colvec)
            rowsite_super = (rowsite, r_row)
            colsite_super = (colsite, r_col)
            irow = get(supercell.siteindices, rowsite_super, -1)
            icol = get(supercell.siteindices, colsite_super, -1)
            push!(nnbonds, ((irow, icol), R_col-R_row, bondsign))
        end
        for (rowvec, rowsite, colvec, colsite, bondsign) in NNN_BOND_TYPES
            R_row, r_row = hypercube.wrap(r .+ rowvec)
            R_col, r_col = hypercube.wrap(r .+ colvec)
            rowsite_super = (rowsite, r_row)
            colsite_super = (colsite, r_col)
            irow = get(supercell.siteindices, rowsite_super, -1)
            icol = get(supercell.siteindices, colsite_super, -1)
            push!(nnnbonds, ((irow, icol), R_col-R_row, bondsign))
        end
    end

    nn_triangles = []
    for r in lattice.bravais_coordinates
      triangle = []
      for (rowvec, roworb, colvec, colorb, bondsign) in NN_BOND_TYPES[1:3]
        R_row, r_row = hypercube.wrap(r .+ rowvec)
        R_col, r_col = hypercube.wrap(r .+ colvec)
        roworb_super = (roworb, r_row)
        colorb_super = (colorb, r_col)
        irow = get(supercell.siteindices, roworb_super, -1)
        icol = get(supercell.siteindices, colorb_super, -1)
        push!(triangle, ((irow, icol), R_col-R_row))
      end
      push!(nn_triangles, (triangle, 1))

      triangle = []
      for (rowvec, roworb, colvec, colorb, bondsign) in NN_BOND_TYPES[4:6]
        R_row, r_row = hypercube.wrap(r .+ rowvec)
        R_col, r_col = hypercube.wrap(r .+ colvec)
        roworb_super = (roworb, r_row)
        colorb_super = (colorb, r_col)
        irow = get(supercell.siteindices, roworb_super, -1)
        icol = get(supercell.siteindices, colorb_super, -1)
        push!(triangle, ((irow, icol), R_col-R_row))
      end
      push!(nn_triangles, (triangle, -1))
    end

    nnn_triangles = []
    for r in lattice.bravais_coordinates
      triangle = []
      for (rowvec, roworb, colvec, colorb, bondsign) in NNN_BOND_TYPES[1:3]
        R_row, r_row = hypercube.wrap(r .+ rowvec)
        R_col, r_col = hypercube.wrap(r .+ colvec)
        roworb_super = (roworb, r_row)
        colorb_super = (colorb, r_col)
        irow = get(supercell.siteindices, roworb_super, -1)
        icol = get(supercell.siteindices, colorb_super, -1)
        push!(triangle, ((irow, icol), R_col-R_row))
      end
      push!(nnn_triangles, (triangle, 1))

      triangle = []
      for (rowvec, roworb, colvec, colorb, bondsign) in NNN_BOND_TYPES[4:6]
        R_row, r_row = hypercube.wrap(r .+ rowvec)
        R_col, r_col = hypercube.wrap(r .+ colvec)
        roworb_super = (roworb, r_row)
        colorb_super = (colorb, r_col)
        irow = get(supercell.siteindices, roworb_super, -1)
        icol = get(supercell.siteindices, colorb_super, -1)
        push!(triangle, ((irow, icol), R_col-R_row))
      end
      push!(nnn_triangles, (triangle, -1))
    end

    return (unitcell=unitcell,
            lattice=lattice,
            space_symmetry_embedding=ssymbed,
            nearest_neighbor_bonds=nnbonds,
            next_nearest_neighbor_bonds=nnnbonds,
            nearest_neighbor_triangles=nn_triangles,
            next_nearest_neighbor_triangles=nnn_triangles,)
end

function make_kagome_obc(sites::AbstractVector{<:Tuple})
    site_set = Set(sites)

    min_R1 = minimum(R[1] for (sitetype, R) in sites)
    max_R1 = maximum(R[1] for (sitetype, R) in sites)
    min_R2 = minimum(R[2] for (sitetype, R) in sites)
    max_R2 = maximum(R[2] for (sitetype, R) in sites)

    n1 = max_R1 - min_R1 + 1
    n2 = max_R2 - min_R2 + 1

    big_kagome = make_kagome_lattice([n1+1 0; 0 n2+1])

    site_folding = Dict{Int, Int}()

    new_kagome_supercell = makeunitcell(
        big_kagome.lattice.supercell.latticevectors,
        SiteType=Tuple{String, Vector{Int}}
    )

    for (siteindex, ((sitetype, siteuc), sitecoord)) in enumerate(big_kagome.lattice.supercell.sites)
        if (sitetype, siteuc + [min_R1, min_R2]) in site_set
            addsite!(new_kagome_supercell, (sitetype, siteuc + [min_R1, min_R2]), sitecoord)
            site_folding[siteindex] = numsite(new_kagome_supercell)
        end
    end

    nn_bonds = Tuple{Int, Int}[]
    for ((i, j), R, _) in big_kagome.nearest_neighbor_bonds
        if !haskey(site_folding, i) || !haskey(site_folding, j) || !iszero(R)
            continue
        end
        ip = site_folding[i]
        jp = site_folding[j]
        if ip < jp
            push!(nn_bonds, (ip, jp))
        else
            push!(nn_bonds, (jp, ip))
        end
    end

    #@show big_kagome.nearest_neighbor_triangles

    nn_triangles = Tuple{Int, Int, Int}[]
    for (tri, sgn) in big_kagome.nearest_neighbor_triangles
        triangle_vertices = Set{Int}()
        complete_triangle = true
        for ((i, j), R) in tri
            if !iszero(R)
                complete_triangle = false
                break
            end
            ip = get(site_folding, i, -1)
            jp = get(site_folding, j, -1)

            if ip <= 0 || jp <= 0
                complete_triangle = false
                break
            end

            push!(triangle_vertices, ip)
            push!(triangle_vertices, jp)
        end

        !complete_triangle && continue

        push!(nn_triangles, tuple(sort(collect(triangle_vertices))...))
    end

    return (unitcell=new_kagome_supercell,
            nearest_neighbor_bonds=nn_bonds,
            nearest_neighbor_triangles=nn_triangles)

end

end # module Kagome


text:
using Test
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

@testset "Viewer.jl" begin

	 # function Viewer(meshes)
   @testset "Viewer" begin
      @test
      @test
      @test
      @test
   end

	 # function releaseGpuResources(viewer::Viewer)
   @testset "releaseGpuResources" begin
      @test
      @test
      @test
      @test
   end

	 # function runViewer(viewer::Viewer)
   @testset "runViewer" begin
      @test
      @test
      @test
      @test
   end

	 # function VIEW(meshes)
   @testset "VIEW" begin
      @test
      @test
      @test
      @test
   end

	 # function getModelview(viewer::Viewer)
   @testset "getModelview" begin
      @test
      @test
      @test
      @test
   end

	 # function getProjection(viewer::Viewer)
   @testset "getProjection" begin
      @test
      @test
      @test
      @test
   end

	 # function projectPoint(viewer::Viewer,pos::Point3d)
   @testset "projectPoint" begin
      @test
      @test
      @test
      @test
   end

	 # function unprojectPoint(viewer::Viewer,x::Float32,y::Float32)
   @testset "unprojectPoint" begin
      @test
      @test
      @test
      @test
   end

	 # function getShader(viewer::Viewer,lighting_enabled,color_attribute_enabled)
   @testset "getShader" begin
      @test
      @test
      @test
      @test
   end

	 # function glRender(viewer::Viewer)
   @testset "glRender" begin
      @test
      @test
      @test
      @test
   end

	 # function redisplay(viewer::Viewer)
   @testset "redisplay" begin
      @test
      @test
      @test
      @test
   end

	 # function handleResizeEvent(viewer)
   @testset "handleResizeEvent" begin
      @test
      @test
      @test
      @test
   end

	 # function handleMouseMoveEvent(viewer,x,y)
   @testset "handleMouseMoveEvent" begin
      @test
      @test
      @test
      @test
   end

	 # function handleMouseWheelEvent(viewer,delta)
   @testset "handleMouseWheelEvent" begin
      @test
      @test
      @test
      @test
   end

	 # function handleKeyPressEvent(viewer,key, scancode, action, mods)
   @testset "handleKeyPressEvent" begin
      @test
      @test
      @test
      @test
   end

end


text:
module Adjacent

export prevFloat, nextFloat, prevprevFloat, nextnextFloat,
       nextAwayFromZero, nextnextAwayFromZero,
       nextNearerToZero, nextnextNearerToZero

# exact for |x| > 8.900295434028806e-308
nextNearerToZero(x::Float64)   = (x-1.1102230246251568e-16*x)-5.0e-324  # do not simplify
nextAwayFromZero(x::Float64)   = (x+1.1102230246251568e-16*x)+5.0e-324  # do not simplify
nextnextNearerToZero(x::Float64)   = nextNearerToZero(nextNearerToZero(x)) # changing the multiplier does not work
nextnextAwayFromZero(x::Float64)   = nextAwayFromZero(nextAwayFromZero(x))

# exact for |x| > 4.814825f-35
nextNearerToZero(x::Float32)   = (x-5.960465f-8*x)-1.435f-42 # do not simplify
nextAwayFromZero(x::Float32)   = (x+5.960465f-8*x)+1.435f-42 # do not simplify
nextnextNearerToZero(x::Float64)   = nextNearerToZero(nextNearerToZero(x)) # changing the multiplier does not work
nextnextAwayFromZero(x::Float64)   = nextAwayFromZero(nextAwayFromZero(x))

# the multiplicative formulation for Float16 is exact for |x| > Float16(0.25)
# which is quite coarse, we do not use that here
nextNearerToZero(x::Float16) = (x < 0) ? nextfloat(x) : prevfloat(x)
nextAwayFromZero(x::Float16) = (x < 0) ? prevfloat(x) : nextfloat(x)
nextFloat(x::Float16) = nextfloat(x)
prevFloat(x::Float16) = prevfloat(x)

nextFloat{T<:AbstractFloat}(x::T) = signbit(x) ? nextNearerToZero(x) : nextAwayFromZero(x)
prevFloat{T<:AbstractFloat}(x::T) = signbit(x) ? nextAwayFromZero(x) : nextNearerToZero(x)

nextnextFloat{T<:AbstractFloat}(x::T) = signbit(x) ? nextnextNearerToZero(x) : nextnextAwayFromZero(x)
prevprevFloat{T<:AbstractFloat}(x::T) = signbit(x) ? nextnextAwayFromZero(x) : nextnextNearerToZero(x)

end # module


text:
module MCMCChains

using AxisArrays
const axes = Base.axes

import AbstractMCMC
import AbstractMCMC: chainscat
using Distributions
using RecipesBase
using SpecialFunctions
using Formatting
import StatsBase: autocov, counts, sem, AbstractWeights,
    autocor, describe, quantile, sample, summarystats, cov
using Requires

using LinearAlgebra: diag
import Serialization: serialize, deserialize
import Random
import Statistics: std, cor, mean, var

export Chains, chains, chainscat
export set_section, get_params, sections, sort_sections, setinfo, set_names
export autocor, describe, sample, summarystats, AbstractWeights, mean, quantile
export ChainDataFrame, DataFrame
export summarize

# Export diagnostics functions
export discretediag, gelmandiag, gewekediag, heideldiag, rafterydiag
export hpd, ess

"""
    Chains

Parameters:

- `value`: An `AxisArray` object with axes `iter` × `var` × `chains`
- `logevidence` : A field containing the logevidence.
- `name_map` : A `NamedTuple` mapping each variable to a section.
- `info` : A `NamedTuple` containing miscellaneous information relevant to the chain.
The `info` field can be set using `setinfo(c::Chains, n::NamedTuple)`.
"""
struct Chains{A, T, K<:NamedTuple, L<:NamedTuple} <: AbstractMCMC.AbstractChains
    value::AxisArray{A,3}
    logevidence::T
    name_map::K
    info::L
end

include("utils.jl")
include("chains.jl")
include("constructors.jl")
include("summarize.jl")
include("discretediag.jl")
include("fileio.jl")
include("gelmandiag.jl")
include("gewekediag.jl")
include("heideldiag.jl")
include("mcse.jl")
include("rafterydiag.jl")
include("sampling.jl")
include("stats.jl")
include("modelstats.jl")
include("plot.jl")

function __init__()
    @require DataFrames="a93c6f00-e57d-5684-b7b6-d8193f3e46c0" include("dataframes-compat.jl")
end

end # module


text:
module GenericBoundaryTriangulationsTests

using Test
using Gridap.TensorValues
using Gridap.Helpers
using Gridap.Arrays
using Gridap.Geometry
using Gridap.Fields
using Gridap.ReferenceFEs
using FillArrays

using Gridap.Geometry: FaceCellCoordinates
using Gridap.Geometry: FaceToCellMapValue
using Gridap.Geometry: CellShapeFunsAtFaces

domain = (0,4,0,4)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)
model = UnstructuredDiscreteModel(model)
topo = get_grid_topology(model)

btrian = GenericBoundaryTriangulation(model,get_isboundary_face(topo,1))
test_triangulation(btrian)

s = CompressedArray([Point{1,Float64}[(0.25,),(0.75,)]],get_cell_type(btrian))

face_to_fvertex_to_qcoors = FaceCellCoordinates(btrian)
r =  Vector{Point{2,Float64}}[
  [(0.0,0.0),(1.0,0.0)],[(0.0,0.0),(0.0,1.0)],
  [(0.0,0.0),(1.0,0.0)],[(1.0,0.0),(1.0,1.0)],
  [(0.0,1.0),(1.0,1.0)],[(0.0,0.0),(0.0,1.0)],
  [(0.0,1.0),(1.0,1.0)],[(1.0,0.0),(1.0,1.0)]]
test_array(face_to_fvertex_to_qcoors,r)

s2q = get_face_to_cell_map(btrian)

q = evaluate(s2q,s)
r = Vector{Point{2,Float64}}[
  [(0.25,0.0),(0.75,0.0)],[(0.0,0.25),(0.0,0.75)],
  [(0.25,0.0),(0.75,0.0)],[(1.0,0.25),(1.0,0.75)],
  [(0.25,1.0),(0.75,1.0)],[(0.0,0.25),(0.0,0.75)],
  [(0.25,1.0),(0.75,1.0)],[(1.0,0.25),(1.0,0.75)]]
test_array(q,r)
@test isa(q,FaceToCellMapValue)

q2x = reindex(get_cell_map(get_grid(model)), btrian.glue.face_to_cell)

s2x = compose(q2x,s2q)

x = evaluate(s2x,s)
r = Vector{Point{2,Float64}}[
  [(0.5,0.0),(1.5,0.0)],[(0.0,0.5),(0.0,1.5)],
  [(2.5,0.0),(3.5,0.0)],[(4.0,0.5),(4.0,1.5)],
  [(0.5,4.0),(1.5,4.0)],[(0.0,2.5),(0.0,3.5)],
  [(2.5,4.0),(3.5,4.0)],[(4.0,2.5),(4.0,3.5)]]
test_array(x,r)

cell_shapefuns_q = reindex(get_cell_shapefuns(get_grid(model)), btrian.glue.face_to_cell)
cell_shapefuns_s = compose(cell_shapefuns_q,s2q)
@test isa(evaluate(cell_shapefuns_s,s), CellShapeFunsAtFaces)

nvec = get_normal_vector(btrian)

nvec_s = evaluate(nvec,s)
r = Vector{Point{2,Float64}}[
  [(0.0,-1.0),(0.0,-1.0)],[(-1.0,0.0),(-1.0,0.0)],
  [(0.0,-1.0),(0.0,-1.0)],[(1.0,-0.0),(1.0,-0.0)],
  [(0.0,1.0),(0.0,1.0)],[(-1.0,0.0),(-1.0,0.0)],
  [(0.0,1.0),(0.0,1.0)],[(1.0,-0.0),(1.0,-0.0)]]
test_array(nvec_s,r)

end # module


text:
struct QPSolveFailure <: Exception
    terminationstatus
    primalstatus
    dualstatus
end

Base.showerror(io::IO, err::QPSolveFailure) = print(io, """
    QP solve unsuccessful.
    Termination status: $(err.terminationstatus)
    Primal status: $(err.primalstatus)
    Dual status: $(err.dualstatus)""")


text:
# This file calculates the UL Decomposition of a Toeplitz operator,
# via the Wiener–Hopf factorization.  If the symbol of
# the Toeplitz operator is G(z), we find a Wiener–Hopf factorization
#
# 	G(z) = φ_+(z)*φ_-(z)
#
# where φ_± are analytic in the interior/exterior of the unit
# circle.  Then we have the UL decomposition
#
#   T[G]=T[φ_+]*T[φ_-]
#
# The Wiener–Hopf factorization is found by rephrasing as a
# Riemann–Hilbert problem:
#
#   φ_+(z)=φ_-(z)^(-1)*G(z)
#
# and writing
#
#   φ_-(z)^(-1) = 1 + C_-V(z) and φ_+(z)=1+C_+V(z)=1+V+C_-V(z)
#
# where C_± are the left/right limits of the Cauchy operator.
# This leads to a singular integral equation
#
#   V + C_-V*(1-G) = G - 1
#
# This approach also extends to the matrix case.

using ApproxFun, SingularIntegralEquations

# Scalar case

G=Fun(z->2+cos(z+1/z),Circle()) # the symbol of the Toeplitz operator
T=ToeplitzOperator(G)

C  = Cauchy(-1)
V  = (I+(1-G)*C)\(G-1)

Φmi = 1+C*V
Φp = V+Φmi

L  = ToeplitzOperator(1/Φmi)
U  = ToeplitzOperator(Φp)

norm((T-U*L)[1:10,1:10])  # check the accuracy

# Matrix case

G=Fun(z->[-1 -3; -3 -1]/z +
         [ 2  2;  1 -3] +
         [ 2 -1;  1  2]*z,Circle())
T=ToeplitzOperator(G)

C  = Cauchy(-1)
V  = V=(I+(I-G)*C)\(G-I)

Φmi = I+C*V
Φp = V+Φmi

L  = ToeplitzOperator(inv(Φmi))
U  = ToeplitzOperator(Φp)

norm((T-U*L)[1:10,1:10])  # check the accuracy


text:
"""
read all WinMF *.mat files in current folder and convert them to log scaled csv
frequency responses
"""

using MAT, Dierckx, DelimitedFiles, myutils

function logResample(data1::Vector, samplerate::Number, fLog=logrange(10, 23e3, points_per_octave=96))
	nData = length(data1)
	fData = range(0, stop=samplerate/2, length=nData)
	a=3
	n=((fLog[end]-fLog[end-1])/(fData[2]-fData[1])*2*a) |>round|>Int
	x = range(-1, stop=1, length=n)
	lanczos = sinc.(x).*sinc.(x*a)
	lanczos ./= sum(lanczos)
	function lanczosInterpolate(f0, d, itp)
		fAxis = range(f0-a*d, stop=f0+a*d, length=n)
		pointValues = itp(fAxis)
		return sum(lanczos .* pointValues)
	end

	logData = Array{Complex{Float64}, 1}(undef, length(fLog))

	itpReal = Spline1D(fData, real.(data1))
	itpImag = Spline1D(fData, imag.(data1))

	for i in 1:length(fLog)
		f0 = fLog[i]
		if i<length(fLog)
			d = fLog[i+1]-f0
		else
			d = -fLog[i-1]+f0
		end
		logData[i] = lanczosInterpolate(f0, d, itpReal) + im*lanczosInterpolate(f0, d, itpImag)
	end

	arr = Float32[fLog 20*log10.(abs.(logData)) angle.(logData)/pi*180]

	return arr
end

function myCSVwrite(arr, filename)
	file = open(filename, "w")
	write(file, "f in Hz; Level in dB; Phase in deg\n")
	writedlm(file, arr, ';')
	close(file)
end

inNames = filter(x -> splitext(x)[2]==".mat", readdir())
for i=1:length(inNames)
	name = inNames[i]
	vars = matread(name)
	myCSVwrite(logResample(vars["data"]|>vec, vars["samplerate"]), string(name[10:end-3], "csv"))
end


text:
function findlocalmaxima(signal)
    # Finds the local maxima of a discrete signal by comparing the
    # value at each point to the neighbouring points.
    #
    # Input:
    # signal: denotes the discrete signal from which we seek the local maxima
    #
    # Output:
    # inds: indices of the local maxima, i.e., signal[inds] are the local maxima

   inds = Int[]
   if length(signal)>1
       if signal[end]<signal[1]>signal[2]
           push!(inds,1)
       end
       i = 2
       while i < length(signal)
           if signal[i-1]<signal[i]>signal[i+1]
               push!(inds,i)
               i += 1
           end
           i += 1
       end
       if signal[1]<signal[end]>signal[end-1]
           push!(inds,length(signal))
       end
   end
   inds
 end


text:
using FileIO
include("julia2ims.jl")
field = ["s1-1" "s11-3" "s12-1" "s23-2" "s27-3" "s29-1" "s32-2" "s34-4" "s36-1" "s3-1" "s11-2" "s23-3" "s27-2" "s32-1" "s33-3" ]
data_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_hamamatsu-60X-TIRF/20200319";
output_dir = "$(data_dir)_visualization"

try
    mkdir("$output_dir")
catch
end

function max_projection(cell)
    println("Processing $cell")
    #raw_img = load(File(format"TIFF", "$data_dir/HE7-11-1-80uw-const_3_$pos.ome.tiff"));
    raw_img = loadims2("$(data_dir)_result/output/d19$cell.ims");
    img_size = size(raw_img);
    img_max = zeros(N0f16, img_size[1], img_size[2], Int(img_size[3]/20));
    for i in 1:Int(img_size[3]/20)
        img_max[:, :, i] = maximum(raw_img[:, :, (i-1)*20+1:i*20], dims=3);
    end
    save("$output_dir/$cell.tiff", img_max);
end

for cell in field
    max_projection(cell)
end


text:
# A Fast Algorithm for Approximate Quantiles in High Speed Data Streams
# Qi Zhang and Wei Wang

struct DataItem{T}
    value::T
    rmax::Int64
end

const QuantileSummary{T} = Array{DataItem{T},1}

function init(values::Array{T,1}) where {T}
    sort!(values)
    summaries = QuantileSummary{T}()
    for i = 1:length(values)
        item = DataItem{T}(values[i], i)
        push!(summaries, item)
    end
    summaries
end

function Base.merge(x::QuantileSummary{T}, y::QuantileSummary{T})::QuantileSummary{T} where {T}

    if length(x) == 0
        return y
    end
    if length(y) == 0
        return x
    end

    idx = 1;
    list = QuantileSummary{T}(undef, length(x) + length(y))

    xi = 1
    yi = 1

    xprevious = nothing
    yprevious = nothing

    while xi <= length(x) || yi <= length(y)
        xitem = xi <= length(x) ? x[xi] : nothing
        yitem = yi <= length(y) ? y[yi] : nothing

        # find the insertion point

        if isnothing(xitem)
            xr = yitem
            ys = xprevious
            yt = nothing

            yprevious = yitem
            yi += 1

        elseif isnothing(yitem)
            xr = xitem
            ys = yprevious
            yt = nothing

            xprevious = xitem
            xi += 1
        else
            if xitem.value < yitem.value
                xr = xitem
                ys = yprevious
                yt = yitem

                xprevious = xitem
                xi += 1
            else

                xr = yitem
                ys = xprevious
                yt = xitem

                yprevious = yitem
                yi += 1

            end
        end

        if isnothing(yt)
            rmax = xr.rmax + ys.rmax
        else
            rmax = xr.rmax + yt.rmax - 1
        end

        list[idx] = DataItem{T}(xr.value, rmax);
        idx += 1;
    end
    list
end

function compress(list::QuantileSummary{T}, sz::Int)::QuantileSummary{T} where {T}

    if length(list) <= sz
        return list
    end

    count = list[end].rmax
    stepsize = count / sz

    newlist = QuantileSummary{T}(undef, sz+1)

    newlist[1] = list[1]
    i = 1
    next = 2;

    for rank = stepsize:stepsize:count
        while list[i].rmax < rank
            i += 1
        end

        newitem = DataItem{T}(
            list[i].value,
            list[i].rmax,
        )

        newlist[next] = newitem;
        next += 1
        i += 1
    end
    newlist
end

function percentile(qs::QuantileSummary, value)
    idx = searchsortedfirst(qs, DataItem(value, 0), by = x -> x.value)
    if idx > length(qs)
        idx = length(qs)
    end

    qs[idx].rmax / Float64(qs[end].rmax)
end

function qindex(qs::QuantileSummary, percentile::Number)
    idx = Int(round(percentile * length(qs))) + 1 # rank is one based
    if idx > length(qs)
        idx = length(qs)
    end
    idx
end

function qvalue(qs::QuantileSummary, percentile::Number)
    qs[qindex(qs, percentile)].value
end

mutable struct FixedSizeSummaryBuilder{T}
    values::Array{T,1}
    levels::Array{QuantileSummary{T},1}

    maxsize::Int64

    FixedSizeSummaryBuilder{T}(eps::Float64, N::Int) where {T} = begin
        maxsize = Int64(floor(log2(eps * N) / eps))
        if (maxsize < 1)
            error("eps $(eps) with N = $N leads to illegal size $(maxsize) ")
        end
        new{T}(Array{T, 1}(), Array{QuantileSummary{T}, 1}(), maxsize)
    end
end

function _fit!(fs::FixedSizeSummaryBuilder{T}, value) where {T}
    push!(fs.values, value)
    if length(fs.values) == fs.maxsize
        summary = init(fs.values)
        inserted = false
        for i = 1:length(fs.levels)
            other = fs.levels[i]
            if isempty(other)
                fs.levels[i] = summary
                inserted = true
                break
            else
                summary = merge(summary, other)
                summary = compress(summary, fs.maxsize)
                fs.levels[i] = QuantileSummary{T}()
            end
        end
        if (!inserted)
            push!(fs.levels, summary)
        end
        empty!(fs.values)
    end
end

function summarize(fs::FixedSizeSummaryBuilder)::QuantileSummary
    summary = init(fs.values)
    for qs in fs.levels
        summary = merge(summary, qs)
    end
    compress(summary, fs.maxsize)
end

function Base.merge(fs::FixedSizeSummaryBuilder)::QuantileSummary
    summary = init(fs.values)
    for qs in fs.levels
        summary = merge(summary, qs)
    end
    summary
end

mutable struct QuantileStat{T} <: OnlineStat{T}
    eps::Float64

    current::Int64
    summary::QuantileSummary{T}
    workingset::FixedSizeSummaryBuilder{T}

    total::Int64
    countSoFar::Int64
    n::Int64

    cachedSummary::Union{Nothing, QuantileSummary{T}}

    QuantileStat{T}(eps::Float64 = 0.01) where {T} = begin
        current = 3
        total = Int64(floor(((1 << current) / eps)))
        n = 0
        workingset = FixedSizeSummaryBuilder{T}(eps / 2, total)
        new{T}(
            eps,
            current,
            QuantileSummary{T}(),
            workingset,
            total,
            n,
            0,
            nothing,
        )
    end
    QuantileStat(eps::Float64 = 0.01) = QuantileStat{Number}(eps)
end

function _fit!(qb::QuantileStat{T}, value) where {T}
    qb.cachedSummary = nothing

    _fit!(qb.workingset, T(value))
    qb.n += 1
    qb.countSoFar += 1
    if qb.countSoFar == qb.total
        qb.current += 1
        qb.total = Int64(floor((1 << qb.current) / qb.eps))
        qb.countSoFar = 1
        qb.summary = merge(
            qb.summary,
            compress(merge(qb.workingset), Int64(floor(2 / qb.eps))),
        )
        qb.workingset = FixedSizeSummaryBuilder{T}(qb.eps / 2, qb.total)
    end
end

targetsize(qb::QuantileStat) = Int64(floor(2/qb.eps))

function _merge!(qb1::QuantileStat{T}, qb2::QuantileStat{T}) where {T}

    # merge current summaries
    qb1.summary = merge(qb1.summary, qb2.summary)

    # merge the current working set of qb2
    qb1.summary = merge(qb1.summary, compress(merge(qb2.workingset), targetsize(qb1)))

    # Make sure nobs is still correct
    qb1.n += qb2.n

    # we keep qb1 in the same state as it was wrt exponentially increasing workingsets

    qb1
end

function summarize(qb::QuantileStat)::QuantileSummary
    if isnothing(qb.cachedSummary)
        qb.cachedSummary = merge(
            qb.summary,
            compress(merge(qb.workingset), Int64(floor(2 / qb.eps))),
        )

        # do a final compression to create a smaller summary that is still an ϵ-approximate summary
        qb.cachedSummary = compress(qb.cachedSummary, Int64(floor(2/qb.eps)))
    end
    qb.cachedSummary
end

ndata(qs::QuantileSummary) = length(qs)
ndata(qb::QuantileStat) = ndata(qb.summary) + ndata(qb.workingset)
ndata(fs::FixedSizeSummaryBuilder) = length(fs.values) + sum(ndata(qs) for qs in fs.levels)

value(qb::QuantileStat) = summarize(qb)

function qvalue(qb::QuantileStat, percentile::Number)
    qvalue(summarize(qb), percentile)
end

function qindex(qb::QuantileStat, percentile::Number)
    qindex(summarize(qb), percentile)
end

function percentile(qs::QuantileStat, value)
    percentile(summarize(qs), value)
end

function median(qb::QuantileStat)
    qvalue(qb, 0.5)
end


text:
#===============================================================================
            Some simple online statistics for the Markov chain
===============================================================================#

"""
    mutable struct GenericChainStats{T} <: ChainStats
        rolling_ar::Vector{Vector{Float64}}
        roll_window::Int64
        mean::Vector{T}
        cov::Matrix{T}
        N::Int64
    end

Simple online statistics for the Markov chain.
"""
mutable struct GenericChainStats{T} <: ChainStats
    rolling_ar::Vector{Vector{Float64}}
    roll_window::Int64
    mean::Vector{T}
    cov::Matrix{T}
    N::Int64

    function GenericChainStats(
            state::Vector{T},
            updt_len,
            num_mcmc_steps,
            roll_window=100
        ) where T
        new{T}(
            [zeros(Float64, updt_len) for _=1:num_mcmc_steps],
            roll_window,
            zero(state),
            zeros(T, (length(state), length(state))),
            1
        )
    end

end

function update_stats!(
        global_ws::GlobalWorkspace,
        local_ws::LocalWorkspace,
        step
    )
    cs = global_ws.sub_ws.stats
    θ = global_ws.sub_ws.state
    old_sum_sq = (cs.N-1)/cs.N * cs.cov + cs.mean * cs.mean'
    cs.mean .= cs.mean .* (cs.N/(cs.N+1)) .+ θ ./ (cs.N+1)
    new_sum_sq = old_sum_sq + (θ * θ')/cs.N
    cs.cov .= new_sum_sq - (cs.N+1)/cs.N*(cs.mean * cs.mean')

    accepted_now = local_ws.acceptance_history[step.mcmciter]
    ra_prev = cs.rolling_ar[max(1,step.mcmciter-1)][step.pidx]
    accepted_outside_window = (
        step.mcmciter > cs.roll_window ?
        local_ws.acceptance_history[step.mcmciter-cs.roll_window] :
        false
    )

    cs.rolling_ar[step.mcmciter][step.pidx] = (
        ra_prev * cs.roll_window
        + (accepted_now - accepted_outside_window)
    )  / min(cs.roll_window, cs.N)
    cs.N += 1
end


text:
# this example shows subplotting and
# use of layout.
using Plots
using Random
y = randn(100,4)
f = scatter(x, y, layout = (4, 1))
display(f)


text:
using Distributed

machines = readlines(ENV["PBS_NODEFILE"])
jobid = ENV["PBS_JOBID"]
addprocs(machines; exename=`/opt/pbs/bin/pbs_attach -j $jobid $(Sys.BINDIR)/julia`)

include(ARGS[2])


text:
############################################################################################
"""
$(TYPEDEF)
Abstract super type for AbstractDifferentiableObjective results.
"""
abstract type ℓObjectiveResult <: BaytesCore.AbstractResult end

"""
$(TYPEDEF)
Stores result for log density and parameter for 'ℓobjective' evaluation at 'parameter'.

# Fields
$(TYPEDFIELDS)
"""
struct ℓDensityResult{T,S} <: ℓObjectiveResult
    "Parameter in unconstrained space."
    θᵤ::T
    "Log density at θᵤ."
    ℓθᵤ::S
    function ℓDensityResult(θᵤ::AbstractVector{S}, ℓθᵤ::S) where {S<:Real}
        return new{typeof(θᵤ),S}(θᵤ, ℓθᵤ)
    end
end
function ℓDensityResult(objective::Objective, θᵤ::AbstractVector{T}) where {T<:Real}
    return ℓDensityResult(θᵤ, T(objective(θᵤ)))
end
function ℓDensityResult(objective::Objective)
    return ℓDensityResult(objective, unconstrain_flatten(objective.model, objective.tagged))
end
ℓDensityResult(diff::DiffObjective) = ℓDensityResult(diff.objective)

"""
$(TYPEDEF)
Stores result for log density, gradient, and parameter for 'ℓobjective' evaluation at 'parameter'.

# Fields
$(TYPEDFIELDS)
"""
struct ℓGradientResult{T,S,G} <: ℓObjectiveResult
    "Parameter in unconstrained space."
    θᵤ::T
    "Log density at θᵤ."
    ℓθᵤ::S
    "Gradient of log density at θᵤ."
    ∇ℓθᵤ::G
    function ℓGradientResult(
        θᵤ::T, ℓθᵤ::S, ∇ℓθᵤ::G
    ) where {T<:AbstractVector,S<:Real,G<:AbstractVector}
        @argcheck length(θᵤ) == length(∇ℓθᵤ)
        return new{T,S,G}(θᵤ, ℓθᵤ, ∇ℓθᵤ)
    end
end
function log_density(
    objective::Objective,
    θᵤ::AbstractVector{T}=unconstrain_flatten(objective.model, objective.tagged),
) where {T<:Real}
    ℓθᵤ = objective(θᵤ)
    if isfinite(ℓθᵤ)
        ℓDensityResult(θᵤ, T(ℓθᵤ))
    else
        ℓDensityResult(θᵤ, T(-Inf))
    end
end
function log_density_and_gradient(
    objective::Objective,
    tune::AbstractDifferentiableTune,
    θᵤ::AbstractVector{T}=unconstrain_flatten(objective.model, objective.tagged),
) where {T<:Real}
    ℓθᵤ, ∇ℓθᵤ = _log_density_and_gradient(objective, tune, θᵤ)
    if isfinite(ℓθᵤ)
        ℓGradientResult(θᵤ, ℓθᵤ, ∇ℓθᵤ)
    else
        #!NOTE: second θᵤ used just as a placeholder
        ℓGradientResult(θᵤ, oftype(ℓθᵤ, -Inf), θᵤ)
    end
end

function log_density(
    diff::DiffObjective,
    θᵤ::AbstractVector{T}=unconstrain_flatten(diff.objective.model, diff.objective.tagged),
) where {T<:Real}
    ℓθᵤ = _log_density(diff, θᵤ)
    if isfinite(ℓθᵤ)
        ℓDensityResult(θᵤ, T(ℓθᵤ))
    else
        ℓDensityResult(θᵤ, T(-Inf)) #oftype(ℓθᵤ, -Inf))
    end
end
function log_density_and_gradient(
    diff::DiffObjective,
    θᵤ::AbstractVector{T}=unconstrain_flatten(diff.objective.model, diff.objective.tagged),
) where {T<:Real}
    ℓθᵤ, ∇ℓθᵤ = _log_density_and_gradient(diff, θᵤ)
    if isfinite(ℓθᵤ)
        ℓGradientResult(θᵤ, ℓθᵤ, ∇ℓθᵤ)
    else
        #!NOTE: second θᵤ used just as a placeholder
        ℓGradientResult(θᵤ, oftype(ℓθᵤ, -Inf), θᵤ)
    end
end
#=
struct ℓHessianResult{T, S} <: ℓObjectiveResult
end
function log_density_and_hessian(objective::L, θᵤ::AbstractVector{T}) where {L<:AbstractDifferentiableObjective, T<:Real}
end
=#

############################################################################################
function checkfinite(result::T) where {T<:ℓObjectiveResult}
    return isfinite(result.ℓθᵤ) && _checkfinite(result.θᵤ) ? true : false
end
function checkfinite(
    result₀::T, result::T, min_Δ::Float64=min_Δ
) where {T<:ℓObjectiveResult}
    checkfinite(result) && ((result.ℓθᵤ - result₀.ℓθᵤ) > min_Δ) || return false
    return true
end
function checkfinite(
    ℓθ₀::R, ℓθ::R, result::T, min_Δ::Float64=min_Δ
) where {R<:Real,T<:ℓObjectiveResult}
    checkfinite(result) && ((ℓθ - ℓθ₀) > min_Δ) || return false
    return true
end

############################################################################################
# Export
export
    ℓObjectiveResult,
    ℓGradientResult,
    ℓDensityResult,
    log_density,
    log_density_and_gradient#,
#log_density_and_hessian


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %% [markdown]
# * https://github.com/JuliaMath/Combinatorics.jl/blob/master/src/permutations.jl#L47
# * https://discourse.julialang.org/t/is-there-a-function-behaving-the-same-as-next-permutation-does-in-c/63451/13

# %%
module O

"""
zero allocation version of AquaIndigo's code
"""
function next_perm_old!(v::AbstractVector)
    length(v) ≤ 1 && return false
    k = findlast(isless(v[i], v[i+1]) for i in firstindex(v):lastindex(v)-1)
    isnothing(k) && (reverse!(v); return false)
    i = k + firstindex(v) - 1
    j = findlast(isless(v[i], v[j]) for j in i+1:lastindex(v)) + i
    v[i], v[j] = v[j], v[i]
    reverse!(v, i + 1)
    return true
end

"""
    next_perm!(v::AbstractVector)

changes `v` into the next permutation of `v` under the lexicographic order and returns `false` if it is the last permutation and `true` otherwise.

This is a variant of `Combinatorics.nextpermutation(m, t, state)`.

Examples

```
julia> v = [1, 2, 3]; println(v); while O.next_perm!(v) println(v) end
[1, 2, 3]
[1, 3, 2]
[2, 1, 3]
[2, 3, 1]
[3, 1, 2]
[3, 2, 1]
```

```
julia> v = [1, 2, 2]; println(v); while O.next_perm!(v) println(v) end
[1, 2, 2]
[2, 1, 2]
[2, 2, 1]
```
"""
function next_perm!(v::AbstractVector)
    length(v) ≤ 1 && return false
    i = lastindex(v) - 1
    @inbounds while i ≥ firstindex(v) && !isless(v[i], v[i+1]) i -= 1 end
    i < firstindex(v) && (reverse!(v); return false)
    j = lastindex(v)
    @inbounds while j > i && !isless(v[i], v[j]) j -= 1 end
    @inbounds v[i], v[j] = v[j], v[i]
    reverse!(v, i + 1)
    return true
end

function collect_perm!(v::AbstractVector)
    a = [copy(v)]
    while next_perm!(v) push!(a, copy(v)) end
    a
end

"""
    Perm(v::AbstractVector)

constructs the lexicographic order iterator of the all permutations of `v`.

Examples

```
julia> for s in O.Perm([1, 2, 3]) println(s) end
[1, 2, 3]
[1, 3, 2]
[2, 1, 3]
[2, 3, 1]
[3, 1, 2]
[3, 2, 1]
```

```
julia> for s in O.Perm([1, 2, 2]) println(s) end
[1, 2, 2]
[2, 1, 2]
[2, 2, 1]
```
"""
struct Perm{T<:AbstractVector}
    v::T
    Perm(v::T) where T<: AbstractVector = new{T}(sort(v))
end

function Base.length(p::Perm)
    v = p.v
    N = factorial(length(v))
    i = firstindex(v)
    for j in eachindex(v)
        if isless(v[i], v[j])
            N = N ÷ factorial(j-i)
            i = j
        end
    end
    N = N ÷ factorial(lastindex(v)-i+1)
    N
end

Base.eltype(p::Perm{T}) where T<:AbstractVector = T

function Base.iterate(p::Perm)
    s = copy(p.v)
    t = copy(s)
    next_perm!(t) || return (s, nothing)
    u = copy(t)
    next_perm!(u) || return (s, (t, nothing))
    return (s, (t, u))
end

function Base.iterate(p::Perm, state)
    isnothing(state) && return nothing
    s, t = state
    isnothing(t) && return (s, nothing)
    u = copy(t)
    next_perm!(u) || return (s, (t, nothing))
    return (s, (t, u))
end

#function Base.iterate(p::Perm, s = copy(p.v))
#    isnothing(s) && return nothing
#    next_perm!(s) || return (s, nothing)
#    return (copy(s), s)
#end

end

# %%
@doc O.next_perm!

# %%
@doc O.Perm

# %%
v = [1, 2, 3]; println(v); while O.next_perm!(v) println(v) end

# %%
v = [1, 2, 2]; println(v); while O.next_perm!(v) println(v) end

# %%
for s in O.Perm([1, 2, 3]) println(s) end

# %%
for s in O.Perm([1, 2, 2]) println(s) end

# %%
v = sort([NaN, NaN, 2, 1])
@show v
while O.next_perm!(v)
    @show v
end

# %%
p = O.Perm([NaN, NaN, 2, 1])
for s in p
    @show s
end

# %%
using BenchmarkTools

# %%
p = O.Perm([NaN, NaN, 2, 1])
_, state = iterate(p)

# %%
iterate(p, state)

# %%
p = O.Perm([NaN, NaN, 2, 1])
@btime collect($p)

# %%
v = sort([NaN, NaN, 2, 1])
@btime O.collect_perm!($v)

# %%
test_next_perm!(v) = (c = 1; while O.next_perm!(v) c += 1 end; c)
@show s = [1, 1, 2, 2, 2, 3, 4, 5, 6, 6]
@btime test_next_perm!($s)

# %%
test_next_perm_old!(v) = (c = 1; while O.next_perm_old!(v) c += 1 end; c)
@show s = [1, 1, 2, 2, 2, 3, 4, 5, 6, 6]
@btime test_next_perm_old!($s)

# %%
c = @btime collect($(O.Perm(s)))
length(c)

# %%
@code_warntype O.next_perm!(s)

# %%
using OffsetArrays
@show t = OffsetArray([1, 1, 2, 2, 2, 3, 4, 5, 6, 6], -5:4)
@btime collect($(O.Perm(t)))

# %%
@btime O.collect_perm!($(sort(t)))

# %%
@code_warntype O.next_perm!(t)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
mealCost = parse(Float64, readline())
tipPercentage = parse(Float64, readline())
taxPercentage = parse(Float64, readline())

totalCost = 0.0
totalCost += mealCost
totalCost += tipPercentage * 0.01 * mealCost
totalCost += taxPercentage * 0.01 * mealCost

println(convert(Int64, round(totalCost) ) )

max_stars_repo_path: Hackerrank/30 Days of Code/Julia/day 02.jl

max_stars_repo_name: Next-Gen-UI/Code-Dynamics

max_issues_repo_name: Next-Gen-UI/Code-Dynamics

max_forks_repo_name: Next-Gen-UI/Code-Dynamics


text:
# Consider all integer combinations of ab for 2 ≤ a ≤ 5 and 2 ≤ b ≤ 5:
#
# 2^2=4,  2^3=8,   2^4=16,  2^5=32
# 3^2=9,  3^3=27,  3^4=81,  3^5=243
# 4^2=16, 4^3=64,  4^4=256, 4^5=1024
# 5^2=25, 5^3=125, 5^4=625, 5^5=3125
#
# If they are then placed in numerical order, with any repeats removed, we get
# the following sequence of 15 distinct terms:
#
# 4, 8, 9, 16, 25, 27, 32, 64, 81, 125, 243, 256, 625, 1024, 3125
#
# How many distinct terms are in the sequence generated by ab for 2 ≤ a ≤ 100
# and 2 ≤ b ≤ 100?

using ProjectEulerSolutions

# Just try all combinations and store unique ones using a set.  BigInt must be
# used because 100^100 is larger than 2^64.  Can easily make it a bit faster
# using some additinoal logic to prevent trying multiples of powers.  Might be
# possible to do it without even performing powers, just looking at counts of
# prime factors, foregoing the need for BigInts, but that would probably be
# slower.
function p029solution_loop(amin::Integer=2, amax::Integer=5,
                      bmin::Integer=2, bmax::Integer=5)::Integer
    powers = Set{BigInt}()
    for a in amin:amax
        for b in bmin:bmax
            push!(powers, BigInt(a)^b)
        end
    end
    return length(powers)
end

function p029solution_listcomprehension(amin::Integer=2, amax::Integer=5,
                                        bmin::Integer=2, bmax::Integer=5)::Integer
    return length(Set{BigInt}([BigInt(a)^b for a in amin:amax, b in bmin:bmax]))
end

p029 = Problems.Problem(Dict("Explicit loop" => p029solution_loop,
                             "List comprehension" => p029solution_listcomprehension))

Problems.benchmark(p029, 2, 100, 2, 100)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
using Khepri
#@static if VERSION < v"0.7.0-DEV.2005"
#    using Base.Test
#else
#    using Test
#end

# write your own tests here
#@test 1 == 2


text:
using ApproxFun, SpecialFunctions, LinearAlgebra, Test
    import DomainSets
    import DomainSets: UnionDomain
    import ApproxFun: ChebyshevDirichlet, Ultraspherical, PiecewiseSegment, ContinuousSpace, space, SpaceOperator,
                        testspace, testbandedoperator, testraggedbelowoperator, testcalculus, testtransforms,
                        testfunctional

@testset "Spaces" begin
    @testset "ChebyshevDirichlet" begin
        testtransforms(ChebyshevDirichlet{1,1}())

        @test Fun(exp,ChebyshevDirichlet{1,1})(.1) ≈ exp(.1)
        @test Fun(Fun(exp,ChebyshevDirichlet{1,1}),Ultraspherical(1))(.1) ≈ exp(.1)

        @test domain(ChebyshevDirichlet{1,1}()) === ChebyshevInterval()
        @test Fun(Fun(exp,ChebyshevDirichlet{1,1}(Segment(-1,1))),Ultraspherical(1))(.1) ≈ exp(.1)

        d = ChebyshevInterval()
        sp = ChebyshevDirichlet{1,1}(d)
        B=Dirichlet(sp)
        D=Derivative(d)
        L=D^2+I
        @test B[1:2,1:4] ≈ [1 -1 0 0; 1 1 0 0]

        @test csc(2)sin(1 - 0.1)  ≈ ([Dirichlet(d);L]\[[1.,0.],0.])(0.1)
        @test csc(2)sin(1 - 0.1)  ≈ ([B;L]\[[1.,0.],0.])(0.1)

        @test norm(([B;L]\[[1.,0],0])-([Dirichlet(d);L]\[[1.,0],0])) <10eps()
    end

    @testset "PiecewiseSpace" begin
        x = Fun(identity,UnionDomain(-1..0, 0..1))
        sp=space(x)
        testtransforms(sp;minpoints=2)

        D = Derivative(sp)
        testbandedoperator(D)

        B = [Dirichlet(sp); continuity(sp,0:1)]
        u = [B; D^2] \ [[1,0],zeros(2),0];
        u2 = [Dirichlet();Derivative(Chebyshev())^2] \ [[1.,0],0]
        @test u(0.) ≈ u2(0.)

        x = Fun(identity, UnionDomain(-10..0, 0..1, 1..15))
        sp=space(x)
        D=Derivative(sp)
        B=Dirichlet(sp)

        u=[B;
            continuity(sp,0:1);
            D^2-x]\[[airyai(-10.),0],zeros(4),0];

        @test u(0.) ≈ airyai(0.)

        s=Fun(sin,-2..2)|>abs
        c=Fun(cos,-2..2)|>abs
        sc=Fun(x -> abs(sin(x))+abs(cos(x)), ContinuousSpace(PiecewiseSegment([-2,(-π/2),0,π/2,2])))
        @test norm(sc-(c+s))<100eps()
    end

    @testset "max/min creates breakpoints" begin
        x = Fun()
        g = 4*(x-0.2)
        f = max(-1,g)
        f2 = min(f,1)

        @test norm(max(x,x)-x)<100eps()
        @test norm(min(x,x)-x)<100eps()

        f3=Fun(x->x<-0.05 ? -1.0 : (x<0.45 ? 4*(x-.2) : 1), (-1..1) \ [-0.05,0.45])
        @test norm(f2.(range(-1,stop=1,length=10))-f3.(range(-1,stop=1,length=10))) < 4eps()
    end

    @testset "Complex piecewise" begin
        x = Fun(identity, Segment(im,0) ∪ Segment(0,1))
        @test x(0.5) ≈ 0.5
        @test x(0.5im) ≈ 0.5im

        @test Fun(Fun(1.0),space(x))(0.5) == 1.0
        @test Fun(Fun(1.0),space(x))(0.5im) == 1.0

        @test (x+1)(0.5) ≈ 1.5
        @test (x-1)(0.5) ≈ -0.5
        @test (1-x)(0.5) ≈ 0.5

        @test sqrt(1-x)(0.2im) ≈ sqrt(1-0.2im)
        @test sqrt(1-x)(0.2) ≈ sqrt(1-0.2)

        w=2/(sqrt(1-x)*sqrt(1+im*x))
        @test w(0.2im) ≈ 2/(sqrt(1-0.2im)*sqrt(1+im*(0.2im)))
        @test w(0.2) ≈ 2/(sqrt(1-0.2)*sqrt(1+im*(0.2)))
    end

    @testset "Equivalent spaces" begin
        @test norm(Fun(cos,Chebyshev)-Fun(cos,Jacobi(-0.5,-0.5)))<100eps()
        @test norm(Fun(cos,Chebyshev)-Fun(cos,JacobiWeight(0,0)))<100eps()
        @test norm(Fun(cos,Jacobi(-0.5,-0.5))-Fun(cos,JacobiWeight(0,0))) < 100eps()
        @test norm(Fun(cos,Chebyshev)-Fun(cos,JacobiWeight(0,0,Jacobi(-0.5,-0.5))))<100eps()
        @test norm(Fun(cos,Jacobi(-0.5,-0.5))-Fun(cos,JacobiWeight(0,0,Jacobi(-0.5,-0.5))))<100eps()
    end

    @testset "ContinuousSpace" begin
        d = PiecewiseSegment(1.,2.,3.,4.)
        S = ContinuousSpace(d)
        testtransforms(S; minpoints=3, invertibletransform=false)

        D = Derivative(S)
        testbandedoperator(D)

        A = [ldirichlet(S); D-I]
        testraggedbelowoperator(A)
        u = [ldirichlet(S); D-I] \ [exp(1.); 0]

        @test u(1.1) ≈ exp(1.1)
        @test u(3.4) ≈ exp(3.4)
        @test last(u) ≈ exp(4)

        d=PiecewiseSegment(0,1.,1+im,im,0.)
        s=Space(d)

        # narrow down bug
        @test s isa ContinuousSpace
        @test ApproxFun.rangetype(s) == Float64
        cs=ApproxFun.canonicalspace(s)
        @test ApproxFun.rangetype(cs) == Float64

        @test conversion_type(s,cs) == s

        @test Fun(exp,d)(.1) ≈ exp(.1)
    end

    @testset "Triple SumSpace" begin
        x=Fun()
        w=log(1-x)+sqrt(1-x^2)
        f=w+x
        @test f(0.1) ≈ (w(0.1)+0.1)
        @test (w+1)(0.1) ≈ (w(0.1)+1)
        @test (w+x+1)(0.1) ≈ (w(0.1)+1.1)
        @test ((w+x)+1)(0.1) ≈ (w(0.1)+1.1)
    end

    @testset "SumSpace bug" begin
        dsp=JacobiWeight(1.,0.,Jacobi(1.,0.,0..1))⊕JacobiWeight(0.5,0.,Jacobi(0.5,-0.5,0..1))
        rsp=Legendre(0..1)⊕JacobiWeight(0.5,0.,Jacobi(0.5,0.5,0..1))

        C=Conversion(dsp,rsp)

        f=Fun(dsp,[1.,2.,3.,4.,5.])
        @test f(0.1) ≈ (C*f)(0.1)
    end

    @testset "Piecewise + Constant" begin
        Γ=Circle() ∪ Circle(0.0,0.4)
        o=ones(Γ)
        @test o(1.) ≈ 1.0
        @test o(0.4) ≈ 1.0

        G=Fun(z->in(z,component(Γ,2)) ? [1 0; -1/z 1] : [z 0; 0 1/z],Γ)
        @test (G-I)(exp(0.1im)) ≈ (G(exp(0.1im))-I)
    end

    @testset "Previoius segfault" begin
        x=Fun(identity,-1..1)
        f=x+sin(2x)*sqrt(1-x^2)
        @test f(0.1) ≈ 0.1+sin(2*0.1)*sqrt(1-0.1^2)
    end

    @testset "Multiple piecewisespace" begin
        x=Fun(identity,-3 .. -2)+Fun(identity,2..3)
        w=sqrt(9-x^2)
        f=w+Fun()
        @test (f+w)(2.5) ≈ 2w(2.5)
        @test (f+w)(.5) ≈ f(.5)
    end

    @testset "Jacobi recurrence bug" begin
        S=Jacobi(0.5,-0.5)
        f=Fun(exp,S)
        @test f(0.1) ≈ exp(0.1)
    end

    @testset "Jacobi evaluation bug" begin
        S=Jacobi(0.5,-0.5)
        f=Fun(exp,S)
        B = ldirichlet(S)
        testfunctional(B)
        @test B*f ≈ exp(-1)
    end

    @testset "cancellation conversion" begin
        x=Fun(0..1)
        f=exp(x)-1
        Fun(f,JacobiWeight(1.,0.,0..1))
    end

    @testset "Hermite" begin
        f=Fun(x->x+x^2,Hermite())
        @test f(1.) ≈ 2.
        @test values(f) ≈ points(f)+points(f).^2

        D = Derivative(Hermite())
        testbandedoperator(D)

        g = D*f
        @test g(1.) ≈ 3.
    end

    @testset "Arc exp" begin
        z=Fun(identity,Arc(0.,.1,0.,π/2))
        @test exp(z)(0.1exp(0.2im)) ≈ exp(0.1exp(0.2im))
    end

    @testset "Extending function" begin
        Γ = Segment(-im,1.0-im) ∪ Curve(Fun(x->exp(0.8im)*(x+x^2-1+im*(x-4x^3+x^4)/6))) ∪ Circle(2.0,0.2)

        @test isempty(component(Γ,1)\component(Γ,1))
        @test Γ \ component(Γ,1) == component(Γ,2) ∪ component(Γ,3)

        @test norm(Fun(ones(component(Γ,1)),Γ) - Fun(x->x ∈ component(Γ,1) ? 1.0 : 0.0,Γ)) == 0
    end

    @testset "Line" begin
        f=Fun(z->2exp(z^2),PeriodicLine(0.,π/2))
        @test f(1.1im) ≈ 2exp(-1.1^2)

        f=Fun(z->2exp(z^2),Line(0.,π/2))
        @test f(1.1im) ≈ 2exp(-1.1^2)
    end

    @testset "Exp Curve" begin
        a = 1+10*im; b = 2-6*im
        d = Curve(Fun(x->1+a*x+x^2+b*x^3))

        x=Fun(d)

        @test exp(x)(1+a*0.1+0.1^2+b*0.1^3) ≈ exp(1+a*0.1+0.1^2+b*0.1^3)
    end

    @testset "ChebyshevDirichlet multiplication" begin
        S=ChebyshevDirichlet()
        x=Fun()
        @test norm((ApproxFun.Recurrence(S)*Fun(exp,S)-Fun(x->x*exp(x),S)).coefficients) < 100eps()
        @test norm((x*Fun(exp,S)-Fun(x->x*exp(x),S)).coefficients) < 100eps()

        @test Integral(S,1)*Fun(S,[1.,2.,3.]) ≈ integrate(Fun(Fun(S,[1.,2.,3.]),Chebyshev()))
    end

    @testset "QuotientSpace" begin
        for (bcs,ret) in ((Dirichlet(Chebyshev()),[1 -1 0 0 0;1 1 0 0 0]),
                          (Neumann(Chebyshev()),[0 1 -4 0 0;0 1 4 0 0]),
                          ([DefiniteIntegral(Chebyshev());SpaceOperator(BasisFunctional(2),Chebyshev(),ConstantSpace())],[2 0 0 0 0;0 1 0 0 0]),
                          (vcat(bvp(Chebyshev(),4)...),[1 -1 1 -1 0;0 1 -4 9 0;1 1 1 1 0;0 1 4 9 0]))
            QS = QuotientSpace(bcs)
            C = Conversion(QS, QS.space)

            norm((bcs*C)[1:size(bcs, 1),1:5] - ret) < 1000eps()
        end
    end

    @testset "Union of ChebyshevDirichlet" begin
        dom = UnionDomain(0..1, 2..3)
        @test components(union(JacobiWeight.(-0.5,-0.5,ChebyshevDirichlet{1,1}.(components(dom)))...)) ==
            (JacobiWeight.(-0.5,-0.5,ChebyshevDirichlet{1,1}.(components(dom)))...,)
    end

    @testset "Ultraspherical special functions" begin
        x = Fun(Ultraspherical(2,0..1))
        sqrt(x)(0.1) ≈ sqrt(0.1)

        f = Fun(x->x*exp(x),Ultraspherical(1,0..1))
        sqrt(f(0.1)) ≈ sqrt(f)(0.1)
    end

    @testset "Fast Ultraspherical*Chebyshev" begin
        f = Fun(exp, Ultraspherical(1))
        g = Fun(exp)

        @test ApproxFun.hasfasttransform(f)
        @test ApproxFun.pointscompatible(f,g)
        @test ApproxFun.hasfasttransformtimes(f,g)
        @test ApproxFun.transformtimes(f,g)(0.1) ≈ exp(0.2)
    end

    @testset "PointSpace" begin
        δ = Fun(ApproxFun.PointSpace([2.0]),[1.0])
        f = Fun(x->cos(50x)) + δ
        g = Fun(x->cos(50x),Ultraspherical(1)) + δ
        @test domain(f) == domain(g)

        @test (f*g)(0.1) ≈ cos(50*0.1)^2
        @test (f*g)(2.0) ≈ 1

        @test exp(f)(0.1) ≈ exp(cos(50*0.1))
        @test exp(f)(2.0) ≈ exp(1)

        @test iszero(sign(-Fun(zero)))
        @test iszero(sign(Fun(zero)))
        @test iszero(abs(-Fun(zero)))
        @test iszero(abs(Fun(zero)))
        @test angle(-Fun(zero)) ≈ Fun(π)
        @test iszero(angle(Fun(zero)))

        @test sign(f)(0.1) ≈ sign(cos(50*0.1))
        @test sign(f)(2.0) ≈ 1
        @test abs(f)(0.1) ≈ abs(cos(50*0.1))
        @test abs(f)(2.0) ≈ 1
        @test angle(f)(0.1) ≈ angle(cos(50*0.1))
        @test angle(f)(2.0) ≈ 0
    end

    @testset "Jump Locations" begin
        x = Fun(UnionDomain(0..1, 2..3))
        @test length(jumplocations(sign(x))) == 0

        x = Fun(Chebyshev(-1..1))
        @test length(jumplocations(x)) == 0
        @test all(jumplocations(sign(x) + sign(x+0.2)) .≈ [-0.2, 0])
    end

    @testset "one for SumSpace" begin
        S = Jacobi(0,1) ⊕ JacobiWeight(1/3,0,Jacobi(1/3,2/3)) ⊕ JacobiWeight(2/3,0,Jacobi(2/3,1/3))
        o = ones(S)
        @test o(0.5) ≈ 1
    end

    @testset "blockbandwidths for FiniteOperator of pointscompatibleace bug" begin
        S = ApproxFun.PointSpace([1.0,2.0])
        @test ApproxFun.blockbandwidths(FiniteOperator([1 2; 3 4],S,S)) == (0,0)
    end

    @testset "SumSpace Conversion" begin
        H = ApproxFun.HeavisideSpace([-1.0,0.0,1.0])
        C = ApproxFun.ContinuousSpace(ApproxFun.PiecewiseSegment([-1.0,0,1]))
        S = H + C
        P = Ultraspherical(1,Segment(-1.0,0.0)) ∪ Ultraspherical(1,Segment(0.0,1.0))
        f = Fun(S, randn(100))
        @test f(0.1) ≈ Fun(f, P)(0.1)

        @test Conversion(S,P)[1,1] == 1.0
        @test Conversion(S,P)[1,2] == 0.5
        @test rangespace(Conversion(S,P)) == Ultraspherical(1,Segment(-1.0,0.0)) ∪ Ultraspherical(1,Segment(0.0,1.0))
        @test domainspace(Conversion(S,P)) ==
            ApproxFun.SumSpace(ApproxFun.HeavisideSpace([-1.0,0.0,1.0]),
                    ApproxFun.ContinuousSpace(ApproxFun.PiecewiseSegment([-1.0,0,1])))

        D = ApproxFun.DiracSpace([-1.0,0.0,1.0])
        S2 = ApproxFun.SumSpace(D , P)
        f = Fun(S, randn(100))
        @test (Conversion(S,S2) * f)(0.1) ≈ f(0.1)
    end

    @testset "Mix Fourier-Chebyshev (#602)" begin
        s = Chebyshev(-π..π)
        a = Fun(t-> 1+sin(cos(2t)), s)
        L = Derivative() + a
        f = Fun(t->exp(sin(10t)), s)
        B = periodic(s,0)
        uChebyshev = [B;L] \ [0.;f]

        s = Fourier(-π..π)
        a = Fun(t-> 1+sin(cos(2t)), s)
        L = Derivative() + a
        f = Fun(t->exp(sin(10t)), s)
        uFourier = L\f

        @test norm(uFourier-uChebyshev) ≤ 100eps()
    end

    @testset "real line integral" begin
        f = x -> x^2/(x^4+1)

        z₁,z₂,z₃,z₄ = exp(im*π/4), exp(3im*π/4), exp(5im*π/4), exp(7im*π/4)

        res₁ = z₁^2 / ((z₁ - z₂)*(z₁ - z₃)*(z₁ - z₄) )
        res₂ = z₂^2 / ((z₂ - z₁)*(z₂ - z₃)*(z₂ - z₄) )

        @test_skip sum(Fun(f, Line())) ≈ 2π*im*(res₁ + res₂)
    end

    @testset "Piecewise Legendre" begin
        f = Fun(x -> exp(-40(x-0.1)^2), Legendre(0..1), 1000)
        @test f(0.2) ≈ exp(-40(0.1)^2)

        f = Fun(x -> exp(-40(x-0.1)^2), Legendre(0..1))
        @test f(0.2) ≈ exp(-40(0.1)^2)

        sp = Legendre(Segment(0 , -1)) ⊕ Legendre(0 .. 1)
        @time f = Fun(x->sign(x)*exp(-40(x-0.1)^2), sp)
        @test f(0.2) ≈ exp(-40(0.1)^2)
        @test f(-0.2) ≈ -exp(-40(0.3)^2)
    end

    @testset "Ambiguous evaluate" begin
        f = Fun(1.0, Chebyshev(NaN..NaN))
        @test f(0.1) ≈ 1.0
    end

    @testset "x + abs(x) (#637)" begin
        x = Fun()
        g = abs(x) + x
        @test ncoefficients(g) == 3
        @test g(0.1) ≈ 0.2
        @test g(-0.2) ≈ 0.0
    end

    @testset "piecewise sample (#635)" begin
        f = abs(Fun(sin, -5..5))
        @test integrate(f)(-4.0) ≈ -(cos(-4.0) - cos(-5.0))
        @test -(cos(-π) - cos(-5.0)) + cos(-3.0) - cos(-π) ≈ integrate(f)(-3.0)
        r = ApproxFun.sample(f,10)
        @test maximum(r) ≤ 5
        @test minimum(r) ≥ -5
    end
end


text:
module OpenTelemetryExporterPrometheus

export PrometheusExporter

using OpenTelemetrySDK
using HTTP
using Sockets

"""
    PrometheusExporter(; host = "127.0.0.1", port = 9966, kw...)

It will setup a http server configured by `host` and `port` on initialization.
Here the extra keyword arguments `kw` will be forwarded to `HTTP.listen`.

## Usage

```julia
r = MetricReader(PrometheusExporter())
```

Note that since `PrometheusExporter` is pull based exporter, which means there's
no need to execute `r()` to update the metrics.
"""
mutable struct PrometheusExporter <: AbstractExporter
    server::Sockets.TCPServer
    provider::Union{MeterProvider,Nothing}
    taskref::Ref{Task}
    function PrometheusExporter(; host = "127.0.0.1", port = 9966, kw...)
        server = Sockets.listen(Sockets.InetAddr(parse(IPAddr, host), port))
        exporter = new(server, nothing, Ref{Task}())
        exporter.taskref[] = @async HTTP.listen(host, port; server = server, kw...) do http
            HTTP.setstatus(http, 200)
            HTTP.setheader(http, "Content-Type" => "text/plain")

            if isnothing(exporter.provider)
                write(http, "MeterProvider is not set yet!!!")
            else
                for ins in keys(exporter.provider.async_instruments)
                    ins()
                end
                text_based_format(http, exporter.provider)
            end
            return
        end
        finalizer(exporter) do x
            shut_down!(x)
        end
        exporter
    end
end

OpenTelemetrySDK.shut_down!(r::PrometheusExporter) = close(r.server)

function (r::MetricReader{<:MeterProvider,<:PrometheusExporter})()
    if isnothing(r.exporter.provider)
        r.exporter.provider = r.provider
    else
        @info "The prometheus exporter is already properly set!"
    end
end

# TODO: support exemplars
function text_based_format(io, provider::MeterProvider)
    for m in metrics(provider)
        write(io, "# HELP $(m.name) $(m.description)\n")
        write(io, "# TYPE $(m.name) $(prometheus_type(m.aggregation))\n")
        for (attrs, point) in m
            if point.value isa OpenTelemetrySDK.HistogramValue
                val = point.value
                for (i, c) in enumerate(Iterators.accumulate(+, val.counts))
                    if length(attrs) > 0
                        if i == length(val.counts)
                            write(io, "$(m.name)_bucket{")
                            # TODO: escape
                            join(io, ("$k=\"$v\"" for (k, v) in pairs(attrs)), ",")
                            write(io, ",le=\"+Inf\"} $c $(point.time_unix_nano ÷ 10^6) \n")

                            write(io, "$(m.name)_count{")
                            # TODO: escape
                            join(io, ("$k=\"$v\"" for (k, v) in pairs(attrs)), ",")
                            write(io, "} $c $(point.time_unix_nano ÷ 10^6) \n")
                        else
                            write(io, "$(m.name)_bucket{")
                            # TODO: escape
                            join(io, ("$k=\"$v\"" for (k, v) in pairs(attrs)), ",")
                            write(
                                io,
                                ",le=\"$(val.boundaries[i])\"} $c $(point.time_unix_nano ÷ 10^6) \n",
                            )
                        end
                    else
                        if i == length(val.counts)
                            write(
                                io,
                                "$(m.name)_bucket{le=\"+Inf\"} $c $(point.time_unix_nano ÷ 10^6) \n",
                            )
                            write(io, "$(m.name)_count $c\n")
                        else
                            write(
                                io,
                                "$(m.name)_bucket{le=\"$(val.boundaries[i])\"} $c $(point.time_unix_nano ÷ 10^6) \n",
                            )
                        end
                    end
                end
                # ???
                if !isnothing(val.sum)
                    if length(attrs) > 0
                        write(io, "$(m.name)_sum{")
                        # TODO: escape
                        join(io, ("$k=\"$v\"" for (k, v) in pairs(attrs)), ",")
                        write(io, "} $(val.sum) $(point.time_unix_nano ÷ 10^6) \n")
                    else
                        write(io, "$(m.name)_sum $(val.sum)\n")
                    end
                end
            else
                write(io, "$(m.name){")
                # TODO: escape
                join(io, ("$k=\"$v\"" for (k, v) in pairs(attrs)), ",")
                write(io, "} $(point.value) $(point.time_unix_nano ÷ 10^6) \n")
            end
        end
        write(io, "\n")
    end
end

prometheus_type(::SumAgg) = "counter"
prometheus_type(::LastValueAgg) = "gauge"
prometheus_type(::HistogramAgg) = "histogram"

end # module


text:
using Base.Threads: @threads

function (nw::Network)(du, u::T, p, t) where {T}
    # for now just one layer
    layer = nw.nl
    dupt = (du, u, p, t)

    # get cache array for accumulator of appropriate type
    _acc = accumulator_cache(layer, T)
    fill!(_acc, zero(eltype(T)))

    # go through all colobatches separatly to avoid writing conflicts to accumulator
    unroll_colorbatches!(nw, layer, layer.colorbatches, dupt, _acc)

    # can be run parallel
    process_batches!(nw, layer, nw.vertexbatches, dupt, _acc)
    return nothing
end

@unroll function unroll_colorbatches!(nw, layer, colorbatches, dupt, _acc)
    @unroll for cbatch in colorbatches
        process_batches!(nw, layer, cbatch.edgebatches, dupt, _acc)
    end
end

@unroll function process_batches!(nw, layer, batches, dupt, acc)
    @unroll for batch in batches
        process_batch!(nw, layer, batch, dupt, acc)
    end
end

function process_batch!(nw::Network{:seq}, layer, batch, dupt, acc)
    for idx in 1:length(batch)
        element_kernel!(layer, batch, dupt, acc, idx)
    end
end

function process_batch!(nw::Network{:threaded}, layer, batch, dupt, acc)
    @threads for idx in 1:length(batch)
        element_kernel!(layer, batch, dupt, acc, idx)
    end
end

@inline function element_kernel!(layer, batch::VertexBatch{F}, dupt, acc, i) where {F}
    @inbounds begin
    du, u, p, t = dupt
    (v_r, acc_r) = vertex_ranges(layer, batch, i)
    pv_r = parameter_range(batch, i)

    vdu = view(du, v_r)
    vu  = view(u, v_r)
    vacc = view(acc, acc_r)
    pv = p===nothing ? nothing : view(p, pv_r)

    batch.fun.f(vdu, vu, vacc, pv, t)
    end
end

@inline function element_kernel!(layer, batch::EdgeBatch{F}, dupt::T, acc, i) where {F<:StaticEdge, T}
    @inbounds begin
    # a cache of size dim per thread
    _, u, p, t = dupt

    # collect all the ranges to index into the data arrays
    (src_r, dst_r, src_acc_r, dst_acc_r, _) = src_dst_ranges(layer, batch, i)
    pe_r = parameter_range(batch, i)

    # create the views into the data & parameters
    vs = view(u, src_r)
    vd = view(u, dst_r)
    pe = p===nothing ? nothing : view(p, pe_r)

    # apply the edge function
    e = batch.fun.f(vs, vd, pe, t)

    apply_accumulation!(coupling(F), layer.accumulator, acc, src_acc_r, dst_acc_r, e)
    end
end

@inline function element_kernel!(layer, batch::EdgeBatch{F}, dupt, acc, i) where {F<:ODEEdge}
    @inbounds begin
    du, u, p, t = dupt
    # for i in 1:length(batch)
    # collect all the ranges to index into the data arrays
    (src_r, dst_r, src_acc_r, dst_acc_r, e_r) = src_dst_ranges(layer, batch, i)
    pe_r = parameter_range(batch, i)

    # create the views into the data & parameters
    de = view(du, e_r)
    e  = view(u, e_r)
    vs = view(u, src_r)
    vd = view(u, dst_r)
    pe = p===nothing ? nothing : view(p, pe_r)

    # apply the edge function
    batch.fun.f(de, e, vs, vd, pe, t)

    apply_accumulation!(coupling(F), layer.accumulator, acc, src_acc_r, dst_acc_r, e)
    end
end

Base.@propagate_inbounds function apply_accumulation!(::AntiSymmetric, f, acc, src_acc_r, dst_acc_r, edge)
    _acc_dst = view(acc, dst_acc_r)
    _acc_src = view(acc, src_acc_r)

    for i in 1:length(dst_acc_r)
        _acc_dst[i] = f(_acc_dst[i],  edge[i])
        _acc_src[i] = f(_acc_src[i], -edge[i])
    end
end

Base.@propagate_inbounds function apply_accumulation!(::Symmetric, f, _acc, src_acc_r, dst_acc_r, edge)
    _acc_dst = view(_acc, dst_acc_r)
    _acc_src = view(_acc, src_acc_r)

    for i in 1:length(dst_acc_r)
        _acc_dst[i] = f(_acc_dst[i], edge[i])
        _acc_src[i] = f(_acc_src[i], edge[i])
    end
end

Base.@propagate_inbounds function apply_accumulation!(::Directed, f, _acc, _, dst_acc_r, edge)
    _acc_dst = view(_acc, dst_acc_r)

    for i in 1:length(dst_acc_r)
        _acc_dst[i] = f(_acc_dst[i], edge[i])
    end
end

Base.@propagate_inbounds function apply_accumulation!(::Fiducial, f, _acc, src_acc_r, dst_acc_r, edge)
    _acc_dst = view(_acc, dst_acc_r)
    _acc_src = view(_acc, src_acc_r)
    offset = Int(length(edge) / 2)

    for i in 1:length(dst_acc_r)
        _acc_dst[i] = f(_acc_dst[i], edge[i])
        _acc_src[i] = f(_acc_src[i], edge[offset + i])
    end
end


text:
function nlp_cvx_102_010(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - linear and nonlinear constraints
    # Variants
    #   010 - constraint intersection
    #   011 - one binding constraint (linear)
    #   012 - one binding constraint (nonlienar)
    #   013 - one binding constraint (quadratic objective)
    #   014 - no binding constraints (quadratic objective)

    model = Model(optimizer)

    @variable(model, x)
    @variable(model, y)

    @objective(model, Min, -x)
    @NLconstraint(model, x^2 + y^2 <= 1.0)
    @constraint(model, x + y >= 1.2)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, -0.974165743715913, tol = objective_tol)
    return check_solution(
        [x, y],
        [0.974165743715913, 0.2258342542139504],
        tol = primal_tol,
    )
end


text:
using Plots
using Printf
using SpectralMethodFFT

function x_k_init(c::ConfigFFT)

    ngrid = c.ngrids[1]
    ix_mid = ngrid÷2 + 1

    x_u = XFunc(zeros(ngrid), c)
    x_u[ix_mid] = 1.
    k_u = k_x(x_u)
    return x_u, k_u

end

function k_develop_k(k_u::KFunc, dt::Real)

    c = k_u.config
    xlen, = xlens_xranges(c.xranges)
    k_K = k_Kgen(c)

    return exp(- (2π*k_K/xlen)^2 * dt) * k_u

end

function main()

    # -------------
    # Parameters
    # -------------
    # *** space ***
    xrange = (-10., 10.)
    ngrid = 128
    c = ConfigFFT((ngrid,), (xrange,))
    # *** time ***
    t_st = 0.
    t_ed = 1.
    nt = 200
    dt = (t_ed - t_st) / nt

    # -------------
    # Time Integration
    # -------------
    x_u, k_u = x_k_init(c)
    x_X = x_Xgen(c)

    anim = Animation()
    for it = 0:nt
        t = @sprintf("%4.3f", t_st + it*dt)
        x_u = x_k(k_u)
        plt = plot(
                x_X, x_u,
                title="t=$t",
                label="heat",
                ylims=(-0.1, 1.1),
                xlabel="x", ylabel="u")
        frame(anim, plt)
        if it < nt
            k_u = k_develop_k(k_u, dt)
        end
    end
    gif(anim, "heat_1d.gif", fps=20)

end

main()

max_stars_repo_name: mat-der-D/SpectralMethodFFTW

max_issues_repo_name: mat-der-D/SpectralMethodFFTW

max_forks_repo_name: mat-der-D/SpectralMethodFFTW


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(w::GeoWeights)
  @series begin
    label --> "weights"
    domain(w), collect(w)
  end
end


text:
using CmdStan, Turing, ContinuousBenchmarks

include(joinpath(ContinuousBenchmarks.STAN_DATA_DIR, "bernoulli-stan.data.jl"))
include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "bernoulli-stan.model.jl"))

include(joinpath(ContinuousBenchmarks.DATA_DIR, "toy-data", "gdemo-stan.data.jl"))
include(joinpath(ContinuousBenchmarks.MODELS_DIR, "toy-models", "gdemo-stan.model.jl"))

stan_model_name = "simplegauss"
simplegaussstan = Stanmodel(name=stan_model_name, model=simplegaussstanmodel, nchains=1);

rc, simple_gauss_stan_sim = stan(simplegaussstan, simplegaussstandata, CmdStanDir=ContinuousBenchmarks.CMDSTAN_HOME, summary=false);

s_stan = simple_gauss_stan_sim[1:1000, ["s"], :].value[:]
m_stan = simple_gauss_stan_sim[1:1000, ["m"], :].value[:]
sg_time = get_stan_time(stan_model_name)