text:
#######################################################################
# Onsite/Hopping selectors
#######################################################################
abstract type Selector{M,S} end

struct OnsiteSelector{M,S} <: Selector{M,S}
    region::M
    sublats::S  # NTuple{N,NameType} (unresolved) or Vector{Int} (resolved on a lattice)
end

struct HoppingSelector{M,S,D,T} <: Selector{M,S}
    region::M
    sublats::S  # NTuple{N,Pair{NameType,NameType}} (unres) or Vector{Pair{Int,Int}} (res)
    dns::D
    range::T
end

onsiteselector(; region = missing, sublats = missing) =
    OnsiteSelector(region, sanitize_sublats(sublats))

hoppingselector(; region = missing, sublats = missing, dn = missing, range = missing) =
    HoppingSelector(region, sanitize_sublatpairs(sublats), sanitize_dn(dn), sanitize_range(range))

sanitize_sublats(s::Missing) = missing
sanitize_sublats(s::Integer) = (nametype(s),)
sanitize_sublats(s::NameType) = (s,)
sanitize_sublats(s::Tuple) = nametype.(s)
sanitize_sublats(s::Tuple{}) = ()
sanitize_sublats(n) = throw(ErrorException(
    "`sublats` for `onsite` must be either `missing`, a sublattice name or a tuple of names, see `onsite` for details"))

sanitize_sublatpairs(s::Missing) = missing
sanitize_sublatpairs(p::Pair{<:Union{NameType,Integer}, <:Union{NameType,Integer}}) =
    (ensurenametype(p),)
sanitize_sublatpairs(pp::Pair) =
    ensurenametype.(tupletopair.(tupleproduct(first(pp), last(pp))))
sanitize_sublatpairs(ps::Tuple) = tuplejoin(sanitize_sublatpairs.(ps)...)
sanitize_sublatpairs(s) = throw(ErrorException(
    "`sublats` for `hopping` must be either `missing` or a series of sublattice `Pair`s, see `hopping` for details"))

ensurenametype((s1, s2)::Pair) = nametype(s1) => nametype(s2)

sanitize_dn(dn::Missing) = missing
sanitize_dn(dn::Tuple{Vararg{Number,N}}) where {N} = (_sanitize_dn(dn),)
sanitize_dn(dn) = (_sanitize_dn(dn),)
sanitize_dn(dn::Tuple) = _sanitize_dn.(dn)
_sanitize_dn(dn::Tuple{Vararg{Number,N}}) where {N} = SVector{N,Int}(dn)
_sanitize_dn(dn::SVector{N}) where {N} = SVector{N,Int}(dn)
_sanitize_dn(dn::Vector) = SVector{length(dn),Int}(dn)

sanitize_range(::Missing) = missing
sanitize_range(range::Real) = isfinite(range) ? float(range) + sqrt(eps(float(range))) : float(range)

sublats(s::OnsiteSelector{<:Any,Missing}, lat::AbstractLattice) = collect(1:nsublats(lat))

function sublats(s::OnsiteSelector{<:Any,<:Tuple}, lat::AbstractLattice)
    names = lat.unitcell.names
    ss = Int[]
    for name in s.sublats
        i = findfirst(isequal(name), names)
        i !== nothing && push!(ss, i)
    end
    return ss
end

function sublats(s::HoppingSelector{<:Any,<:Tuple}, lat::AbstractLattice)
    names = lat.unitcell.names
    ss = Pair{Int,Int}[]
    for (nsrc, ndst) in s.sublats
        isrc = findfirst(isequal(nsrc), names)
        idst = findfirst(isequal(ndst), names)
        isrc !== nothing && idst !== nothing && push!(ss, isrc => idst)
    end
    return ss
end

sublats(s::HoppingSelector{<:Any,Missing}, lat::AbstractLattice) =
    tupletopair.(vec(collect(Iterators.product(1:nsublats(lat), 1:nsublats(lat)))))

# selector already resolved for a lattice
sublats(s::Selector{<:Any,<:Vector}, lat) = s.sublats

# API

resolve(s::HoppingSelector, lat::AbstractLattice) =
    HoppingSelector(s.region, sublats(s, lat), _checkdims(s.dns, lat), s.range)
resolve(s::OnsiteSelector, lat::AbstractLattice) = OnsiteSelector(s.region, sublats(s, lat))

_checkdims(dns::Missing, lat::AbstractLattice{E,L}) where {E,L} = dns
_checkdims(dns::Tuple{Vararg{SVector{L,Int}}}, lat::AbstractLattice{E,L}) where {E,L} = dns
_checkdims(dns, lat::AbstractLattice{E,L}) where {E,L} =
    throw(DimensionMismatch("Specified cell distance `dn` does not match lattice dimension $L"))

# are sites at (i,j) and (dni, dnj) or (dn, 0) selected?
(s::OnsiteSelector)(lat::AbstractLattice, (i, j)::Tuple, (dni, dnj)::Tuple{SVector, SVector}) =
    isonsite((i, j), (dni, dnj)) && isinregion(i, dni, s.region, lat) && isinsublats(sublat(lat, i), s.sublats)

(s::HoppingSelector)(lat::AbstractLattice, inds, dns) =
    !isonsite(inds, dns) && isinregion(inds, dns, s.region, lat) && isindns(dns, s.dns) &&
    isinrange(inds, s.range, lat) && isinsublats(tupletopair(sublat.(Ref(lat), inds)), s.sublats)

isonsite((i, j), (dni, dnj)) = i == j && dni == dnj

isinregion(i::Int, dn, ::Missing, lat) = true
isinregion(i::Int, dn, region::Function, lat) = region(sites(lat)[i] + bravais(lat) * dn)

isinregion(is::Tuple{Int,Int}, dns, ::Missing, lat) = true
function isinregion((row, col)::Tuple{Int,Int}, (dnrow, dncol), region::Function, lat)
    br = bravais(lat)
    r, dr = _rdr(sites(lat)[col] + br * dncol, sites(lat)[row] + br * dnrow)
    return region(r, dr)
end

isinsublats(s::Int, ::Missing) = true
isinsublats(s::Int, sublats::Vector{Int}) = s in sublats
isinsublats(ss::Pair{Int,Int}, ::Missing) = true
isinsublats(ss::Pair{Int,Int}, sublats::Vector{Pair{Int,Int}}) = ss in sublats
isinsublats(s, sublats) =
    throw(ArgumentError("Sublattices $sublats in selector are not resolved."))

isindns((dnrow, dncol)::Tuple{SVector,SVector}, dns) = isindns(dnrow - dncol, dns)
isindns(dn::SVector{L,Int}, dns::Tuple{Vararg{SVector{L,Int}}}) where {L} = dn in dns
isindns(dn::SVector, dns::Missing) = true
isindns(dn, dns) =
    throw(ArgumentError("Cell distance dn in selector is incompatible with Lattice."))

isinrange(inds, ::Missing, lat) = true
isinrange((row, col)::Tuple{Int,Int}, range::Number, lat) =
    norm(sites(lat)[col] - sites(lat)[row]) <= range

# merge non-missing fields of s´ into s
merge_non_missing(s::OnsiteSelector, s´::OnsiteSelector) =
    OnsiteSelector(merge_non_missing.(
        (s.region,  s.sublats),
        (s´.region, s´.sublats))...)
merge_non_missing(s::HoppingSelector, s´::HoppingSelector) =
    HoppingSelector(merge_non_missing.(
        (s.region,  s.sublats,  s.dns,  s.range),
        (s´.region, s´.sublats, s´.dns, s´.range))...)
merge_non_missing(o, o´::Missing) = o
merge_non_missing(o, o´) = o´

Base.adjoint(s::OnsiteSelector) = s
function Base.adjoint(s::HoppingSelector)
    # is_unconstrained_selector(s) &&
    #     @warn("Taking the adjoint of an unconstrained hopping is likely unintended")
    region´ = _adjoint(s.region)
    sublats´ = _adjoint.(s.sublats)
    dns´ = _adjoint.(s.dns)
    range = s.range
    return HoppingSelector(region´, sublats´, dns´, range)
end

_adjoint(::Missing) = missing
_adjoint(f::Function) = (r, dr) -> f(r, -dr)
_adjoint(t::Pair) = reverse(t)
_adjoint(t::SVector) = -t

# is_unconstrained_selector(s::HoppingSelector{Missing,Missing,Missing}) = true
# is_unconstrained_selector(s::HoppingSelector) = false

#######################################################################
# Tightbinding types
#######################################################################
abstract type TightbindingModelTerm end
abstract type AbstractOnsiteTerm <: TightbindingModelTerm end
abstract type AbstractHoppingTerm <: TightbindingModelTerm end

struct TightbindingModel{N,T<:Tuple{Vararg{TightbindingModelTerm,N}}}
    terms::T
end

struct OnsiteTerm{F,S<:OnsiteSelector,C} <: AbstractOnsiteTerm
    o::F
    selector::S
    coefficient::C
end

struct HoppingTerm{F,S<:HoppingSelector,C} <: AbstractHoppingTerm
    t::F
    selector::S
    coefficient::C
end

#######################################################################
# TightbindingModel API
#######################################################################
terms(t::TightbindingModel) = t.terms

TightbindingModel(ts::TightbindingModelTerm...) = TightbindingModel(ts)

(m::TightbindingModel)(r, dr) = sum(t -> t(r, dr), m.terms)

# External API #

Base.:*(x, m::TightbindingModel) = TightbindingModel(x .* m.terms)
Base.:*(m::TightbindingModel, x) = x * m
Base.:-(m::TightbindingModel) = TightbindingModel((-1) .* m.terms)

Base.:+(m::TightbindingModel, t::TightbindingModel) = TightbindingModel((m.terms..., t.terms...))
Base.:-(m::TightbindingModel, t::TightbindingModel) = m + (-t)

function Base.show(io::IO, m::TightbindingModel{N}) where {N}
    ioindent = IOContext(io, :indent => "  ")
    print(io, "TightbindingModel{$N}: model with $N terms", "\n")
    foreach(t -> print(ioindent, t, "\n"), m.terms)
end

#######################################################################
# TightbindingModelTerm API
#######################################################################
OnsiteTerm(t::OnsiteTerm, os::OnsiteSelector) =
    OnsiteTerm(t.o, os, t.coefficient)

(o::OnsiteTerm{<:Function})(r,dr) = o.coefficient * o.o(r)
(o::OnsiteTerm)(r,dr) = o.coefficient * o.o

HoppingTerm(t::HoppingTerm, os::HoppingSelector) =
    HoppingTerm(t.t, os, t.coefficient)

(h::HoppingTerm{<:Function})(r, dr) = h.coefficient * h.t(r, dr)
(h::HoppingTerm)(r, dr) = h.coefficient * h.t

sublats(t::TightbindingModelTerm, lat) = sublats(t.selector, lat)

sublats(m::TightbindingModel) = (t -> t.selector.sublats).(terms(m))

displayparameter(::Type{<:Function}) = "Function"
displayparameter(::Type{T}) where {T} = "$T"

function Base.show(io::IO, o::OnsiteTerm{F}) where {F}
    i = get(io, :indent, "")
    print(io,
"$(i)OnsiteTerm{$(displayparameter(F))}:
$(i)  Sublattices      : $(o.selector.sublats === missing ? "any" : o.selector.sublats)
$(i)  Coefficient      : $(o.coefficient)")
end

function Base.show(io::IO, h::HoppingTerm{F}) where {F}
    i = get(io, :indent, "")
    print(io,
"$(i)HoppingTerm{$(displayparameter(F))}:
$(i)  Sublattice pairs : $(h.selector.sublats === missing ? "any" : h.selector.sublats)
$(i)  dn cell distance : $(h.selector.dns === missing ? "any" : h.selector.dns)
$(i)  Hopping range    : $(round(h.selector.range, digits = 6))
$(i)  Coefficient      : $(h.coefficient)")
end

# External API #
"""
    onsite(o; region = missing, sublats = missing)

Create an `TightbindingModel` with a single `OnsiteTerm` that applies an onsite energy `o`
to a `Lattice` when creating a `Hamiltonian` with `hamiltonian`.

The onsite energy `o` can be a number, a matrix (preferably `SMatrix`), a `UniformScaling`
(e.g. `3*I`) or a function of the form `r -> ...` for a position-dependent onsite energy.

The dimension of `o::AbstractMatrix` must match the orbital dimension of applicable
sublattices (see also `orbitals` option for `hamiltonian`). If `o::UniformScaling` it will
be converted to an identity matrix of the appropriate size when applied to multiorbital
sublattices. Similarly, if `o::SMatrix` it will be truncated or padded to the appropriate
size.

    onsite(model::TightbindingModel; kw...)

Return a `TightbindingModel` with only the onsite terms of `model`. Any non-missing `kw` is
applied to all such terms.

# Keyword arguments

Only sites at position `r` in sublattice with name `s::NameType` will be selected if
`region(r) && s in sublats` is true. Any missing `region` or `sublat` will not be used to
constraint the selection.

The keyword `sublats` allows the following formats:

    sublats = :A           # Onsite on sublat :A only
    sublats = (:A,)        # Same as above
    sublats = (:A, :B)     # Onsite on sublat :A and :B

# Combining models

`OnsiteTerm`s and `HoppingTerm`s created with `onsite` or `hopping` can added or substracted
together or be multiplied by scalars to build more complicated `TightbindingModel`s, e.g.
`onsite(1) - 3 * hopping(2)`

# Examples

```jldoctest
julia> model = onsite(1, sublats = (:A,:B)) - 2 * hopping(2, sublats = :A=>:A)
TightbindingModel{2}: model with 2 terms
  OnsiteTerm{Int64}:
    Sublattices      : (:A, :B)
    Coefficient      : 1
  HoppingTerm{Int64}:
    Sublattice pairs : (:A => :A,)
    dn cell distance : any
    Hopping range    : 1.0
    Coefficient      : -2

julia> newmodel = onsite(model; sublats = :A) + hopping(model)
TightbindingModel{2}: model with 2 terms
  OnsiteTerm{Int64}:
    Sublattices      : (:A,)
    Coefficient      : 1
  HoppingTerm{Int64}:
    Sublattice pairs : (:A => :A,)
    dn cell distance : any
    Hopping range    : 1.0
    Coefficient      : -2

julia> LatticePresets.honeycomb() |> hamiltonian(onsite(r->@SMatrix[1 2; 3 4]), orbitals = Val(2))
Hamiltonian{<:Lattice} : Hamiltonian on a 2D Lattice in 2D space
  Bloch harmonics  : 1 (SparseMatrixCSC, sparse)
  Harmonic size    : 2 × 2
  Orbitals         : ((:a, :a), (:a, :a))
  Element type     : 2 × 2 blocks (Complex{Float64})
  Onsites          : 2
  Hoppings         : 0
  Coordination     : 0.0
```

# See also:
    `hopping`
"""
onsite(o; kw...) = onsite(o, onsiteselector(; kw...))

onsite(o, selector::Selector) = TightbindingModel(OnsiteTerm(o, selector, 1))

onsite(m::TightbindingModel, selector::Selector) =
    TightbindingModel(_onlyonsites(selector, m.terms...))

_onlyonsites(s, t::OnsiteTerm, args...) =
    (OnsiteTerm(t, merge_non_missing(t.selector, s)), _onlyonsites(s, args...)...)
_onlyonsites(s, t::HoppingTerm, args...) = (_onlyonsites(s, args...)...,)
_onlyonsites(s) = ()

"""
    hopping(t; region = missing, sublats = missing, dn = missing, range = 1, plusadjoint = false)

Create an `TightbindingModel` with a single `HoppingTerm` that applies a hopping `t` to a
`Lattice` when creating a `Hamiltonian` with `hamiltonian`.

The hopping amplitude `t` can be a number, a matrix (preferably `SMatrix`), a
`UniformScaling` (e.g. `3*I`) or a function of the form `(r, dr) -> ...` for a
position-dependent hopping (`r` is the bond center, and `dr` the bond vector). If `sublats`
is specified as a sublattice name pair, or tuple thereof, `hopping` is only applied between
sublattices with said names.

The dimension of `t::AbstractMatrix` must match the orbital dimension of applicable
sublattices (see also `orbitals` option for `hamiltonian`). If `t::UniformScaling` it will
be converted to a (possibly rectangular) identity matrix of the appropriate size when
applied to multiorbital sublattices. Similarly, if `t::SMatrix` it will be truncated or
padded to the appropriate size.

`OnsiteTerm`s and `HoppingTerm`s created with `onsite` or `hopping` can be added or
substracted together to build more complicated `TightbindingModel`s.

    hopping(model::TightbindingModel; kw...)

Return a `TightbindingModel` with only the hopping terms of `model`. Any non-missing `kw` is
applied to all such terms.

# Keyword arguments

Only hoppings between two sites at positions `r₁ = r - dr/2` and `r₂ = r + dr`, belonging to
unit cells at integer distance `dn´` and to sublattices `s₁` and `s₂` will be selected if:
`region(r, dr) && s in sublats && dn´ in dn && norm(dr) <= range`. If any of these is
`missing` it will not be used to constraint the selection. Note that the default `range` is
1, not `missing`.

The keyword `dn` can be a `Tuple`/`Vector`/`SVector` of `Int`s, or a tuple thereof. The
keyword `sublats` allows the following formats:

    sublats = :A => :B                 # Hopping from :A to :B sublattices
    sublats = (:A => :B,)              # Same as above
    sublats = (:A => :B, :C => :D)     # Hopping from :A to :B or :C to :D
    sublats = (:A, :C) .=> (:B, :D)    # Broadcasted pairs, same as above
    sublats = (:A, :C) => (:B, :D)     # Direct product, (:A=>:B, :A=:D, :C=>:B, :C=>D)

The keyword `plusadjoint` produces a model with the input hopping term plus its adjoint.
Note that this substitution is made before multiplying by any coefficient, so that
`im*hopping(..., plusadjoint = true) == im*(hopping(...) + hopping(...)')`.

# Combining models

`OnsiteTerm`s and `HoppingTerm`s created with `onsite` or `hopping` can added or substracted
together or be multiplied by scalars to build more complicated `TightbindingModel`s, e.g.
`onsite(1) - 3 * hopping(2)`

# Examples

```jldoctest
julia> model = 3 * onsite(1) - hopping(2, dn = ((1,2), (0,0)), sublats = :A=>:B)
TightbindingModel{2}: model with 2 terms
  OnsiteTerm{Int64}:
    Sublattices      : any
    Coefficient      : 3
  HoppingTerm{Int64}:
    Sublattice pairs : (:A => :B,)
    dn cell distance : ([1, 2], [0, 0])
    Hopping range    : 1.0
    Coefficient      : -1

julia> newmodel = onsite(model) + hopping(model, range = 2)
TightbindingModel{2}: model with 2 terms
  OnsiteTerm{Int64}:
    Sublattices      : any
    Coefficient      : 3
  HoppingTerm{Int64}:
    Sublattice pairs : (:A => :B,)
    dn cell distance : ([1, 2], [0, 0])
    Hopping range    : 2.0
    Coefficient      : -1

julia> LatticePresets.honeycomb() |> hamiltonian(hopping((r,dr) -> cos(r[1]), sublats = (:A,:B) => (:A,:B)))
Hamiltonian{<:Lattice} : Hamiltonian on a 2D Lattice in 2D space
  Bloch harmonics  : 7 (SparseMatrixCSC, sparse)
  Harmonic size    : 2 × 2
  Orbitals         : ((:a,), (:a,))
  Element type     : scalar (Complex{Float64})
  Onsites          : 0
  Hoppings         : 18
  Coordination     : 9.0
```

# See also:
    `onsite`
"""
function hopping(t; plusadjoint = false, range = 1, kw...)
    hop = hopping(t, hoppingselector(; range = range, kw...))
    return plusadjoint ? hop + hop' : hop
end
hopping(t, selector) = TightbindingModel(HoppingTerm(t, selector, 1))

hopping(m::TightbindingModel, selector::Selector) =
    TightbindingModel(_onlyhoppings(selector, m.terms...))

_onlyhoppings(s, t::OnsiteTerm, args...) = (_onlyhoppings(s, args...)...,)
_onlyhoppings(s, t::HoppingTerm, args...) =
    (HoppingTerm(t, merge_non_missing(t.selector, s)), _onlyhoppings(s, args...)...)
_onlyhoppings(s) = ()

Base.:*(x, o::OnsiteTerm) =
    OnsiteTerm(o.o, o.selector, x * o.coefficient)
Base.:*(x, t::HoppingTerm) = HoppingTerm(t.t, t.selector, x * t.coefficient)
Base.:*(t::TightbindingModelTerm, x) = x * t
Base.:-(t::TightbindingModelTerm) = (-1) * t

Base.adjoint(t::TightbindingModel) = TightbindingModel(adjoint.(terms(t)))
Base.adjoint(t::OnsiteTerm{Function}) = OnsiteTerm(r -> t.o(r)', t.selector, t.coefficient')
Base.adjoint(t::OnsiteTerm) = OnsiteTerm(t.o', t.selector, t.coefficient')
Base.adjoint(t::HoppingTerm{Function}) = HoppingTerm((r, dr) -> t.t(r, -dr)', t.selector', t.coefficient')
Base.adjoint(t::HoppingTerm) = HoppingTerm(t.t', t.selector', t.coefficient')

#######################################################################
# offdiagonal
#######################################################################
"""
    offdiagonal(model, lat, nsublats::NTuple{N,Int})

Build a restricted version of `model` that applies only to off-diagonal blocks formed by
sublattice groups of size `nsublats`.
"""
offdiagonal(m::TightbindingModel, lat, nsublats) =
    TightbindingModel(offdiagonal.(m.terms, Ref(lat), Ref(nsublats)))

offdiagonal(o::OnsiteTerm, lat, nsublats) =
    throw(ArgumentError("No onsite terms allowed in off-diagonal coupling"))

function offdiagonal(t::HoppingTerm, lat, nsublats)
    selector´ = resolve(t.selector, lat)
    s = selector´.sublats
    sr = sublatranges(nsublats...)
    filter!(spair ->  findblock(first(spair), sr) != findblock(last(spair), sr), s)
    return HoppingTerm(t.t, selector´, t.coefficient)
end

sublatranges(i::Int, is::Int...) = _sublatranges((1:i,), is...)
_sublatranges(rs::Tuple, i::Int, is...) = _sublatranges((rs..., last(last(rs)) + 1: last(last(rs)) + i), is...)
_sublatranges(rs::Tuple) = rs

findblock(s, sr) = findfirst(r -> s in r, sr)

#######################################################################
# @onsite! and @hopping!
#######################################################################
abstract type ElementModifier{N,F,S} end

struct ParametricFunction{N,F,P<:Val}
    f::F
    params::P
end

ParametricFunction{N}(f::F, p::P) where {N,F,P} = ParametricFunction{N,F,P}(f, p)

(pf::ParametricFunction)(args...; kw...) = pf.f(args...; kw...)

struct OnsiteModifier{N,F<:ParametricFunction{N},S<:Selector} <: ElementModifier{N,F,S}
    f::F
    selector::S
end

OnsiteModifier(f, selector) = OnsiteModifier(f, selector, false)

struct HoppingModifier{N,F<:ParametricFunction{N},S<:Selector} <: ElementModifier{N,F,S}
    f::F
    selector::S
end

HoppingModifier(f, selector) = HoppingModifier(f, selector, false)

const UniformModifier = ElementModifier{1}
const UniformHoppingModifier = HoppingModifier{1}
const UniformOnsiteModifier = OnsiteModifier{1}

"""
    parameters(p::ElementModifier...)

Return the parameter names for one or several  `ElementModifier` created with `@onsite!` or
`@hopping!`
"""
parameters(ms::ElementModifier...) = mergetuples(_parameters.(ms)...)
_parameters(m::ElementModifier) = _parameters(m.f)
_parameters(pf::ParametricFunction) = pf.params

"""
    @onsite!(args -> body; kw...)

Create an `ElementModifier`, to be used with `parametric`, that applies `f = args -> body`
to onsite energies specified by `kw` (see `onsite` for details  on possible `kw`s). The form
of `args -> body` may be `(o; params...) -> ...` or `(o, r; params...) -> ...` if the
modification is position (`r`) dependent. Keyword arguments `params` are optional, and
include any parameters that `body` depends on that the user may want to tune.

# See also:
    `@hopping!`, `parametric`
"""
macro onsite!(kw, f)
    f, N, params = get_f_N_params(f, "Only @onsite!(args -> body; kw...) syntax supported")
    return esc(:(Quantica.OnsiteModifier(Quantica.ParametricFunction{$N}($f, $(Val(params))), Quantica.onsiteselector($kw))))
end

macro onsite!(f)
    f, N, params = get_f_N_params(f, "Only @onsite!(args -> body; kw...) syntax supported")
    return esc(:(Quantica.OnsiteModifier(Quantica.ParametricFunction{$N}($f, $(Val(params))), Quantica.onsiteselector())))
end

"""
    @hopping!(args -> body; kw...)

Create an `ElementModifier`, to be used with `parametric`, that applies `f = args -> body`
to hoppings energies specified by `kw` (see `hopping` for details on possible `kw`s). The
form of `args -> body` may be `(t; params...) -> ...` or `(t, r, dr; params...) -> ...` if
the modification is position (`r`, `dr`) dependent. Keyword arguments `params` are optional,
and include any parameters that `body` depends on that the user may want to tune.

# See also:
    `@onsite!`, `parametric`
"""
macro hopping!(kw, f)
    f, N, params = get_f_N_params(f, "Only @hopping!(args -> body; kw...) syntax supported")
    return esc(:(Quantica.HoppingModifier(Quantica.ParametricFunction{$N}($f, $(Val(params))), Quantica.hoppingselector($kw))))
end

macro hopping!(f)
    f, N, params = get_f_N_params(f, "Only @hopping!(args -> body; kw...) syntax supported")
    return esc(:(Quantica.HoppingModifier(Quantica.ParametricFunction{$N}($f, $(Val(params))), Quantica.hoppingselector())))
end

# Extracts normalized f, number of arguments and kwarg names from an anonymous function f
function get_f_N_params(f, msg)
    (f isa Expr && f.head == :->) || throw(ArgumentError(msg))
    d = ExprTools.splitdef(f)
    kwargs = convert(Vector{Any}, get!(d, :kwargs, []))
    d[:kwargs] = kwargs  # in case it wasn't Vector{Any} originally
    if isempty(kwargs)
        params = ()
        push!(kwargs, :(_...))  # normalization : append _... to kwargs
    else
        params = get_kwname.(kwargs) |> Tuple
        if !isempty(params) && last(params) == :...
            params = tuplemost(params)  # drop _... kwarg from params
        else
            push!(kwargs, :(_...))  # normalization : append _... to kwargs
        end
    end
    N = haskey(d, :args) ? length(d[:args]) : 0
    f´ = ExprTools.combinedef(d)
    return f´, N, params
end

get_kwname(x::Symbol) = x
get_kwname(x::Expr) = x.head === :kw ? x.args[1] : x.head  # x.head == :...

resolve(o::OnsiteModifier, lat) = OnsiteModifier(o.f, resolve(o.selector, lat))

resolve(h::HoppingModifier, lat) = HoppingModifier(h.f, resolve(h.selector, lat))

resolve(t::Tuple, lat) = _resolve(lat, t...)
_resolve(lat, t, ts...) = (resolve(t, lat), _resolve(lat, ts...)...)
_resolve(lat) = ()

# Intended for resolved ElementModifier{N} only. The N is the number of arguments accepted.
@inline (o!::UniformOnsiteModifier)(o, r; kw...) = o!(o; kw...)
@inline (o!::UniformOnsiteModifier)(o, r, dr; kw...) = o!(o; kw...)
@inline (o!::UniformOnsiteModifier)(o; kw...) = o!.f(o; kw...)
@inline (o!::OnsiteModifier)(o, r, dr; kw...) = o!(o, r; kw...)
@inline (o!::OnsiteModifier)(o, r; kw...) = o!.f(o, r; kw...)

@inline (h!::UniformHoppingModifier)(t, r, dr; kw...) = h!(t; kw...)
@inline (h!::UniformHoppingModifier)(t; kw...) = h!.f(t; kw...)
@inline (h!::HoppingModifier)(t, r, dr; kw...) = h!.f(t, r, dr; kw...)


text:
function edge_shape(shape::Symbol=:line; center_offset::Real = 3, end_offsets::Tuple{Real, Real} = (0, 0))
    draw = (video, object, frame; shape=shape, center_offset=center_offset, end_offsets=end_offsets, p1, p2, from_vertex_bbx, to_vertex_bbx, kwargs...) -> begin
        # Calculate  edge segment outside of vertex boundaries
        # ToDo:
        #   1. For curved edges choose _p1 and _p2 using intersection circle-arc
        #   2. For node objects other than circle find a intersection of polygon and line
        #   3. [FIX] When ther are self loops sometimes self loop passes through the node
        # Known bug: If the gap between nodes are too less neither of the intersection
        #   points lie within the line connecting p1 and p2 hence ispointline fails
        # Temporary fix: Make sure there is enough gaps within 2 node objects
        self_loop = get(Dict(collect(kwargs)...), :self_loop, false)
        _, ip1, ip2 = intersectionlinecircle(p1, p2, p1, distance(O, (from_vertex_bbx[2]-from_vertex_bbx[1])/2))
        _p1 = ispointonline(ip1, ip2, (p1+p2)/2; extended=false) ? ip1 : ip2
        _, ip1, ip2 = intersectionlinecircle(p1, p2, p2, distance(O, (to_vertex_bbx[2]-to_vertex_bbx[1])/2))
        _p2 = ispointonline(ip1, ip2, (p1+p2)/2; extended=false) ? ip1 : ip2
        if self_loop
            neighbor_vertices = neighbors(GRAPHS[object.opts[:_graph_idx]], GRAPH_EDGES[object.opts[:_edge_idx]].from_vertex.vertex_id, strict=true)
            neighbor_pos = Vector{Point}(map(nv -> nv.object.opts[:position], neighbor_vertices))
            idx, sector = largest_circle_sector(p1, neighbor_pos)
            angle = slope(p1, neighbor_pos[idx])+sector/2
            direction = Point(cos(angle), sin(angle))
            outline = draw_self_loop_edge(p1, from_vertex_bbx, center_offset, end_offsets, direction)
        elseif shape == :curved
            outline = draw_curved_edge(_p1, _p2, center_offset, end_offsets)
        else
            outline = draw_straight_edge(_p1, _p2, center_offset, end_offsets)
        end
        poly(outline, :stroke)
        @register_style_opts object shape center_offset end_offsets outline
    end
    return :edge_shape=>draw
end

function edge_style(; color="red", linewidth::Real=2, dash::String="dotted", blend="on")
    styling = true
    draw = (video, object, frame; linewidth=linewidth, dash=dash, color=color, blend=blend, styling=styling, kwargs...) -> begin
        # ToDo: Allow a list of colors, linewidths, dash to be specified.
        # Interpolate in between when blend option is turned on
        sethue(color)
        setline(linewidth)
        setdash(dash)
        @register_style_opts object color linewidth dash blend styling
    end
    return :edge_style=>draw
end

function edge_arrow(; start=nothing, finish=nothing, color="black")
    draw = (video, object, frame; arrow_start=start, arrow_finish=finish, arrow_color=color, outline, kwargs...) ->  begin
        sethue(arrow_color)
        # Skip the default case
        if outline == [O, O]
            return
        end
        if is_directed(GRAPHS[object.opts[:_graph_idx]].graph.adjacency_graph) && arrow_finish === nothing
            arrow_finish = true
        end
        l = length(outline)
        if arrow_start !== nothing
            Luxor.arrow(reverse(outline[1:2])...)
        end
        if arrow_finish !== nothing
            Luxor.arrow(outline[l-1:l]...)
        end
        @register_style_opts object arrow_start arrow_finish arrow_color
    end
    return :edge_arrow=>draw
end

function edge_text_style(; color="black", fonttype="Times Roman", fontsize=10, opacity=0.8)
    draw = (video, object, frame; text_color=color, fonttype=fonttype, fontsize=fontsize, text_opacity=opacity, kwargs...) -> begin
        sethue(text_color)
        setfont(fonttype, fontsize)
        setopacity(text_opacity)
        @register_style_opts object text_color fonttype fontsize text_opacity
    end
    return :edge_text_style=>draw
end

function edge_text(text_content; position::Real = 0.5, offset::Real = 2)
    draw = (video, object, frame; text_content=text_content, position=position, offset=offset, outline, kwargs...) ->  begin
        # ToDo: What to do for bezier curved edges?
        # ToDo: Add edge label rotation
        # Align text along edge
        if outline == [O, O]
            return
        end
        l = length(outline)
        idx = floor(Int, l * position)
        if l%2 == 0
            pt = (outline[idx] + outline[idx+1])/2
            sl = slope(outline[idx], outline[idx+1])
        else
            pt = outline[idx]
            sl = idx == l ? slope(outline[idx%l+1], outline[idx]) : slope(outline[idx], outline[idx+1])
        end
        translate(pt)
        rotate(sl)
        label(text_content, :N, O, leader=false, offset=offset)
        sl == 0 ? rotate(1/sl) : nothing
        translate(-pt)
        @register_style_opts object text_content position offset
    end
    return :edge_text=>draw
end

function draw_self_loop_edge(p1, from_vertex_bbx, center_offset, end_offsets, direction)
    r = max((from_vertex_bbx[2]-from_vertex_bbx[1])/2..., center_offset)
    offset = min((from_vertex_bbx[2]-from_vertex_bbx[1])/2..., center_offset)
    c = p1 + r * direction
    # Might replace by a clipping region to remove the extra edge remnants appearing inside a vertex
    # pt1, pt2 = intersectioncircleboundingbox(c, r, from_vertex_bbx...)
    # A temporary solution unless intersectioncircleboundingbox is fixed
    _, pt1, pt2 = intersectioncirclecircle(c, r, p1, distance(O, (from_vertex_bbx[2]-from_vertex_bbx[1])/2))
    pt1, pt2 = isarcclockwise(c, pt1, pt2) ? (pt1, pt2) : (pt2, pt1)
    angles = slope(c, pt2), slope(c, pt1)
    angle_offsets = end_offsets./r
    new_angles = angles .+ (angle_offsets[1], -angle_offsets[2])
    Luxor.newpath()
    Luxor.arc(c, r, new_angles..., :path)
    Luxor.closepath()
    return first(pathtopoly())
end

function draw_curved_edge(p1, p2, center_offset, end_offsets)
    center_pt = Luxor.perpendicular((p1+p2)/2, p1, center_offset)
    c, r = center3pts(p1, center_pt, p2)
    angles = slope(c, p1), slope(c, p2)
    angle_offsets = end_offsets./r
    new_angles = angles .+ (angle_offsets[1], -angle_offsets[2])
    Luxor.newpath()
    Luxor.arc(c, r, new_angles..., :path)
    Luxor.closepath()
    return first(pathtopoly())
end

function draw_straight_edge(p1, p2, center_offset, end_offsets)
    d = Luxor.distance(p1, p2)
    t1 = end_offsets[1]/d
    t2 = (d-end_offsets[2])/d
    new_p1 = t1*p1 + (1-t1)*p2
    new_p2 = t2*p1 + (1-t2)*p2
    Luxor.newpath()
    line(new_p1, new_p2, :path)
    Luxor.closepath()
    return first(pathtopoly())
end

function largest_circle_sector(c::Point, pts::Vector{Point})
    idx = 1
    sector = 0
    n = length(pts)
    if n==1
        return 1, 2*pi
    end
    slopes = sort(Vector{Tuple}([(slope(c, pts[i]), i) for i in 1:n]))
    for i in 1:n
        s = slopes[i%n+1][1]-slopes[i][1]+(i==n)*2*pi
        if sector < s
            sector = s
            idx = slopes[i][2]
        end
    end
    return idx, sector
end

# # ToDo: Needs fixes
# Add this to Luxor
# function intersectioncircleboundingbox(c::Point, radius::Real, upperleft::Point, lowerright::Point)
#     # Find the closest corner of box to the center of the circle
#     px = abs(lowerright[1] - c[1]) < abs(upperleft[1] - c[1]) ? lowerright[1] : upperleft[1]
#     py = abs(lowerright[2] - c[2]) < abs(upperleft[2] - c[2]) ? lowerright[2] : upperleft[2]
#     # Find the other 2 corner points joining this point
#     o1 = Point(upperleft[1]+lowerright[1]-px, py)
#     o2 = Point(px, upperleft[2]+lowerright[2]-py)
#     # Calculate intersection points
#     nints, ip1, ip2 = zip(intersectionlinecircle.(Point(px, py), [o1, o2], c, radius)...)
#     p = Point(px, py)
#     intersection_points = ispointonline.(ip1, [p, p], [o1, o2]; extended=false).*ip1 + ispointonline.(ip2, [p, p], [o1, o2]; extended=false).*ip2
#     # throw(ErrorException("$nints, $ip1, $ip2, $intersection_points"))
#     return intersection_points
# end

ext: jl


text:
module TestParameterTypes

using Mimi
using Test

import Mimi:
    external_params, external_param, TimestepMatrix, TimestepVector,
    ArrayModelParameter, ScalarModelParameter, FixedTimestep, build, import_params!

#
# Test that parameter type mismatches are caught
#
expr = :(
    @defcomp BadComp1 begin
        a = Parameter(index=[time, regions], default=[10, 11, 12])  # should be 2D default
        function run_timestep(p, v, d, t)
        end
    end
)
@test_throws LoadError eval(expr)

expr = :(
    @defcomp BadComp2 begin
        a = Parameter(default=[10, 11, 12])  # should be scalar default
        function run_timestep(p, v, d, t)
        end
    end
)
@test_throws LoadError eval(expr)

#
# Test that the old type parameterization syntax errors
#
expr = :(
    @defcomp BadComp3 begin
        a::Int = Parameter()
        function run_timestep(p, v, d, t)
        end
    end
)
eval(expr)  # Just a deprecation warning for v0.10, then will change to error in v1.0
# @test_throws LoadError eval(expr)

@defcomp MyComp begin
    a = Parameter(index=[time, regions], default=ones(101,3))
    b = Parameter(index=[time], default=1:101)
    c = Parameter(index=[regions])
    d = Parameter()
    e = Parameter(index=[four])
    f = Parameter{Array{Float64, 2}}()
    g = Parameter{Int}(default=10.0)    # value should be Int despite Float64 default
    h = Parameter(default=10)           # should be "numtype", despite Int default
    j = Parameter{Int}(index = [regions])

    function run_timestep(p, v, d, t)
    end
end

# Check that explicit number type for model works as expected
numtype = Float32
arrtype = Union{Missing, numtype}

m = Model(numtype)

set_dimension!(m, :time, 2000:2100)
set_dimension!(m, :regions, 3)
set_dimension!(m, :four, 4)

add_comp!(m, MyComp)
set_param!(m, :MyComp, :c, [4,5,6])
set_param!(m, :MyComp, :d, 0.5)   # 32-bit float constant
set_param!(m, :MyComp, :e, [1,2,3,4])
set_param!(m, :MyComp, :f, reshape(1:16, 4, 4))
set_param!(m, :MyComp, :j, [1,2,3])

build(m)    # applies defaults, creating external params
extpars = external_params(m)

# TBD: These are not (yet) external params. Defaults are applied at build time.
@test isa(extpars[:a], ArrayModelParameter)
@test isa(extpars[:b], ArrayModelParameter)

@test isa(extpars[:c], ArrayModelParameter)
@test isa(extpars[:d], ScalarModelParameter)
@test isa(extpars[:e], ArrayModelParameter)
@test isa(extpars[:f], ScalarModelParameter) # note that :f is stored as a scalar parameter even though its values are an array

@test typeof(extpars[:a].values) == TimestepMatrix{FixedTimestep{2000, 1}, arrtype, 1}
@test typeof(extpars[:b].values) == TimestepVector{FixedTimestep{2000, 1}, arrtype}

@test typeof(extpars[:c].values) == Array{arrtype, 1}
@test typeof(extpars[:d].value) == numtype
@test typeof(extpars[:e].values) == Array{arrtype, 1}
@test typeof(extpars[:f].value) == Array{Float64, 2}
@test typeof(extpars[:g].value) <: Int
@test typeof(extpars[:h].value) == numtype

# test updating parameters
@test_throws ErrorException update_param!(m, :a, 5) # expects an array
@test_throws ErrorException update_param!(m, :a, ones(101)) # wrong size
@test_throws ErrorException update_param!(m, :a, fill("hi", 101, 3)) # wrong type

update_param!(m, :a, Array{Int,2}(zeros(101, 3))) # should be able to convert from Int to Float

@test_throws ErrorException update_param!(m, :d, ones(5)) # wrong type; should be scalar
update_param!(m, :d, 5) # should work, will convert to float
@test extpars[:d].value == 5
@test_throws ErrorException update_param!(m, :e, 5) # wrong type; should be array
@test_throws ErrorException update_param!(m, :e, ones(10)) # wrong size
update_param!(m, :e, [4,5,6,7])

@test length(extpars) == 9
@test typeof(extpars[:a].values) == TimestepMatrix{FixedTimestep{2000, 1}, arrtype, 1}

@test typeof(extpars[:d].value) == numtype
@test typeof(extpars[:e].values) == Array{arrtype, 1}

#------------------------------------------------------------------------------
# Test updating TimestepArrays with update_param!
#------------------------------------------------------------------------------

@defcomp MyComp2 begin
    x=Parameter(index=[time])
    y=Variable(index=[time])
    function run_timestep(p,v,d,t)
        v.y[t]=p.x[t]
    end
end

# 1. Test with Fixed Timesteps

m = Model()
set_dimension!(m, :time, 2000:2002)
add_comp!(m, MyComp2) # ; first=2000, last=2002)
set_param!(m, :MyComp2, :x, [1, 2, 3])

# N.B. `first` and `last` are now disabled.
# Can't move last beyond last for a component
# @test_throws ErrorException set_dimension!(m, :time, 2001:2003)

set_dimension!(m, :time, 2001:2002)

update_param!(m, :x, [4, 5, 6], update_timesteps = false)
x = external_param(m.md, :x)
@test x.values isa Mimi.TimestepArray{Mimi.FixedTimestep{2000, 1, LAST} where LAST, Union{Missing,Float64}, 1}
@test x.values.data == [4., 5., 6.]
# TBD: this fails, but I'm not sure how it's supposed to behave. It says:
# (ERROR: BoundsError: attempt to access 3-element Array{Float64,1} at index [4])
# run(m)
# @test m[:MyComp2, :y][1] == 5   # 2001
# @test m[:MyComp2, :y][2] == 6   # 2002

update_param!(m, :x, [2, 3], update_timesteps = true)
x = external_param(m.md, :x)
@test x.values isa Mimi.TimestepArray{Mimi.FixedTimestep{2001, 1, LAST} where LAST, Union{Missing,Float64}, 1}
@test x.values.data == [2., 3.]
run(m)
@test m[:MyComp2, :y][1] == 2   # 2001
@test m[:MyComp2, :y][2] == 3   # 2002

# 2. Test with Variable Timesteps

m = Model()
set_dimension!(m, :time, [2000, 2005, 2020])
add_comp!(m, MyComp2)
set_param!(m, :MyComp2, :x, [1, 2, 3])

set_dimension!(m, :time, [2005, 2020, 2050])

update_param!(m, :x, [4, 5, 6], update_timesteps = false)
x = external_param(m.md, :x)
@test x.values isa Mimi.TimestepArray{Mimi.VariableTimestep{(2000, 2005, 2020)}, Union{Missing,Float64}, 1}
@test x.values.data == [4., 5., 6.]
#run(m)
#@test m[:MyComp2, :y][1] == 5   # 2005
#@test m[:MyComp2, :y][2] == 6   # 2020

update_param!(m, :x, [2, 3, 4], update_timesteps = true)
x = external_param(m.md, :x)
@test x.values isa Mimi.TimestepArray{Mimi.VariableTimestep{(2005, 2020, 2050)}, Union{Missing,Float64}, 1}
@test x.values.data == [2., 3., 4.]
run(m)
@test m[:MyComp2, :y][1] == 2   # 2005
@test m[:MyComp2, :y][2] == 3   # 2020

# 3. Test updating from a dictionary

m = Model()
set_dimension!(m, :time, [2000, 2005, 2020])
add_comp!(m, MyComp2)
set_param!(m, :MyComp2, :x, [1, 2, 3])

set_dimension!(m, :time, [2005, 2020, 2050])

update_params!(m, Dict(:x=>[2, 3, 4]), update_timesteps = true)
x = external_param(m.md, :x)
@test x.values isa Mimi.TimestepArray{Mimi.VariableTimestep{(2005, 2020, 2050)}, Union{Missing,Float64}, 1}
@test x.values.data == [2., 3., 4.]
run(m)
@test m[:MyComp2, :y][1] == 2   # 2005
@test m[:MyComp2, :y][2] == 3   # 2020
@test m[:MyComp2, :y][3] == 4   # 2050

# 4. Test updating the time index to a different length

m = Model()
set_dimension!(m, :time, 2000:2002)     # length 3
add_comp!(m, MyComp2)
set_param!(m, :MyComp2, :x, [1, 2, 3])

set_dimension!(m, :time, 1999:2003)     # length 5

@test_throws ErrorException update_param!(m, :x, [2, 3, 4, 5, 6], update_timesteps = false)
update_param!(m, :x, [2, 3, 4, 5, 6], update_timesteps = true)
x = external_param(m.md, :x)
@test x.values isa Mimi.TimestepArray{Mimi.FixedTimestep{1999, 1, LAST} where LAST, Union{Missing,Float64}, 1}
@test x.values.data == [2., 3., 4., 5., 6.]

run(m)
@test m[:MyComp2, :y] == [2., 3., 4., 5., 6.]

# 5. Test all the warning and error cases

@defcomp MyComp3 begin
    regions=Index()
    x=Parameter(index=[time])       # One timestep array parameter
    y=Parameter(index=[regions])    # One non-timestep array parameter
    z=Parameter()                   # One scalar parameter
end

m = Model()                             # Build the model
set_dimension!(m, :time, 2000:2002)     # Set the time dimension
set_dimension!(m, :regions, [:A, :B])
add_comp!(m, MyComp3)
set_param!(m, :MyComp3, :x, [1, 2, 3])
set_param!(m, :MyComp3, :y, [10, 20])
set_param!(m, :MyComp3, :z, 0)

@test_throws ErrorException update_param!(m, :x, [1, 2, 3, 4]) # Will throw an error because size
@test_throws ErrorException update_param!(m, :y, [10, 15], update_timesteps=true) # Not a timestep array
update_param!(m, :y, [10, 15])
@test external_param(m.md, :y).values == [10., 15.]
@test_throws ErrorException update_param!(m, :z, 1, update_timesteps=true) # Scalar parameter
update_param!(m, :z, 1)
@test external_param(m.md, :z).value == 1

# Reset the time dimensions
set_dimension!(m, :time, 2005:2007)

update_params!(m, Dict(:x=>[3,4,5], :y=>[10,20], :z=>0), update_timesteps=true) # Won't error when updating from a dictionary

@test external_param(m.md, :x).values isa Mimi.TimestepArray{Mimi.FixedTimestep{2005,1},Union{Missing,Float64},1}
@test external_param(m.md, :x).values.data == [3.,4.,5.]
@test external_param(m.md, :y).values == [10.,20.]
@test external_param(m.md, :z).value == 0

#------------------------------------------------------------------------------
# Test the three different set_param! methods for a Symbol type parameter
#------------------------------------------------------------------------------

@defcomp A begin
    p1 = Parameter{Symbol}()
end

function _get_model()
    m = Model()
    set_dimension!(m, :time, 10)
    add_comp!(m, A)
    return m
end

# Test the 3-argument version of set_param!
m = _get_model()
@test_throws MethodError set_param!(m, :p1, 3)  # Can't set it with an Int

set_param!(m, :p1, :foo)    # Set it with a Symbol
run(m)
@test m[:A, :p1] == :foo

# Test the 4-argument version of set_param!
m = _get_model()
@test_throws MethodError set_param!(m, :A, :p1, 3)

set_param!(m, :A, :p1, :foo)
run(m)
@test m[:A, :p1] == :foo

# Test the 5-argument version of set_param!
m = _get_model()
@test_throws MethodError set_param!(m, :A, :p1, :A_p1, 3)

set_param!(m, :A, :p1, :A_p1, :foo)
run(m)
@test m[:A, :p1] == :foo

#------------------------------------------------------------------------------
# Test that if set_param! errors in the connection step,
#       the created param doesn't remain in the model's list of params
#------------------------------------------------------------------------------

@defcomp A begin
    p1 = Parameter(index = [time])
end

@defcomp B begin
    p1 = Parameter(index = [time])
end

m = Model()
set_dimension!(m, :time, 10)
add_comp!(m, A)
add_comp!(m, B)

@test_throws ErrorException set_param!(m, :p1, 1:5)     # this will error because the provided data is the wrong size
@test isempty(m.md.external_params)                     # But it should not be added to the model's dictionary

end #module


text:
## Old site

cd("/home/danielc/Documentos/GitHub/Julia-Para-Economistas/Julia Para Economistas")

using JuDoc

serve()

## New site

cd("/home/danielc/Documentos/GitHub/Julia-Para-Economistas/Julia P Economistas F")

using Franklin

serve()

max_stars_repo_name: danmrc/Julia-Para-Economistas

max_issues_repo_name: danmrc/Julia-Para-Economistas

max_forks_repo_name: danmrc/Julia-Para-Economistas


text:
## Differentiation
# Only k = 1 for now...
Derivative(S::OscLaurent{DD,RR},k::Integer) where {DD,RR} = ConcreteDerivative(S,k)
bandwidths(D::ConcreteDerivative{OscLaurent{DD,RR}})  where {DD,RR} = (3,3)
rangespace(D::ConcreteDerivative{S}) where {S<:OscLaurent}=D.space
function mob_derivative_getindex(d::PeriodicLine,m,j::Integer,k::Integer)
    # j is the row index, i.e. (j,k) not (k,j) as other ApproxFun implementations
    if j == 1 && 2 <= k <= 3
        -(-1)^k*0.5im/d.L
    elseif iseven(j)
        l = j÷2
        if k == j
            -1.0im*l/d.L
        elseif k == j+2
            -1.0im*(l+1)/2*1/d.L
        elseif k == j-2 && j > 2
            -1.0im*(l-1)/2*1/d.L
        else
            0.0im
        end
    else
        l = j÷2
        if k == j
            1.0im*l/d.L
        elseif k == j+2
            1.0im*(l+1)/2*1/d.L
        elseif k == j-2 && j > 3
            1.0im*(l-1)/2*1/d.L
        else
            0.0im
        end
    end
end

getindex(D::ConcreteDerivative{OscLaurent{DD,RR},OT,T},k::Integer,j::Integer) where {DD,RR,OT,T} =
convert(T,(k == j ? domainspace(D).exp*1im : 0.) + mob_derivative_getindex(domain(D),D.order,k,j))

## Multiplication, same as Laurent
Multiplication(f::Fun{OscLaurent{DD,RR}},sp::OscLaurent{DD,RR}) where {DD,RR} = ConcreteMultiplication(cfstype(f),f,sp)

function laurent_getindex(negative::AbstractVector{T},nonnegative::AbstractVector{T},k::Integer,j::Integer) where T
    # switch to double-infinite indices
    k=iseven(k) ? -k÷2 : (k-1)÷2
    j=iseven(j) ? -j÷2 : (j-1)÷2

    if 0<k-j≤length(negative)
        negative[k-j]
    elseif 0≤j-k≤length(nonnegative)-1
        nonnegative[j-k+1]
    else
        zero(T)
    end
end

rangespace(T::ConcreteMultiplication{OscLaurent{DD,RR},OscLaurent{DD,RR}}) where {DD,RR} = OscLaurent(domainspace(T).domain,T.space.exp + space(T.f).exp)
function getindex(T::ConcreteMultiplication{OscLaurent{DD,RR},OscLaurent{DD,RR}},k::Integer,j::Integer) where {DD,RR}
    isempty(T.f.coefficients) && return zero(eltype(T))
    laurent_getindex(T.f.coefficients[3:2:end],T.f.coefficients[[1;2:2:end]],k,j)

end

function bandwidths(T::ConcreteMultiplication{OscLaurent{DD,RR},OscLaurent{DD,RR}}) where {DD,RR}
    bbi = blockbandwidths(T)
    (2bbi[1]+1,2bbi[2]+1)
end

function blockbandwidths(T::ConcreteMultiplication{OscLaurent{DD,RR},OscLaurent{DD,RR}}) where {DD,RR}
    m = ncoefficients(T.f)÷2
    (m,m)
end

# old code for testing
function fouriertransform(f::Fun{T}) where T <: OscLaurent
    α = f.space.exp
    L = f.space.domain.L
    sp1 = Laguerre(1.0,Ray(α,0.0,1/(2.0*L),true))
    sp2 = Laguerre(1.0,Ray(α,π,1/(2.0*L),true))
    c1 = f.coefficients[3:2:end]
    c2 = f.coefficients[2:2:end]
    c1 = c1.*[4*π*L*(-1)^(i+1) for i=1:length(c1)]
    c2 = c2.*[4*π*L*(-1)^(i+1) for i=1:length(c2)]
    Fun(LaguerreWeight(0.0,0.5,sp1),c1) + Fun(LaguerreWeight(0.0,0.5,sp2),c2)
end

#Op class = FourierOperator
#Op = FourierTransform
#ConcOp = ConcreteFourierTransform
#WrappOp = FourierTransformWrapper

### Fourier transform

abstract type FourierOperator{S,OT,T} <: Operator{T} end

macro fourier_operator(Op)
    ConcOp = Meta.parse("Concrete"*string(Op))
    OpWrap = Meta.parse(string(Op)*"Wrapper")
    return esc(quote
        abstract type $Op{SSS,OT,TTT} <: FourierOperator{SSS,OT,TTT} end

        struct $ConcOp{S<:Space,OT,T} <: $Op{S,OT,T}
            space::S
            sign::OT
        end

        struct $OpWrap{BT<:Operator,S<:Space,OT,T} <: $Op{S,OT,T}
            op::BT
            sign::OT
        end

        @wrapper $OpWrap

        $ConcOp(sp::Space,k) = $ConcOp{typeof(sp),typeof(k),prectype(sp)}(sp,k)
        $ConcOp(sp::Space) = $ConcOp(sp,1.0)

        $Op(sp::UnsetSpace,k) = $ConcOp(sp,k)
        $Op(sp::UnsetSpace,k::Real) = $ConcOp(sp,k)

        # not needed yet
        function Base.convert(::Type{Operator{T}},D::$ConcOp) where T
            if T==eltype(D)
            D
            else
                $ConcOp{typeof(D.space),typeof(D.sign),T}(D.space,D.sign)
            end
        end

        function Base.convert(::Type{Operator{T}},D::$OpWrap) where T
            if T==eltype(D)
            D
            else
                # work around typeinfernece bug
                op=convert(Operator{T},D.op)
                $OpWrap{typeof(op),typeof(domainspace(op)),T}(op)
            end
        end

        ApproxFunBase.domain(D::$ConcOp) = domain(D.space)
        ApproxFunBase.domainspace(D::$ConcOp) = D.space
        Base.getindex(::$ConcOp{UnsetSpace,T},k::Integer,j::Integer) where {OT,T} =
            error("Spaces cannot be inferred for operator")
        ApproxFunBase.rangespace(D::$ConcOp{UnsetSpace,OT,T}) where {OT,T} = UnsetSpace()

        ApproxFunBase.promotedomainspace(D::$Op,sp::UnsetSpace) = D
        function ApproxFunBase.promotedomainspace(D::$Op,sp::Space)
            if isambiguous(domain(sp))
                $Op(typeof(sp)(domain(D)),D.sign)
            else
                $Op(sp,D.sign)
            end
        end
    end)
end
hasconversion(a::LaguerreWeight,b::DiracSpace) = false
hasconversion(b::DiracSpace,a::LaguerreWeight) = false
#maxspace(a::SumSpace,b::DiracSpace) = SumSpace(a,b)
#maxspace(b::DiracSpace,a::SumSpace) = SumSpace(b,a)

@fourier_operator(δFourierTransform)
@fourier_operator(SFourierTransform)

abstract type AbstractFourierTransform end
struct FourierTransform <:AbstractFourierTransform
    A
    B
end
FourierTransform(k) = FourierTransform(SFourierTransform(k),δFourierTransform(k))

*(F::AbstractFourierTransform,f::Fun{T}) where {T <: OscLaurent} = F.A*f + F.B*f
*(F::AbstractFourierTransform,f::Fun{T}) where {T <: LaguerreWeight} = F.A*f + F.B*f
*(F::AbstractFourierTransform,f::Fun{T}) where {T <: DiracSpace} = F.A*f + F.B*f
*(F::AbstractFourierTransform,f::Fun{T}) where {T <: PiecewiseSpace} = sum(map(x ->F*x,components(f)))
*(F::AbstractFourierTransform,f::Fun{T}) where {T <: SumSpace} = sum(map(x ->F*x,components(f)))
*(F::AbstractFourierTransform,f::SumFun) = +(map(x -> SumFun(F*x),f.funs)...)
*(F::AbstractFourierTransform,f::Array{T}) where T = map(x -> F*x,f)

# BEGIN: FourierTransform of OscLaurent space

SFourierTransform(sp::OscLaurent,k) = ConcreteSFourierTransform(sp,k)
SFourierTransform(k::Float64) = ConcreteSFourierTransform(UnsetSpace(),k)
SFourierTransform() = ConcreteSFourierTransform(UnsetSpace(),1.0)

δFourierTransform(sp::OscLaurent,k) = ConcreteδFourierTransform(sp,k)
δFourierTransform(k::Float64) = ConcreteδFourierTransform(UnsetSpace(),k)
δFourierTransform() = ConcreteδFourierTransform(UnsetSpace(),1.0)

choosedomainspace(M::FourierOperator{UnsetSpace},sp::Space) =
    iswrapper(M) ? choosedomainspace(M.op,sp) : sp

function rangespace(D::ConcreteSFourierTransform{S,OT,T}) where {S<:OscLaurent,OT,T}
    α = D.space.exp/D.sign
    L = D.space.domain.L*abs(D.sign)
    sp1 = Laguerre(1.0,Ray(α,0.0,1/(2.0*L),true))
    sp2 = Laguerre(1.0,Ray(α,π,1/(2.0*L),true))
    sp0 = DiracSpace(D.space.exp/D.sign)
    if D.sign > 0.
        PiecewiseSpace(sp0,LaguerreWeight(0.0,0.5,sp2),LaguerreWeight(0.0,0.5,sp1))
    else
        PiecewiseSpace(sp0,LaguerreWeight(0.0,0.5,sp1),LaguerreWeight(0.0,0.5,sp2))
    end
end

function rangespace(D::ConcreteδFourierTransform{S,OT,T}) where {S<:OscLaurent,OT,T}
    α = D.space.exp/D.sign
    L = D.space.domain.L*abs(D.sign)
    sp1 = Laguerre(1.0,Ray(α,0.0,1/(2.0*L),true))
    sp2 = Laguerre(1.0,Ray(α,π,1/(2.0*L),true))
    sp0 = DiracSpace(D.space.exp/D.sign)
    if D.sign > 0.
        PiecewiseSpace(sp0,LaguerreWeight(0.0,0.5,sp2),LaguerreWeight(0.0,0.5,sp1))
    else
        PiecewiseSpace(sp0,LaguerreWeight(0.0,0.5,sp1),LaguerreWeight(0.0,0.5,sp2))
    end
end

function getindex(D::ConcreteδFourierTransform{S,OT,T},k::Integer,j::Integer) where {S<:OscLaurent,OT,T}
    if k == 1
        j = iseven(j) ? -j÷2 : (j-1)÷2
        return (-1)^j
    else
        return zero(T)
    end
 end

# Base.size(A::ConcreteδFourierTransform{S,OT,T},k::Integer)  where {S<:OscLaurent,OT,T} = k==1 ? 1 : dimension(domainspace(A))
 bandwidths(D::ConcreteδFourierTransform{OscLaurent{DD,RR},OT,T}) where {DD,RR,OT,T} = 0,∞
 function osclaurent_ft_getindex(k::Integer,j::Integer,L::Float64,T::Type)
     if j == 1 || j != k
         return zero(T)
     else
         # switch to double-infinite indices
         # k=iseven(k) ? -k÷2 : (k-1)÷2
         j = iseven(j) ? -j÷2 : (j-1)÷2
         return -4*π*L*(-1)^(j)
     end
 end
 getindex(D::ConcreteSFourierTransform{OscLaurent{DD,RR},OT,T},k::Integer,j::Integer) where {DD,RR,OT,T} = osclaurent_ft_getindex(k,j,D.space.domain.L,T)
 bandwidths(D::ConcreteSFourierTransform{OscLaurent{DD,RR},OT,T}) where {DD,RR,OT,T} = 0,0

# END: FourierTransform of OscLaurent space

# BEGIN:  FourierTransform of DiracSpace

δFourierTransform(sp::DiracSpace,k) = ConcreteδFourierTransform(sp,k)
SFourierTransform(sp::DiracSpace,k) = ConcreteSFourierTransform(sp,k)

function rangespace(D::ConcreteδFourierTransform{S,OT,T}) where {S<:DiracSpace,OT,T}
    @assert length(D.space.points) == 1
    α = -D.space.points[1]*D.sign
    L = 1.0
    OscConstantSpace(α,L)
end

function rangespace(D::ConcreteSFourierTransform{S,OT,T}) where {S<:DiracSpace,OT,T}
    @assert length(D.space.points) == 1
    α = -D.space.points[1]*D.sign
    L = 1.0
    OscConstantSpace(α,L)
end

function getindex(D::ConcreteδFourierTransform{S,OT,T},k::Integer,j::Integer) where {S<:DiracSpace,OT,T}
     if k == j == 1
         return one(T)
     else
         return zero(T)
     end
end

function getindex(D::ConcreteSFourierTransform{S,OT,T},k::Integer,j::Integer) where {S<:DiracSpace,OT,T}
     return zero(T)
end

bandwidths(D::ConcreteSFourierTransform{S,OT,T}) where {S<:DiracSpace,OT,T} = 0,0
bandwidths(D::ConcreteδFourierTransform{S,OT,T}) where {S<:DiracSpace,OT,T} = 0,0
# END: FourierTransform of DiracSpace

# BEGIN: FourierTransform of LaguerreWeight
SFourierTransform(sp::LaguerreWeight,k) = ConcreteSFourierTransform(sp,k)
δFourierTransform(sp::LaguerreWeight,k) = ConcreteδFourierTransform(sp,k)

function rangespace(D::ConcreteδFourierTransform{S,OT,T}) where {S<:LaguerreWeight,OT,T}
    sp = D.space
    @assert sp.α == 0.0 && sp.L == 0.5
    α = -sp.space.domain.center*D.sign
    L = abs(1/(D.sign*2*sp.space.domain.L))
    OscLaurent(α,L)
end

function rangespace(D::ConcreteSFourierTransform{S,OT,T}) where {S<:LaguerreWeight,OT,T}
    sp = D.space
    @assert sp.α == 0.0 && sp.L == 0.5
    α = -sp.space.domain.center*D.sign
    L = abs(1/(D.sign*2*sp.space.domain.L))
    OscLaurent(α,L)
end

function getindex(D::ConcreteSFourierTransform{S,OT,T},k::Integer,j::Integer) where {S<:LaguerreWeight,OT,T}
    sp = D.space
    L = 1/(D.sign*2*sp.space.domain.L)
    ang = angle(sp.space.domain)
    @assert ang ≈ 0.0 || ang ≈ π
    if (ang ≈ 0.0 && D.sign < 0) || (ang ≈ π && D.sign > 0)
        σ = 1
    else
        σ = 0
    end

    δ = D.sign < 0 ? 0 : 1

    if k > 1 && k != 2j + σ#(ang ≈ 0.0 ? 1 : 0)
        return zero(T)
    elseif k == 1
        return -1/(4*π*L*(-1)^(j+δ))
    else
        return 1/(4*π*L*(-1)^(j+δ))
    end
end

function getindex(D::ConcreteδFourierTransform{S,OT,T},k::Integer,j::Integer) where {S<:LaguerreWeight,OT,T}
     return zero(T)
end

bandwidths(D::ConcreteSFourierTransform{S,OT,T}) where {S<:LaguerreWeight,OT,T} = ∞,3
bandwidths(D::ConcreteδFourierTransform{S,OT,T}) where {S<:LaguerreWeight,OT,T} = 0,0
israggedbelow(D::ConcreteSFourierTransform{S,OT,T}) where {S<:LaguerreWeight,OT,T} = true
colstop(D::ConcreteSFourierTransform{S,OT,T},j::Integer) where {S<:LaguerreWeight,OT,T} = isodd(j) ? 2*j+1 : 2*j+1 # can be refined to include the zero columns


text:
# SPDX-License-Identifier: MIT

# Construct a block diagonal matrix with CUDA CSR matrix.
# Adapted from:
#   https://github.com/JuliaLang/julia/blob/1b93d53fc4bb59350ada898038ed4de2994cce33/stdlib/SparseArrays/src/sparsematrix.jl#L3396
function SparseArrays.blockdiag(X::CuSparseMatrixCSR{T}...) where T
    num = length(X)
    mX = Int[ size(x, 1) for x in X ]
    nX = Int[ size(x, 2) for x in X ]
    m = sum(mX) # number of rows
    n = sum(nX) # number of cols

    rowPtr = CuVector{Cint}(undef, m+1)
    nnzX = Int[ nnz(x) for x in X ]
    nnz_res = sum(nnzX)
    colVal = CuVector{Cint}(undef, nnz_res)
    nzVal = CuVector{T}(undef, nnz_res)

    nnz_sofar = 0
    nX_sofar = 0
    mX_sofar = 0
    for i = 1:num
        rowPtr[ (1:mX[i]+1) .+ mX_sofar ] = X[i].rowPtr .+ nnz_sofar
        colVal[ (1:nnzX[i]) .+ nnz_sofar ] = X[i].colVal .+ nX_sofar
        nzVal[ (1:nnzX[i]) .+ nnz_sofar ] = nonzeros(X[i])
        nnz_sofar += nnzX[i]
        nX_sofar += nX[i]
        mX_sofar += mX[i]

    end
    CUDA.@allowscalar rowPtr[m+1] = nnz_sofar + 1

    CuSparseMatrixCSR{T}(rowPtr, colVal, nzVal, (m, n))
end

function Base.vcat(X::CuSparseVector{T}...) where T
    num = length(X)
    nX = Int[ length(x) for x in X ]
    n = sum(nX)

    nnzX = Int[ nnz(x) for x in X ]
    nnz_total = sum(nnzX)
    iPtr = CuVector{Cint}(undef, nnz_total)
    nzVal = CuVector{T}(undef, nnz_total)

    nX_sofar = 0
    nnz_sofar = 0
    for i = 1:num
        iPtr[ (1:nnzX[i]) .+ nnz_sofar ] = X[i].iPtr .+ nX_sofar
        nzVal[ (1:nnzX[i]) .+ nnz_sofar ] = nonzeros(X[i])
        nX_sofar += nX[i]
        nnz_sofar += nnzX[i]
    end

    CuSparseVector{T}(iPtr, nzVal, n)
end

max_forks_repo_name: BUTSpeechFIT/ProbabilisticFSMs


text:
module ContextVariablesX

# Re-exporting `Base` functions so that Documenter knows what's public:
export @contextvar, ContextVar, get, getindex, snapshot_context, with_context

using Logging: AbstractLogger, Logging
using UUIDs: UUID, UUIDs, uuid4

if isdefined(UUIDs, :uuid5)
    using UUIDs: uuid5
else
    using Compat: uuid5
end

include("payloadlogger.jl")

function _ContextVar end

# `ContextVar` object itself does not hold any data (except the
# default value).  It is actually just a key into the task-local
# context storage.
"""
    ContextVar{T}

Context variable type.  This is the type of the object `var` created by
[`@contextvar var`](@ref @contextvar).  This acts as a reference to the
value stored in a task-local context.  The macro `@contextvar` is the only
public API to construct this object.

!!! warning

    It is unspecified if this type is concrete or not. It may be
    changed to an abstract type and/or include more type parameters in
    the future.
"""
struct ContextVar{T}
    name::Symbol
    _module::Module
    key::UUID
    has_default::Bool
    default::T

    global _ContextVar
    _ContextVar(name, _module, key, ::Type{T}, default) where {T} =
        new{T}(name, _module, key, true, default)
    _ContextVar(name, _module, key, ::Type{T}) where {T} = new{T}(name, _module, key, false)
end

_ContextVar(name, _module, key, ::Nothing, default) =
    _ContextVar(name, _module, key, typeof(default), default)
_ContextVar(name, _module, key, ::Nothing) = _ContextVar(name, _module, key, Any)

#=
baremodule _EmptyModule end

ContextVar{T}(name::Symbol, default) where {T} =
    _ContextVar(name, _EmptyModule, uuid4(), T, default)
ContextVar{T}(name::Symbol) where {T} = _ContextVar(name, _EmptyModule, uuid4(), T)

ContextVar(name::Symbol, default) = ContextVar{typeof(default)}(name, default)
ContextVar(name::Symbol) = ContextVar{Any}(name)

"""
    ContextVar{T}(name::Symbol [, default])
    ContextVar(name::Symbol [, default])

Define *local* context variable.

!!! warning

    Using this constructor to define context variable in the global name space
    (i.e., `const var = ContextVar(:var)`) is not recommended because this is not
    safe to be used with Distribute.jl.  Use `@contextvar var` instead.
"""
ContextVar
=#

Base.eltype(::Type{ContextVar{T}}) where {T} = T

function Base.show(io::IO, var::ContextVar)
    print(io, ContextVar)
    if eltype(var) !== Any && !(var.has_default && typeof(var.default) === eltype(var))
        print(io, '{', eltype(var), '}')
    end
    print(io, '(', repr(var.name))
    if var.has_default
        print(io, ", ")
        show(io, var.default)
    end
    print(io, ')')
end

function Base.show(io::IO, ::MIME"text/plain", var::ContextVar)
    print(io, var._module, '.', var.name, " :: ContextVar")
    if get(io, :compact, false) === false
        print(io, " [", var.key, ']')
        if get(var) === nothing
            print(io, " (not assigned)")
        else
            print(io, " => ")
            show(IOContext(io, :compact => true), MIME"text/plain"(), var[])
        end
    end
end

# The primitives that can be monkey-patched to play with new context
# variable storage types:
"""
    merge_ctxvars(ctx::Union{Nothing,T}, kvs) -> ctx′:: Union{Nothing,T}

!!! warning

    This is not a public API.  This documentation is for making it easier to
    experiment with different implementations of the context variable storage
    backend, by monkey-patching it at run-time.  When this function is
    monkey-patched, `ctxvars_type` should also be monkey-patched to return
    the type `T`.

The first argument `ctx` is either `nothing` or a dict-like object of type `T` where
its `keytype` is `UUID` and `valtype` is `Any`.  The second argument `kvs` is an
iterable of `Pair{UUID,<:Union{Some,Nothing}}` values.  Iterable `kvs` must have
length.

If `ctx` is `nothing` and `kvs` is non-empty, `merge_ctxvars` creates a new
instance of `T`. If `ctx` is not `nothing`, it returns a shallow-copy `ctx′` of
`ctx` where `k => v` is inserted to `ctx′` for each `k => Some(v)` in `kvs`
and `k` is deleted from `ctx′` for each `k => nothing` in `kvs`.
"""
function merge_ctxvars(ctx, kvs)
    # Assumption: eltype(kvs) <: Pair{UUID,<:Union{Some,Nothing}}
    if isempty(kvs)
        return ctx
    else
        # Copy-or-create-on-write:
        vars = ctx === nothing ? ctxvars_type()() : copy(ctx)
        for (k, v) in kvs
            if v === nothing
                delete!(vars, k)
            else
                vars[k] = something(v)
            end
        end
        isempty(vars) && return nothing  # should we?
        return vars
    end
end

ctxvars_type() = Dict{UUID,Any}
_ctxvars_type() = Union{Nothing,ctxvars_type()}
get_task_ctxvars() = _get_task_ctxvars()::_ctxvars_type()

new_merged_ctxvars(kvs) =
    merge_ctxvars(
        get_task_ctxvars(),
        (
            k.key => v === nothing ? v : Some(convert(eltype(k), something(v)))
            for (k, v) in kvs
        ),
    )::_ctxvars_type()

struct _NoValue end

"""
    get(var::ContextVar{T}) -> Union{Some{T},Nothing}

Return `Some(value)` if `value` is assigned to `var`.  Return `nothing` if
unassigned.
"""
function Base.get(var::ContextVar{T}) where {T}
    ctx = get_task_ctxvars()
    if ctx === nothing
        var.has_default && return Some(var.default)
        return nothing
    end
    if var.has_default
        return Some(get(ctx, var.key, var.default)::T)
    else
        y = get(ctx, var.key, _NoValue())
        y isa _NoValue || return Some(ctx[var.key]::T)
    end
    return nothing
end

"""
    getindex(var::ContextVar{T}) -> value::T

Return the `value` assigned to `var`.  Throw a `KeyError` if unassigned.
"""
function Base.getindex(var::ContextVar{T}) where {T}
    maybe = get(var)
    maybe === nothing && throw(KeyError(var))
    return something(maybe)::T
end

"""
    genkey(__module__::Module, varname::Symbol) -> Union{UUID,Nothing}.

Generate a stable UUID for a context variable `__module__.\$varname`.
"""
function genkey(__module__::Module, varname::Symbol)
    fullpath = push!(collect(fullname(__module__)), varname)
    if any(x -> occursin(".", string(x)), fullpath)
        throw(ArgumentError(
            "Modules and variable names must not contain a dot:\n" * join(fullpath, "\n"),
        ))
    end
    pkgid = Base.PkgId(__module__)
    if pkgid.uuid === nothing
        return nothing
    end
    return uuid5(pkgid.uuid, join(fullpath, '.'))
end

"""
    @contextvar var[::T] [= default]

Declare a context variable named `var`.  The type constraint `::T` and
the default value `= default` are optional.  If the default value is given
without the type constraint `::T`, its type `T = typeof(default)` is used.

!!! warning

    Context variables defined outside a proper package does not work
    with `Distributed`.

# Examples

Top-level context variables needs to be declared in a package:

```julia
module MyPackage
@contextvar cvar1
@contextvar cvar2 = 1
@contextvar cvar3::Int
end
```
"""
macro contextvar(ex0)
    ex = ex0
    if Meta.isexpr(ex, :(=))
        length(ex.args) != 2 && throw(ArgumentError("Unsupported syntax:\n$ex0"))
        ex, default = ex.args
        args = Any[esc(default)]
    else
        args = []
    end
    if Meta.isexpr(ex, :(::))
        length(ex.args) != 2 && throw(ArgumentError("Malformed input:\n$ex0"))
        ex, vartype = ex.args
        pushfirst!(args, esc(vartype))
    else
        pushfirst!(args, nothing)
    end
    if !(ex isa Symbol)
        if ex === ex0
            throw(ArgumentError("Unsupported syntax:\n$ex0"))
        else
            throw(ArgumentError("""
                Not a variable name:
                $ex
                Input:
                $ex0
                """))
        end
    end
    varname = QuoteNode(ex)
    key = genkey(__module__, ex)
    if key === nothing
        if isinteractive() && _WARN_DYNAMIC_KEY[]
            @warn (
                "Using `@contextvar` outside a package is discouraged except" *
                " for interactive exploration or quick scripts.  Context" *
                " variables created outside a proper package do not have " *
                " consistent ID and cannot be used with `Distributed`."
            ) maxlog=1
        end
        # Creating a UUID at macro expansion time because:
        # * It would be a memory leak it were created at run-time because
        #   context variable storage can be filled with UUIDs created at
        #   run-time.
        # * Creating it at run-time is doable with function-based interface like
        #   `ContextVar(:name, default)`.
        key = uuid4()
    end
    return Expr(
        :const,
        :($(esc(ex)) = _ContextVar($varname, $__module__, $key, $(args...))),
    )
end

# A hack to make doctest easier
const _WARN_DYNAMIC_KEY = Ref(true)

"""
    with_context(f, var1 => value1, var2 => value2, ...)
    with_context(f, pairs)

Run `f` in a context with given values set to the context variables.  Variables
specified in this form are rolled back to the original value when `with_context`
returns.  It act like a dynamically scoped `let`.  If `nothing` is passed as
a value, corresponding context variable is cleared; i.e., it is unassigned or
takes the default value.  Use `Some(value)` to set `value` if `value` can be
`nothing`.

    with_context(f, nothing)

Run `f` in a new empty context.  All variables are rewind to the original values
when `with_context` returns.

Note that

```julia
var2[] = value2
with_context(var1 => value1) do
    @show var2[]  # shows value2
    var3[] = value3
end
@show var3[]  # shows value3
```

and

```julia
var2[] = value2
with_context(nothing) do
    var1[] = value1
    @show var2[]  # shows default (or throws)
    var3[] = value3
end
@show var3[]  # does not show value3
```

are not equivalent.
"""
with_context
with_context(f, kvs::Pair{<:ContextVar}...) = with_task_ctxvars(f, new_merged_ctxvars(kvs))
with_context(f, ::Nothing) = with_task_ctxvars(f, nothing)

struct ContextSnapshot{T}
    vars::T
end

# TODO: Do we need to implement `Dict{ContextVar}(::ContextSnapshot)`?
# This requires storing UUID-to-ContextVar mapping somewhere.

"""
    snapshot_context() -> snapshot::ContextSnapshot

Get a snapshot of a context that can be passed to [`with_context`](@ref) to
run a function inside the current context at later time.
"""
snapshot_context() = ContextSnapshot(get_task_ctxvars())

with_context(f, snapshot::ContextSnapshot) = with_task_ctxvars(f, snapshot.vars)

end # module


text:
"

 Himmelblau function

 author: Atsushi Sakai
"

function himmelblau_function(x, y)
    return (x^2 + y - 11)^2 + (x + y^2 - 7)^2
end

x = parse(Float64,ARGS[1])
y = parse(Float64,ARGS[2])
println(himmelblau_function(x, y))


text:
function layoutable(::Type{Textbox}, fig_or_scene; bbox = nothing, kwargs...)

    topscene = get_topscene(fig_or_scene)

    attrs = merge!(
        Attributes(kwargs),
        default_attributes(Textbox, topscene).attributes)

    @extract attrs (halign, valign, textsize, stored_string, placeholder,
        textcolor, textcolor_placeholder, displayed_string,
        boxcolor, boxcolor_focused_invalid, boxcolor_focused, boxcolor_hover,
        bordercolor, textpadding, bordercolor_focused, bordercolor_hover, focused,
        bordercolor_focused_invalid,
        borderwidth, cornerradius, cornersegments, boxcolor_focused,
        validator, restriction, cursorcolor, reset_on_defocus, defocus_on_submit)

    decorations = Dict{Symbol, Any}()

    layoutobservables = LayoutObservables{Textbox}(attrs.width, attrs.height,
        attrs.tellwidth, attrs.tellheight,
        halign, valign, attrs.alignmode; suggestedbbox = bbox)

    scenearea = lift(layoutobservables.computedbbox) do bb
        Rect(round.(Int, bb.origin), round.(Int, bb.widths))
    end

    scene = Scene(topscene, scenearea, raw = true, camera = campixel!)

    cursorindex = Node(0)
    ltextbox = Textbox(fig_or_scene, layoutobservables, attrs, decorations, cursorindex, nothing)

    bbox = lift(FRect2D ∘ AbstractPlotting.zero_origin, scenearea)

    roundedrectpoints = lift(roundedrectvertices, scenearea, cornerradius, cornersegments)

    displayed_string[] = isnothing(stored_string[]) ? placeholder[] : stored_string[]

    displayed_is_valid = lift(displayed_string, validator) do str, validator
        valid::Bool = validate_textbox(str, validator)
    end

    hovering = Node(false)

    realbordercolor = lift(bordercolor, bordercolor_focused,
        bordercolor_focused_invalid, bordercolor_hover, focused, displayed_is_valid, hovering,
        typ = Any) do bc, bcf, bcfi, bch, focused, valid, hovering

        if focused
            valid ? bcf : bcfi
        else
            hovering ? bch : bc
        end
    end

    realboxcolor = lift(boxcolor, boxcolor_focused,
        boxcolor_focused_invalid, boxcolor_hover, focused, displayed_is_valid, hovering,
        typ = Any) do bc, bcf, bcfi, bch, focused, valid, hovering

        if focused
            valid ? bcf : bcfi
        else
            hovering ? bch : bc
        end
    end

    box = poly!(topscene, roundedrectpoints, strokewidth = borderwidth,
        strokecolor = realbordercolor,
        color = realboxcolor, raw = true)
    decorations[:box] = box

    displayed_chars = @lift([c for c in $displayed_string])

    realtextcolor = lift(textcolor, textcolor_placeholder, focused, stored_string, displayed_string, typ = Any) do tc, tcph, foc, cont, disp
        # the textbox has normal text color if it's focused
        # if it's defocused, the displayed text has to match the stored text in order
        # to be normal colored
        if foc || cont == disp
            tc
        else
            tcph
        end
    end

    t = Label(scene, text = displayed_string, bbox = bbox, halign = :left, valign = :top,
        width = Auto(true), height = Auto(true), color = realtextcolor,
        textsize = textsize, padding = textpadding)

    displayed_charbbs = lift(t.layoutobservables.reportedsize) do sz
        charbbs(t.elements[:text])
    end

    cursorpoints = lift(cursorindex, displayed_charbbs) do ci, bbs

        if (ci) > length(bbs)
            # correct cursorindex if it's outside of the displayed charbbs range
            cursorindex[] = length(bbs)
            return
        end

        if ci == 0
            [leftline(bbs[1])...]
        else
            [rightline(bbs[ci])...]
        end
    end

    cursor = linesegments!(scene, cursorpoints, color = cursorcolor, linewidth = 2)

    cursoranimtask = nothing

    on(t.layoutobservables.reportedsize) do sz
        layoutobservables.autosize[] = sz
    end

    # trigger text for autosize
    t.text = displayed_string[]

    # trigger bbox
    layoutobservables.suggestedbbox[] = layoutobservables.suggestedbbox[]

    mousestate = addmouseevents!(scene)

    onmouseleftdown(mousestate) do state
        focus!(ltextbox)

        if displayed_string[] == placeholder[]
            displayed_string[] = " "
        end

        pos = state.data
        closest_charindex = argmin(
            [sum((pos .- center(bb)).^2) for bb in displayed_charbbs[]]
        )
        # set cursor to index of closest char if right of center, or previous char if left of center
        cursorindex[] = if (pos .- center(displayed_charbbs[][closest_charindex]))[1] > 0
            closest_charindex
        else
            closest_charindex - 1
        end
    end

    onmouseover(mousestate) do state
        hovering[] = true
    end

    onmouseout(mousestate) do state
        hovering[] = false
    end

    onmousedownoutside(mousestate) do state
        if reset_on_defocus[]
            reset_to_stored()
        end
        defocus!(ltextbox)
    end

    function insertchar!(c, index)
        if displayed_chars[] == [' ']
            empty!(displayed_chars[])
            index = 1
        end
        newchars = [displayed_chars[][1:index-1]; c; displayed_chars[][index:end]]
        displayed_string[] = join(newchars)
        cursorindex[] = index
    end

    function appendchar!(c)
        insertchar!(c, length(displayed_string[]))
    end

    function removechar!(index)
        newchars = [displayed_chars[][1:index-1]; displayed_chars[][index+1:end]]

        if isempty(newchars)
            newchars = [' ']
        end

        if cursorindex[] >= index
            cursorindex[] = max(0, cursorindex[] - 1)
        end

        displayed_string[] = join(newchars)
    end

    on(events(scene).unicode_input) do char_array
        if !focused[] || isempty(char_array)
            return
        end

        for c in char_array
            if is_allowed(c, restriction[])
                insertchar!(c, cursorindex[] + 1)
            end
        end
    end

    function submit()
        if displayed_is_valid[]
            stored_string[] = displayed_string[]
        end
    end

    function reset_to_stored()
        cursorindex[] = 0
        if isnothing(stored_string[])
            displayed_string[] = placeholder[]
        else
            displayed_string[] = stored_string[]
        end
    end

    function cursor_forward()
        cursorindex[] = min(length(displayed_string[]), cursorindex[] + 1)
    end

    function cursor_backward()
        cursorindex[] = max(0, cursorindex[] - 1)
    end

    on(events(scene).keyboardbuttons) do button_set
        if !focused[] || isempty(button_set)
            return
        end

        for key in button_set
            if key == Keyboard.backspace
                removechar!(cursorindex[])
            elseif key == Keyboard.delete
                removechar!(cursorindex[] + 1)
            elseif key == Keyboard.enter
                submit()
                if defocus_on_submit[]
                    defocus!(ltextbox)
                end
            elseif key == Keyboard.escape
                if reset_on_defocus[]
                    reset_to_stored()
                end
                defocus!(ltextbox)
            elseif key == Keyboard.right
                cursor_forward()
            elseif key == Keyboard.left
                cursor_backward()
            end
        end
    end

    ltextbox
end

function charbbs(text)
    positions = AbstractPlotting.layout_text(text[1][], text.position[], text.textsize[],
        text.font[], text.align[], text.rotation[], text.model[],
        text.justification[], text.lineheight[])

    font = AbstractPlotting.to_font(text.font[])

    bbs = [
        AbstractPlotting.FreeTypeAbstraction.height_insensitive_boundingbox(
            AbstractPlotting.FreeTypeAbstraction.get_extent(font, char),
            font
        )
        for char in text[1][]
    ]

    bbs_shifted_scaled = [Rect2D(
            (AbstractPlotting.origin(bb) .* text.textsize[] .+ pos[1:2])...,
            (widths(bb) .* text.textsize[])...)
        for (bb, pos) in zip(bbs, positions)]
end

function validate_textbox(str, validator::Function)
    validator(str)
end

function validate_textbox(str, T::Type)
    !isnothing(tryparse(T, str))
end

function validate_textbox(str, validator::Regex)
    m = match(validator, str)
    # check that the validator matches the whole string
    !isnothing(m) && m.match == str
end

function is_allowed(char, restriction::Nothing)
    true
end

function is_allowed(char, restriction::Function)
    allowed::Bool = restriction(char)
end

"""
    reset!(tb::Textbox)
Resets the stored_string of the given `Textbox` to `nothing` without triggering listeners, and resets the `Textbox` to the `placeholder` text.
"""
function reset!(tb::Textbox)
    tb.stored_string.val = nothing
    tb.displayed_string = tb.placeholder[]
    defocus!(tb)
    nothing
end

"""
    set!(tb::Textbox, string::String)
Sets the stored_string of the given `Textbox` to `string`, triggering listeners of `tb.stored_string`.
"""
function set!(tb::Textbox, string::String)
    if !validate_textbox(string, tb.validator[])
        error("Invalid string \"$(string)\" for textbox.")
    end

    tb.displayed_string = string
    tb.stored_string = string
    nothing
end

"""
    focus!(tb::Textbox)
Focuses an `Textbox` and makes it ready to receive keyboard input.
"""
function focus!(tb::Textbox)
    if !tb.focused[]
        tb.focused = true

        cursoranim = Animations.Loop(
            Animations.Animation(
                [0, 1.0],
                [Colors.alphacolor(COLOR_ACCENT[], 0), Colors.alphacolor(COLOR_ACCENT[], 1)],
                Animations.sineio(n = 2, yoyo = true, postwait = 0.2)),
                0.0, 0.0, 1000)

        if !isnothing(tb.cursoranimtask)
            Animations.stop(tb.cursoranimtask)
            tb.cursoranimtask = nothing
        end

        tb.cursoranimtask = Animations.animate_async(cursoranim; fps = 30) do t, color
            tb.cursorcolor = color
        end
    end
    nothing
end

"""
    defocus!(tb::Textbox)
Defocuses a `Textbox` so it doesn't receive keyboard input.
"""
function defocus!(tb::Textbox)

    if tb.displayed_string[] in (" ", "")
        tb.displayed_string[] = tb.placeholder[]
    end

    if !isnothing(tb.cursoranimtask)
        Animations.stop(tb.cursoranimtask)
        tb.cursoranimtask = nothing
    end
    tb.cursorcolor = :transparent
    tb.focused = false
    nothing
end


text:
"""
    PointVortex(x,y,q)

Construct a point vortex of charge `q` at `x`,`y`.

See also: [`rawData`](@ref), [`ScalarVortex`](@ref), [`vortex!`](@ref), [`randPointVortex`](@ref)
"""
function PointVortex(v::Array{Float64,2})
    ~isempty(v) ? (return PointVortex.(v[:,1],v[:,2],v[:,3])) : (return Array{PointVortex}(undef,0))
end
function PointVortex(vort::Array{ScalarVortex{T},1}) where T
    pv = PointVortex[]
    for j in eachindex(vort)
        push!(pv,vort[j].vort)
    end
    return pv
end

"""
    varray = rawData(v::PointVortex)

Return `varray<:Array{Float64,2}` of vortex data.

See also: [`PointVortex`](@ref)
"""
rawData(v::PointVortex) = [v.xv v.yv v.qv]
function rawData(v::Array{PointVortex,1})
    ~isempty(v) > 0 ? (return reduce(vcat,rawData.(v))) : (return Array{Float64}(undef, 0, 0))
end

function uniform(a,b)
    @assert a<b
    return a + (b-a)*rand()
end
function uniform(a,b,n)
    return uniform.(a*ones(n),b)
end
uniform() = uniform(-.5,.5)
uniform(n) = uniform(-.5,.5,n)

randcharge() = rand([-1 1],)
randcharge(n) = rand([-1 1],n)

"""
    randPointVortex()
    randPointVortex(n::Int)
    randPointVortex(ψ::Field)
    randPointVortex(n::Int,ψ::Field)

Sample `n::Int` random point vortices on spatial domain specified by `ψ::Field`.

See also: [`PointVortex`](@ref), [`randVortex`](@ref)
"""
randPointVortex() = PointVortex(uniform(),uniform(),randcharge())
randPointVortex(n) = PointVortex.(uniform(n), uniform(n), randcharge(n))
function randPointVortex(psi::Field)
    @unpack x,y = psi; a,b = first(x),last(x); c,d = first(y), last(y)
    return PointVortex(uniform(a,b,),uniform(c,d,),randcharge())
end
function randPointVortex(n,psi::Field)
    @unpack x,y = psi;
    dx = x[2]-x[1]; dy = y[2] - y[1]
    a,b = first(x)+2dx,last(x)-2dx; c,d = first(y)+2dy, last(y)-2dy
    return PointVortex.(uniform(a,b,n),uniform(c,d,n),randcharge(n))
end


text:
# Make sure these build and run without errors
@testset "Testing Asap4 From Sim Constructor" begin
    input_file = joinpath(@__DIR__, "..", "taskgraphs", "profile.json")
    c = AsapMapper.SimConstructor(input_file, (architecture = asap4,))
    m = AsapMapper.build_map(c)
    m = AsapMapper.asap_pnr(m)
    AsapMapper.report_routing_stats(m)
end

@testset "Testing ASAP 3" begin
    # Build Architecture
    input_file = joinpath(@__DIR__, "mapper_in.json")
    c = AsapMapper.PMConstructor(input_file)
    m = AsapMapper.build_map(c)
    # Placement
    m = AsapMapper.asap_pnr(m)
    AsapMapper.report_routing_stats(m)

    # Try plotting
    plt = plot(m)
    @test !isnothing(plt)

    AsapMapper.dump_map(m, "mapper_out.json")
    rm("mapper_out.json")
end

# Make sure these build and run without errors
@testset "Testing ASAP 3 profiled" begin
    # Build Architecture
    input_file = joinpath(@__DIR__, "mapper_in.json")
    c = AsapMapper.PMConstructor(input_file, (use_profiled_links = true,))
    m = AsapMapper.build_map(c)
    # Placement
    m = AsapMapper.asap_pnr(m)
    AsapMapper.report_routing_stats(m)

    # Test output generation
    AsapMapper.dump_map(m, "mapper_out.json")
    rm("mapper_out.json")
end

@testset "Testing Asap3 Hexagonal" begin
    input_file = joinpath(@__DIR__, "mapper_in.json")
    c = AsapMapper.PMConstructor(
        input_file,
        (architecture = () -> asap3(Hexagonal(2, 1)),)
    )

    m = AsapMapper.build_map(c)
    m = AsapMapper.asap_pnr(m)
    AsapMapper.report_routing_stats(m)
end

@testset "Testing Asap3 Hexagonal Profiled" begin
    input_file = joinpath(@__DIR__, "mapper_in.json")
    c = AsapMapper.PMConstructor(
        input_file,
        (architecture = () -> asap3(Hexagonal(2, 1)),
         use_profiled_links = true)
    )

    m = AsapMapper.build_map(c)
    m = AsapMapper.asap_pnr(m)
    AsapMapper.report_routing_stats(m)
end


text:
function basis_rels(b::Array{nf_elem, 1}, c; bd::fmpz = fmpz(10^35), no_p::Int = 4, no_rel::Int = 10000, no_coeff::Int = 4, no_id::fmpz = fmpz(0) )
  a = b[1].parent()
  t = b[1].parent()
  nb = length(b)
  one = fmpz(1)
  for i=1:no_rel
    zero!(a)
    for j=1:no_coeff
      cf = rand([-1, 1])
      p  = rand(1:nb)
      if cf==1
        Nemo.add!(a, a, b[p])
      else
        Nemo.sub!(a, a, b[p])
      end
    end
    iszero(a) && continue
    n = norm_div(a, one, no_p)
    if cmpabs(num(n), bd) <= 0
      if no_id != 0
        g = gcd(no_id, num(n))
        if g==1 || gcd(div(num(n), g), g) == 1
          if class_group_add_relation(c, a, n, one)
            a = b[1].parent()
          end
        end
      else
        if class_group_add_relation(c, a, n, one)
          a = b[1].parent()
        end
      end
    end
  end
end

function basis_rels_2(b::Array{nf_elem, 1}, bd::fmpz = fmpz(10^35), no_p::Int = 4, no_rel::Int = 10000, no_coeff::Int = 4, no_id::fmpz = fmpz(0), smooth = 0 )
  a = b[1].parent()
  t = b[1].parent()
  nb = length(b)
  one = fmpz(1)
  rels = Array(Hecke.nf_elem, no_rel)
  i = 1
  l = 0
  while i < no_rel + 1
    l = l + 1
    #println(l)
    if l % 1000 == 0
      println("so far $l tries")
    end
    zero!(a)
    for j=1:no_coeff
      cf = rand([-1, 1])
      p  = rand(1:nb)
      if cf==1
        Nemo.add!(a, a, b[p])
      else
        Hecke.sub!(a, a, b[p])
      end
    end
    iszero(a) && continue
    if no_id != 0
      n = norm_div(a, no_id, no_p)
    else
      n = norm_div(a, one, no_p)
    end
    if smooth != 0
      !is_smooth(smooth, num(n))[1] && continue
    end
    if cmpabs(num(n), bd) <= 0
      if no_id != 0
        g = gcd(no_id, num(n))
        if g==1 || gcd(div(num(n), g), g) == 1
          rels[i] = deepcopy(a)
          i = i +1
          println(i)
        end
      else
        rels[i] = deepcopy(a)
        i = i + 1
        println(i)
      end
    end
  end
  return rels
end

function improve(c::Hecke.ClassGrpCtx)
  H = sub(c.M, 1:rows(c.M), 1:cols(c.M))
  Hecke.upper_triangular(H, mod = 17)
  p = setdiff(Set(1:cols(H)), Set([x.pos[1] for x=H.rows]))
  p = maximum(p)
  b = Hecke.bkz_basis(c.FB.ideals[p]);
#  b = rels_stat(b, ...)
  for x=b
    class_group_add_relation(c, b, n, one)
  end
end

function rels_stat(b::Array{Hecke.nf_elem, 1}; no_p = 4, no_rel::Int = 10000, no_coeff::Int = 4, fixed = 0, smooth=0 )
  a = b[1].parent()
  t = b[1].parent()
  nb = length(b)
  one = fmpz(1)

  stat = Dict{Int, Int}()
  if fixed != 0
    stat[-1] = 0
  end
  if smooth != 0
    stat[-2] = 0
  end
  all_g = Array{Any, 1}()
  for i=1:no_rel
    zero!(a)
    for j=1:no_coeff
      cf = rand([-1, 1])
      p  = rand(1:nb)
      if cf==1
        Nemo.add!(a, a, b[p])
      else
        Nemo.sub!(a, a, b[p])
      end
    end
    if iszero(a) continue; end;
    if fixed!=0 && fixed(a) == a
      stat[-1] += 1
    end
    n = norm_div(a, one, no_p)
    k = div(ndigits(num(n), 2), 5)
    if haskey(stat, k)
      stat[k] += 1
    else
      stat[k] = 1
    end
    if smooth != 0 && is_smooth(smooth, num(n))[1]
      stat[-2] += 1
      push!(all_g, a)
      a =   a = b[1].parent()
    end
  end
  return stat, all_g
end

function find_rels(b::Array{Hecke.nf_elem, 1}; no_p = 4, no_rel::Int = 10000, no_coeff::Int = 4, fixed = 0, smooth=0 )

  for i=10:50
    st, re = rels_stat(b, no_p = no_p, no_rel = no_rel, no_coeff = no_coeff, smooth =smooth)
    toNemo("/home/fieker/Rels128/rels128.$i", re, name="R$i");
  end
end

function find_rels2(b::Array{Hecke.nf_elem, 1}; no_p = 4, no_rel::Int = 10000, no_coeff::Int = 4, fixed = 0, smooth=0 )

  for i=100:150
    st, re = rels_stat(b, no_p = no_p, no_rel = no_rel, no_coeff = no_coeff, smooth =smooth)
    toNemo("/home/fieker/Rels128/rels128.$i", re, name="R$i");
  end
end

function int_fb_max_real(f::Int, B::Int)
  lp = filter(isprime, 2:B)
  l1 = filter(x -> (x % f) ==1 || (x %f ) == f-1, lp) # the deg 1 primes
  filter!(x -> !(x in l1), lp)
  lr = filter(x -> f % x == 0, lp)  # ramified primes
  filter!(x -> !(x in lr), lp)
  lu = Array{Int, 1}()

  d = euler_phi(f)
  for r =2:d
    if d % r != 0
      continue
    end
    for p = lp
      if fmpz(p)^r  > B
        break
      end
      if powermod(p, r, f)==1 || powermod(p, r, f) == f-1
        push!(lu, p)
      end
    end
    filter!(x -> !(x in lu), lp)
  end

  return l1, lr, lu
end

function statistic{T}(st::Dict{Int,T})
  s = T(0)
  s2 = T(0)
  n = 0
  for k in keys(st)
    if k<1
      continue
    end
    t = st[k]
    s += k*t
    s2 += k^2*t
    n += t
  end
  av = s/n
  si = s2/n-av^2

  return av, si
end

function statistic{T}(st::Array{T, 1})
  s = T(0)
  s2 = T(0)
  n = length(st)
  for t in st
    s += t
    s2 += t^2
  end
  av = s/n
  si = s2/n-av^2

  return av, si
end

# If M is upper-triangular with more columns then rows,
# this function returns the non-pivot column indices.
function _find_missing_pivot(M::Smat)
  return setdiff(Set(1:cols(M)), Set([y.pos[1] for y = M.rows ]))
end

function res_degree_in_max_real(p::Int, n::Int)
  r = valuation(n, 2)[2]
  r != 1 && error("Not yet implemented")
  f = Hecke.modord(p, n)
  return (powermod(p, div(f,2), n) == n-1) ? (return div(f, 2)) : (return f)
end

#function elem_type(::AnticNumberField)
#  return Nemo.nf_elem
#end


text:
using Langevin
using ProgressMeter
using Plots

D = Dict{Symbol,Any}()
begin
    D[:matter] = "radiation"
    D[:case] = "lattice"
    D[:space] = "2d1f1v"
    D[:flux] = "kfvs"
    D[:collision] = "bgk"
    D[:nSpecies] = 1
    D[:interpOrder] = 1
    D[:limiter] = "vanleer"
    D[:boundary] = "fix"
    D[:cfl] = 0.5
    D[:maxTime] = 0.3

    D[:x0] = -3.5
    D[:x1] = 3.5
    D[:nx] = 50
    D[:y0] = -3.5
    D[:y1] = 3.5
    D[:ny] = 50
    D[:pMeshType] = "uniform"
    D[:nxg] = 0
    D[:nyg] = 0

    D[:quadrature] = "legendre"
    D[:nq] = 8

    D[:knudsen] = 1.0
    D[:sigmaS] = 1.0
    D[:sigmaA] = 0.0
end

set = set_setup(D)
ps = set_geometry(D)
points, weights = legendre_quadrature(D[:nq])
vs = UnstructVSpace(-1.0, 1.0, length(weights), points, weights)
radiation = Radiation(D[:knudsen], D[:sigmaS], D[:sigmaA])
ib = nothing

ks = SolverSet(set, ps, vs, radiation, ib, pwd())

uq = UQ1D(3, 6, -0.1, 0.1)

function is_absorb(x::T, y::T) where {T<:Real}
    cds = Array{Bool}(undef, 11) # conditions

    cds[1] = -2.5 < x < -1.5 && 1.5 < y < 2.5
    cds[2] = -2.5 < x < -1.5 && -0.5 < y < 0.5
    cds[3] = -2.5 < x < -1.5 && -2.5 < y < -1.5
    cds[4] = -1.5 < x < -0.5 && 0.5 < y < 1.5
    cds[5] = -1.5 < x < -0.5 && -1.5 < y < -0.5
    cds[6] = -0.5 < x < 0.5 && -2.5 < y < -1.5
    cds[7] = 0.5 < x < 1.5 && 0.5 < y < 1.5
    cds[8] = 0.5 < x < 1.5 && -1.5 < y < -0.5
    cds[9] = 1.5 < x < 2.5 && 1.5 < y < 2.5
    cds[10] = 1.5 < x < 2.5 && -0.5 < y < 0.5
    cds[11] = 1.5 < x < 2.5 && -2.5 < y < -1.5

    if any(cds) == true
        return true
    else
        return false
    end
end

begin
    σs = zeros(Float64, ps.nx, ps.ny)
    σa = zeros(Float64, ps.nx, ps.ny)
    for i = 1:ps.nx, j = 1:ps.ny
        if is_absorb(ps.x[i, j], ps.y[i, j])
            σs[i, j] = 0.0
            σa[i, j] = 10.0
        else
            σs[i, j] = 1.0
            σa[i, j] = 0.0
        end
    end
    σt = σs + σa
    σq = zeros(Float64, ps.nx, ps.ny)
    for i = 1:ps.nx, j = 1:ps.ny
        if -0.5 < ps.x[i, j] < 0.5 && -0.5 < ps.y[i, j] < 0.5
            σq[i, j] = 30.0 / (4.0 * π)
        else
            σq[i, j] = 0.0
        end
    end
end

phi = zeros(vs.nu, ps.nx, ps.ny)
for j = 1:ps.nx
    for i = 1:ps.ny
        phi[:, i, j] .= 1e-4
    end
end

t = 0.0
dt = 1.2 / 150 * ks.set.cfl
flux1 = zeros(vs.nu, ps.nx + 1, ps.ny)
flux2 = zeros(vs.nu, ps.nx, ps.ny + 1)

@showprogress for iter = 1:200
    for i = 2:ps.nx, j = 1:ps.ny
        tmp = @view flux1[:, i, j]
        flux_kfvs!(tmp, phi[:, i-1, j], phi[:, i, j], points[:, 1], dt)
    end
    for i = 1:ps.nx, j = 2:ps.ny
        tmp = @view flux2[:, i, j]
        flux_kfvs!(tmp, phi[:, i, j-1], phi[:, i, j], points[:, 2], dt)
    end

    for j = 1:ps.ny, i = 1:ps.nx
        integral = discrete_moments(phi[:, i, j], weights)
        integral /= 4.0 * π

        for q = 1:vs.nu
            phi[q, i, j] =
                phi[q, i, j] +
                (flux1[q, i, j] - flux1[q, i+1, j]) / ps.dx[i, j] +
                (flux2[q, i, j] - flux2[q, i, j+1]) / ps.dy[i, j] +
                dt * σs[i, j] * integral - dt * σt[i, j] * phi[q, i, j] + σq[i, j] * dt
        end
    end

    global t += dt
end

ρ = zeros(ps.nx, ps.ny)
for i = 1:ps.nx, j = 1:ps.ny
    ρ[i, j] = discrete_moments(phi[:, i, j], weights)
end
contourf(ps.x[1:end, 1], ps.y[1, 1:end], ρ', color = :PiYG_3)


text:
module Extension
using MLStyle.Err

Support = Dict{Symbol, Vector{Module}}(
    :GADT      => []
)

export use, @use, used

function used(ext :: Symbol, mod :: Module) :: Bool
    get(Support, ext) do
        throw(UnknownExtension(ext))
    end |> mods -> mod in mods
end

function use(ext :: Symbol, mod :: Module)
    mods = get(Support, ext) do
        Support[ext] = []
    end
    push!(mods, mod)
end

macro use(exts...)
    mod = __module__
    for each in exts
        use(each, mod)
    end
end

end


text:
#=

Example demonstrating sparsity control by means of the "K"-factor passed to the lazy AFW variant

A larger K >= 1 favors sparsity by favoring optimization over the current active set rather than
adding a new FW vertex.

The default for AFW is K = 2.0

=#

using FrankWolfe
using ProgressMeter
using Arpack
using Plots
using DoubleFloats
using ReverseDiff

using LinearAlgebra
using Random

n = Int(1e4)
k = 1000

s = rand(1:100)
@info "Seed $s"
Random.seed!(s)

xpi = rand(n);
total = sum(xpi);

# here the optimal solution lies in the interior if you want an optimal solution on a face and not the interior use:
# const xp = xpi;

const xp = xpi ./ total;

f(x) = norm(x - xp)^2
function grad!(storage, x)
    @. storage = 2 * (x - xp)
end

const lmo = FrankWolfe.KSparseLMO(100, 1.0)

## other LMOs to try
# lmo_big = FrankWolfe.KSparseLMO(100, big"1.0")
# lmo = FrankWolfe.LpNormLMO{Float64,5}(1.0)
# lmo = FrankWolfe.ProbabilitySimplexOracle(1.0);
# lmo = FrankWolfe.UnitSimplexOracle(1.0);

const x00 = FrankWolfe.compute_extreme_point(lmo, rand(n))

## example with BirkhoffPolytopeLMO - uses square matrix.
# const lmo = FrankWolfe.BirkhoffPolytopeLMO()
# cost = rand(n, n)
# const x00 = FrankWolfe.compute_extreme_point(lmo, cost)

function build_callback(trajectory_arr)
    return function callback(state)
        return push!(trajectory_arr, (Tuple(state)[1:5]..., length(state.active_set)))
    end
end

FrankWolfe.benchmark_oracles(f, grad!, () -> randn(n), lmo; k=100)

x0 = deepcopy(x00)
@time x, v, primal, dual_gap, trajectory_shortstep = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
);

trajectory_adaptive = []
callback = build_callback(trajectory_adaptive)

x0 = deepcopy(x00)
@time x, v, primal, dual_gap, _ = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
    callback=callback,
);

println("\n==> Lazy AFW.\n")

trajectory_adaptiveLoc15 = []
callback = build_callback(trajectory_adaptiveLoc15)

x0 = deepcopy(x00)
@time x, v, primal, dual_gap, _ = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    lazy=true,
    K=1.5,
    trajectory=true,
    callback=callback,
);

trajectory_adaptiveLoc2 = []
callback = build_callback(trajectory_adaptiveLoc2)

x0 = deepcopy(x00)
@time x, v, primal, dual_gap, trajectory_adaptiveLoc2 = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    lazy=true,
    K=2.0,
    trajectory=true,
    callback=callback,
);

trajectory_adaptiveLoc4 = []
callback = build_callback(trajectory_adaptiveLoc4)

x0 = deepcopy(x00)
@time x, v, primal, dual_gap, _ = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    K=4.0,
    lazy=true,
    trajectory=true,
    callback=callback,
);

trajectory_adaptiveLoc10 = []
callback = build_callback(trajectory_adaptiveLoc10)

x0 = deepcopy(x00)
@time x, v, primal, dual_gap, _ = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    lazy=true,
    K=10.0,
    verbose=true,
    trajectory=true,
    callback=callback,
);

# Reduction primal/dual error vs. sparsity of solution

dataSparsity =
    [trajectory_adaptive, trajectory_adaptiveLoc15, trajectory_adaptiveLoc2, trajectory_adaptiveLoc4, trajectory_adaptiveLoc10]
labelSparsity = ["AFW", "LAFW-K-1.5", "LAFW-K-2.0", "LAFW-K-4.0", "LAFW-K-10.0"]

FrankWolfe.plot_sparsity(dataSparsity, labelSparsity, legend_position=:topright)


text:
#using Requests
#ecoliModel = Requests.get("http://gcrg.ucsd.edu/sites/default/files/Attachments/Images/downloads/Ecoli_core/ecoli_core_model.mat")
#save(ecoliModel, "ecoli_core_model.mat")

n_cobra = MetabolicFluxAnalysis.import_cobra("testdata/ecoli_core_model.mat","model")
#(myN ,myN_p , myN_n) = getN( n_cobra )
(myN ,myN_p , myN_n) = MetabolicFluxAnalysis.get_N( n_cobra ; sparse_n=true )

# benchmark stoich. matrix generation..:
sns = Array{Any,1}(10_000)
tic()
for zi=1:10_000
    ( sns_i , snxx ,  snyy ) = MetabolicFluxAnalysis.get_N(n_cobra, sparse_n=rand()>0.5 )
    sns[zi] = sns_i
end
toc()