text:
""" 2-category of finitely presented categories.

This module is for the 2-category **Cat** what the module [`FinSets](@ref) is
for the category **Set**: a finitary, combinatorial setting where explicit
calculations can be carried out. We emphasize that the prefix `Fin` means
"finitely presented," not "finite," as finite categories are too restrictive a
notion for many purposes. For example, the free category on a graph is finite if
and only if the graph is DAG, which is a fairly special condition. This usage of
`Fin` is also consistent with `FinSet` because for sets, being finite and being
finitely presented are equivalent.
"""
module FinCats
export FinCat, FinCatGraph, Path, ob_generators, hom_generators, equations,
  is_discrete, is_free, graph, edges, src, tgt, presentation,
  FinFunctor, FinDomFunctor, is_functorial, collect_ob, collect_hom, force,
  FinTransformation, components, is_natural, is_initial

using AutoHashEquals
using Reexport
using StaticArrays: SVector
using DataStructures: IntDisjointSets, in_same_set, num_groups

@reexport using ..Categories
using ...GAT, ...Present, ...Syntax
import ...Present: equations
using ...Theories: Category, Schema, ObExpr, HomExpr, AttrExpr, AttrTypeExpr
import ...Theories: dom, codom, id, compose, ⋅, ∘
using ...CSetDataStructures, ...Graphs
import ...Graphs: edges, src, tgt, enumerate_paths
import ..Categories: ob, hom, ob_map, hom_map, component

# Categories
############

""" Abstract type for finitely presented category.
"""
abstract type FinCat{Ob,Hom} <: Cat{Ob,Hom} end

FinCat(g::HasGraph, args...; kw...) = FinCatGraph(g, args...; kw...)
FinCat(pres::Presentation, args...; kw...) =
  FinCatPresentation(pres, args...; kw...)

""" Object generators of finitely presented category.

The object generators are almost always the same as the objects. In principle,
however, it is possible to have equations between objects, so that there are
fewer objects than object generators.
"""
function ob_generators end

""" Morphism generators of finitely presented category.
"""
function hom_generators end

""" Is the category discrete?

A category is *discrete* if it is has no non-identity morphisms.
"""
is_discrete(C::FinCat) = isempty(hom_generators(C))

""" Is the category freely generated?
"""
is_free(C::FinCat) = isempty(equations(C))

# Categories on graphs
######################

""" Abstract type for category backed by finite generating graph.
"""
abstract type FinCatGraph{G,Ob,Hom} <: FinCat{Ob,Hom} end

""" Generating graph for a finitely presented category.
"""
graph(C::FinCatGraph) = C.graph

ob_generators(C::FinCatGraph) = vertices(graph(C))
hom_generators(C::FinCatGraph) = edges(graph(C))

function Base.show(io::IO, C::FinCatGraph)
  print(io, "FinCat(")
  show(io, graph(C))
  print(io, ", [")
  join(io, equations(C), ", ")
  print(io, "])")
end

# Paths in graphs
#----------------

""" Path in a graph.

The path is allowed to be empty but always has definite start and end points
(source and target vertices).
"""
@auto_hash_equals struct Path{V,E,Edges<:AbstractVector{E}}
  edges::Edges
  src::V
  tgt::V
end
edges(path::Path) = path.edges
src(path::Path) = path.src
tgt(path::Path) = path.tgt

function Path(g::HasGraph, es::AbstractVector)
  !isempty(es) || error("Nonempty edge list needed for nontrivial path")
  all(e -> has_edge(g, e), es) || error("Path has edges not contained in graph")
  Path(es, src(g, first(es)), tgt(g, last(es)))
end

function Path(g::HasGraph, e)
  has_edge(g, e) || error("Edge $e not contained in graph $g")
  Path(SVector(e), src(g,e), tgt(g,e))
end

function Base.empty(::Type{Path}, g::HasGraph, v::T) where T
  has_vertex(g, v) || error("Vertex $v not contained in graph $g")
  Path(SVector{0,T}(), v, v)
end

function Base.vcat(p1::Path, p2::Path)
  tgt(p1) == src(p2) ||
    error("Path start/end points do not match: $(tgt(p1)) != $(src(p2))")
  Path(vcat(edges(p1), edges(p2)), src(p1), tgt(p2))
end

function Base.show(io::IO, path::Path)
  print(io, "Path(")
  show(io, edges(path))
  print(io, ": $(src(path)) → $(tgt(path)))")
end

""" Abstract type for category whose morphisms are paths in a graph.

(Or equivalence classes of paths in a graph, but we compute with
"""
const FinCatPathGraph{G,V,E} = FinCatGraph{G,V,Path{V,E}}

dom(C::FinCatPathGraph, e) = src(graph(C), e)
dom(C::FinCatPathGraph, path::Path) = src(path)
codom(C::FinCatPathGraph, e) = tgt(graph(C), e)
codom(C::FinCatPathGraph, path::Path) = tgt(path)

id(C::FinCatPathGraph, x) = empty(Path, graph(C), x)
compose(C::FinCatPathGraph, fs...) =
  reduce(vcat, coerce_path(graph(C), f) for f in fs)

ob(C::FinCatPathGraph, x) = has_vertex(graph(C), x) ? x :
  error("Vertex $x not contained in graph $(graph(C))")
hom(C::FinCatPathGraph, f) = coerce_path(graph(C), f)

coerce_path(g::HasGraph, path::Path) = path
coerce_path(g::HasGraph, x) = Path(g, x)

# Free category on graph
#-----------------------

""" Free category generated by a finite graph.

The objects of the free category are vertices in the graph and the morphisms are
(possibly empty) paths.
"""
@auto_hash_equals struct FreeCatGraph{G<:HasGraph} <: FinCatPathGraph{G,Int,Int}
  graph::G
end

FinCatGraph(g::HasGraph) = FreeCatGraph(g)

is_free(::FreeCatGraph) = true

function Base.show(io::IO, C::FreeCatGraph)
  print(io, "FinCat(")
  show(io, graph(C))
  print(io, ")")
end

# Category on graph with equations
#---------------------------------

""" Category presented by a finite graph together with path equations.

The objects of the category are vertices in the graph and the morphisms are
paths, quotiented by the congruence relation generated by the path equations.
See (Spivak, 2014, *Category theory for the sciences*, §4.5).
"""
@auto_hash_equals struct FinCatGraphEq{G<:HasGraph,Eqs<:AbstractVector{<:Pair}} <:
    FinCatPathGraph{G,Int,Int}
  graph::G
  equations::Eqs
end

equations(C::FinCatGraphEq) = C.equations

function FinCatGraph(g::HasGraph, eqs::AbstractVector)
  eqs = map(eqs) do eq
    lhs, rhs = coerce_path(g, first(eq)), coerce_path(g, last(eq))
    (src(lhs) == src(rhs) && tgt(lhs) == tgt(rhs)) ||
      error("Paths $lhs and $rhs in equation do not have equal (co)domains")
    lhs => rhs
  end
  FinCatGraphEq(g, eqs)
end

# Symbolic categories
#####################

""" Category defined by a `Presentation` object.

The presentation type can, of course, be a category (`Theories.Category`). It
can also be a schema (`Theories.Schema`). In this case, the schema's objects and
attribute types are regarded as the category's objects and the schema's
morphisms, attributes, and attribute types as the category's morphisms (where
the attribute types are identity morphisms). When the schema is formalized as a
profunctor whose codomain category is discrete, this amounts to taking the
collage of the profunctor.
"""
@auto_hash_equals struct FinCatPresentation{T,Ob,Hom} <: FinCat{Ob,Hom}
  presentation::Presentation{T}
end

function FinCatPresentation(pres::Presentation{T}) where T
  S = pres.syntax
  FinCatPresentation{T,S.Ob,S.Hom}(pres)
end
function FinCatPresentation(pres::Presentation{Schema})
  S = pres.syntax
  Ob = Union{S.Ob, S.AttrType}
  Hom = Union{S.Hom, S.Attr, S.AttrType}
  FinCatPresentation{Schema,Ob,Hom}(pres)
end

presentation(C::FinCatPresentation) = C.presentation

ob_generators(C::FinCatPresentation) = generators(presentation(C), :Ob)
ob_generators(C::FinCatPresentation{Schema}) = let P = presentation(C)
  vcat(generators(P, :Ob), generators(P, :AttrType))
end

hom_generators(C::FinCatPresentation) = generators(presentation(C), :Hom)
hom_generators(C::FinCatPresentation{Schema}) = let P = presentation(C)
  vcat(generators(P, :Hom), generators(P, :Attr))
end

equations(C::FinCatPresentation) = equations(presentation(C))

ob(C::FinCatPresentation, x) = ob(C, presentation(C)[x])
ob(C::FinCatPresentation, x::GATExpr) =
  gat_typeof(x) == :Ob ? x : error("Expression $x is not an object")
ob(C::FinCatPresentation{Schema}, x::GATExpr) =
  gat_typeof(x) ∈ (:Ob, :AttrType) ? x :
    error("Expression $x is not an object or attribute type")

hom(C::FinCatPresentation, f) = hom(C, presentation(C)[f])
hom(C::FinCatPresentation, fs::AbstractVector) =
  mapreduce(f -> hom(C, f), compose, fs)
hom(C::FinCatPresentation, f::GATExpr) =
  gat_typeof(f) == :Hom ? f : error("Expression $f is not a morphism")
hom(C::FinCatPresentation{Schema}, f::GATExpr) =
  gat_typeof(f) ∈ (:Hom, :Attr, :AttrType) ? f :
    error("Expression $f is not a morphism or attribute")

id(C::FinCatPresentation{Schema}, x::AttrTypeExpr) = x
compose(C::FinCatPresentation{Schema}, f::AttrTypeExpr, g::AttrTypeExpr) =
  (f == g) ? f : error("Invalid composite of attribute type identities: $f != $g")

function Base.show(io::IO, C::FinCatPresentation)
  print(io, "FinCat(")
  show(io, presentation(C))
  print(io, ")")
end

# Functors
##########

""" A functor out of a finitely presented category.
"""
const FinDomFunctor{Dom<:FinCat,Codom<:Cat} = Functor{Dom,Codom}

FinDomFunctor(maps::NamedTuple{(:V,:E)}, dom::FinCatGraph, codom::Cat) =
  FinDomFunctor(maps.V, maps.E, dom, codom)

function FinDomFunctor(ob_map, dom::FinCat, codom::Cat{Ob,Hom}) where {Ob,Hom}
  is_discrete(dom) ||
    error("Morphism map omitted by domain category is not discrete: $dom")
  FinDomFunctor(ob_map, empty(ob_map, Hom), dom, codom)
end
FinDomFunctor(ob_map, ::Nothing, dom::FinCat, codom::Cat) =
  FinDomFunctor(ob_map, dom, codom)

function hom_map(F::FinDomFunctor{<:FinCatPathGraph}, path::Path)
  D = codom(F)
  mapreduce(e -> hom_map(F, e), (gs...) -> compose(D, gs...),
            edges(path), init=id(D, ob_map(F, src(path))))
end

ob_map(F::FinDomFunctor, x::GATExpr{:generator}) = ob_map(F, first(x))
hom_map(F::FinDomFunctor, f::GATExpr{:generator}) = hom_map(F, first(f))
hom_map(F::FinDomFunctor, f::GATExpr{:id}) = id(codom(F), ob_map(F, dom(f)))

function hom_map(F::FinDomFunctor, f::GATExpr{:compose})
  D = codom(F)
  mapreduce(f -> hom_map(F, f), (gs...) -> compose(D, gs...), args(f))
end

(F::FinDomFunctor)(expr::ObExpr) = ob_map(F, expr)
(F::FinDomFunctor)(expr::HomExpr) = hom_map(F, expr)

Categories.show_type_constructor(io::IO, ::Type{<:FinDomFunctor}) =
  print(io, "FinDomFunctor")

""" Collect assignments of functor's object map as a vector.
"""
collect_ob(F::FinDomFunctor) = map(x -> ob_map(F, x), ob_generators(dom(F)))

""" Collect assignments of functor's morphism map as a vector.
"""
collect_hom(F::FinDomFunctor) = map(f -> hom_map(F, f), hom_generators(dom(F)))

""" Is the purported functor on a presented category functorial?

This function checks that functor preserves domains and codomains. When
`check_equations` is `true` (the default), it also checks that the functor
preserves all equations in the domain category. Note that in some cases this may
not be possible.

See also: [`is_natural`](@ref).
"""
function is_functorial(F::FinDomFunctor; check_equations::Bool=false)
  C, D = dom(F), codom(F)
  all(hom_generators(C)) do f
    g = hom_map(F, f)
    dom(D,g) == ob_map(F, dom(C,f)) && codom(D,g) == ob_map(F, codom(C,f))
  end || return false

  if check_equations
    all(equations(C)) do (lhs, rhs)
      is_hom_equal(D, hom_map(F, lhs), hom_map(F, rhs))
    end || return false
  end

  true
end

""" A functor between finitely presented categories.
"""
const FinFunctor{Dom<:FinCat,Codom<:FinCat} = FinDomFunctor{Dom,Codom}

FinFunctor(maps, dom::FinCat, codom::FinCat) = FinDomFunctor(maps, dom, codom)
FinFunctor(ob_map, hom_map, dom::FinCat, codom::FinCat) =
  FinDomFunctor(ob_map, hom_map, dom, codom)
FinFunctor(ob_map, hom_map, dom::Presentation, codom::Presentation) =
  FinDomFunctor(ob_map, hom_map, FinCat(dom), FinCat(codom))

Categories.show_type_constructor(io::IO, ::Type{<:FinFunctor}) =
  print(io, "FinFunctor")

"""
Dual to a ["final-functor"](https://ncatlab.org/nlab/show/final+functor), an
initial functor is one for which pulling back diagrams along it does not change
the limits of these diagrams.

This amounts to checking, for a functor C->D, that, for every object d in
Ob(D), the comma category (F/d) is connected.
"""
function is_initial(F::FinFunctor)::Bool
  Gₛ, Gₜ = graph(dom(F)), graph(codom(F))
  pathₛ, pathₜ = enumerate_paths.([Gₛ, Gₜ])

  function connected_nonempty_slice(t::Int)::Bool
    paths_into_t = incident(pathₜ, t, :tgt)
    # Generate slice objects
    ob_slice = Pair{Int,Vector{Int}}[] # s ∈Ob(S) and a path ∈ T(F(s), t)
    for s in vertices(Gₛ)
      paths_s_to_t = incident(pathₜ, ob_map(F,s), :src) ∩ paths_into_t
      append!(ob_slice, [s => pathₜ[p, :eprops] for p in paths_s_to_t])
    end

    # Empty case
    if isempty(ob_slice)
      return false
    end

    """
    For two slice objects (m,pₘ) and (n,pₙ) check for a morphism f ∈ S(M,N) such
    that there is a commutative triangle pₘ = f;pₙ
    """
    function check_pair(i::Int, j::Int)::Bool
      (m,pₘ), (n,pₙ) = ob_slice[i], ob_slice[j]
      es = incident(pathₛ, m, :src) ∩ incident(pathₛ, n, :tgt)
      paths = pathₛ[es, :eprops]
      return any(f -> pₘ == vcat(edges.(hom_map(F,f))..., pₙ), paths)
    end

    # Use check_pair to determine pairwise connectivity
    connected = IntDisjointSets(length(ob_slice)) # sym/trans/refl closure
    obs = 1:length(ob_slice)
    for (i,j) in Base.Iterators.product(obs, obs)
      if !in_same_set(connected, i, j) && check_pair(i,j)
        union!(connected, i, j)
      end
    end

    return num_groups(connected) == 1
  end

  # Check for each t ∈ T whether F/t is connected
  return all(connected_nonempty_slice, 1:nv(Gₜ))
end

# Mapping-based functors
#-----------------------

""" Functor out of a finitely presented category given by explicit mappings.
"""
@auto_hash_equals struct FinDomFunctorMap{Dom<:FinCat, Codom<:Cat,
    ObMap, HomMap} <: FinDomFunctor{Dom,Codom}
  ob_map::ObMap
  hom_map::HomMap
  dom::Dom
  codom::Codom
end

function FinDomFunctor(ob_map::Union{AbstractVector,AbstractDict},
                       hom_map::Union{AbstractVector,AbstractDict},
                       dom::FinCat, codom::Cat)
  length(ob_map) == length(ob_generators(dom)) ||
    error("Length of object map $ob_map does not match domain $dom")
  length(hom_map) == length(hom_generators(dom)) ||
    error("Length of morphism map $hom_map does not match domain $dom")

  ob_map = mappairs(x -> functor_key(dom, x), y -> ob(codom, y), ob_map)
  hom_map = mappairs(f -> functor_key(dom, f), g -> hom(codom, g), hom_map)
  FinDomFunctorMap(ob_map, hom_map, dom, codom)
end
FinDomFunctor(ob_map::Union{AbstractVector{Ob},AbstractDict{<:Any,Ob}},
              hom_map::Union{AbstractVector{Hom},AbstractDict{<:Any,Hom}},
              dom::FinCat) where {Ob,Hom} =
  FinDomFunctor(ob_map, hom_map, dom, TypeCat(Ob, Hom))

functor_key(C::FinCat, x) = x
functor_key(C::FinCat, expr::GATExpr) = head(expr) == :generator ?
  first(expr) : error("Functor must be defined on generators")

Categories.do_ob_map(F::FinDomFunctorMap, x) = F.ob_map[x]
Categories.do_hom_map(F::FinDomFunctorMap, f) = F.hom_map[f]

collect_ob(F::FinDomFunctorMap) = values(F.ob_map)
collect_hom(F::FinDomFunctorMap) = values(F.hom_map)

""" Force evaluation of lazily defined function or functor.
"""
function force(F::FinDomFunctor)
  C = dom(F)
  FinDomFunctor(make_map(x -> ob_map(F, x), ob_generators(C)),
                make_map(f -> hom_map(F, f), hom_generators(C)), C)
end
force(F::FinDomFunctorMap) = F

function Categories.do_compose(F::FinDomFunctorMap, G::FinDomFunctorMap)
  FinDomFunctorMap(mapvals(x -> ob_map(G, x), F.ob_map),
                   mapvals(f -> hom_map(G, f), F.hom_map), dom(F), codom(G))
end

function Base.show(io::IO, F::T) where T <: FinDomFunctorMap
  Categories.show_type_constructor(io, T); print(io, "(")
  show(io, F.ob_map)
  print(io, ", ")
  show(io, F.hom_map)
  print(io, ", ")
  Categories.show_domains(io, F)
  print(io, ")")
end

# Natural transformations
#########################

""" A natural transformation whose domain category is finitely generated.

This type is for natural transformations ``α: F ⇒ G: C → D`` such that the
domain category ``C`` is finitely generated. Such a natural transformation is
given by a finite amount of data (one morphism in ``D`` for each generating
object of ``C``) and its naturality is verified by finitely many equations (one
equation for each generating morphism of ``C``).
"""
const FinTransformation{C<:FinCat,D<:Cat,Dom<:FinDomFunctor,Codom<:FinDomFunctor} =
  Transformation{C,D,Dom,Codom}

FinTransformation(F, G; components...) = FinTransformation(components, F, G)

""" Components of a natural transformation.
"""
components(α::FinTransformation) =
  make_map(x -> component(α, x), ob_generators(dom_ob(α)))

""" Is the transformation between `FinDomFunctors` a natural transformation?

This function uses the fact that to check whether a transformation is natural,
it suffices to check the naturality equations on a generating set of morphisms
of the domain category. In some cases, checking the equations may be expensive
or impossible. When the keyword argument `check_equations` is `false`, only the
domains and codomains of the components are checked.

See also: [`is_functorial`](@ref).
"""
function is_natural(α::FinTransformation; check_equations::Bool=true)
  F, G = dom(α), codom(α)
  C, D = dom(F), codom(F) # == dom(G), codom(G)
  all(ob_generators(C)) do c
    α_c = α[c]
    dom(D, α_c) == ob_map(F,c) && codom(D, α_c) == ob_map(G,c)
  end || return false

  if check_equations
    all(hom_generators(C)) do f
      Ff, Gf = hom_map(F,f), hom_map(G,f)
      α_c, α_d = α[dom(C,f)], α[codom(C,f)]
      is_hom_equal(D, compose(D, α_c, Gf), compose(D, Ff, α_d))
    end || return false
  end

  true
end

function check_transformation_domains(F::Functor, G::Functor)
  # XXX: Equality of TypeCats is too strict, so for now we are punting on
  # (co)domain checks in that case.
  C, C′ = dom(F), dom(G)
  (C isa TypeCat && C′ isa TypeCat) || C == C′ ||
    error("Mismatched domains in functors $F and $G")
  D, D′ = codom(F), codom(G)
  (D isa TypeCat && D′ isa TypeCat) || D == D′ ||
    error("Mismatched codomains in functors $F and $G")
  (C, D)
end

# Mapping-based transformations
#------------------------------

""" Natural transformation with components given by explicit mapping.
"""
@auto_hash_equals struct FinTransformationMap{C<:FinCat,D<:Cat,
    Dom<:FinDomFunctor{C,D},Codom<:FinDomFunctor,Comp} <: FinTransformation{C,D,Dom,Codom}
  components::Comp
  dom::Dom
  codom::Codom
end

function FinTransformation(components::Union{AbstractVector,AbstractDict},
                           F::FinDomFunctor, G::FinDomFunctor)
  C, D = check_transformation_domains(F, G)
  length(components) == length(ob_generators(C)) ||
    error("Incorrect number of components in $components for domain category $C")
  components = mappairs(x -> transformation_key(C,x), f -> hom(D,f), components)
  FinTransformationMap(components, F, G)
end

transformation_key(C::FinCat, x) = x
transformation_key(C::FinCat, expr::GATExpr) = head(expr) == :generator ?
  first(expr) : error("Natural transformation must be defined on generators")

component(α::FinTransformationMap, x) = α.components[x]
component(α::FinTransformationMap, expr::GATExpr{:generator}) =
  component(α, first(expr))

components(α::FinTransformationMap) = α.components

function Categories.do_compose(α::FinTransformationMap, β::FinTransformation)
  F = dom(α)
  D = codom(F)
  FinTransformationMap(mapvals(α.components, keys=true) do c, f
                         compose(D, f, component(β, c))
                       end, F, codom(β))
end

function Categories.do_composeH(F::FinDomFunctorMap, β::Transformation)
  G, H = dom(β), codom(β)
  FinTransformationMap(mapvals(c -> component(β, c), F.ob_map),
                       compose(F, G, strict=false), compose(F, H, strict=false))
end

function Categories.do_composeH(α::FinTransformationMap, H::Functor)
  F, G = dom(α), codom(α)
  FinTransformationMap(mapvals(f -> hom_map(H, f), α.components),
                       compose(F, H, strict=false), compose(G, H, strict=false))
end

function Base.show(io::IO, α::FinTransformationMap)
  print(io, "FinTransformation(")
  show(io, components(α))
  print(io, ", ")
  Categories.show_domains(io, α, recurse=false)
  print(io, ", ")
  Categories.show_domains(io, dom(α))
  print(io, ")")
end

# Dict utilities
################

# Something like this should be built into Julia, but unfortunately isn't.

function mappairs(kmap, vmap, pairs::T) where T<:AbstractDict
  (dicttype(T))(kmap(k) => vmap(v) for (k,v) in pairs)
end
mappairs(kmap, vmap, vec::AbstractVector) = map(vmap, vec)

function mapvals(f, pairs::T; keys::Bool=false, iter::Bool=false) where T<:AbstractDict
  (iter ? identity : dicttype(T))(if keys
    (k => f(k,v) for (k,v) in pairs)
  else
    (k => f(v) for (k,v) in pairs)
  end)
end
function mapvals(f, collection; keys::Bool=false, iter::Bool=false)
  do_map = iter ? Iterators.map : map
  if keys
    do_map(f, eachindex(collection), values(collection))
  else
    do_map(f, values(collection))
  end
end

dicttype(::Type{T}) where T <: AbstractDict = T.name.wrapper
dicttype(::Type{<:Iterators.Pairs}) = Dict

make_map(f, xs::UnitRange{Int}) = map(f, xs)
make_map(f, xs) = Dict(x => f(x) for x in xs)

end


text:
module XNAColors

macro colornames(s::Symbol, expr::Expr)
    return Expr(:block, esc(:($s = $expr)), [esc(:($(Meta.parse(k)) = $v)) for (k, v) in eval(expr)]...)
end

@colornames(colors, Dict{String, Tuple}(
    "Transparent" => (0, 0, 0, 0) ./ 255,
    "AliceBlue" => (240, 248, 255, 255) ./ 255,
    "AntiqueWhite" => (250, 235, 215, 255) ./ 255,
    "Aqua" => (0, 255, 255, 255) ./ 255,
    "Aquamarine" => (127, 255, 212, 255) ./ 255,
    "Azure" => (240, 255, 255, 255) ./ 255,
    "Beige" => (245, 245, 220, 255) ./ 255,
    "Bisque" => (255, 228, 196, 255) ./ 255,
    "Black" => (0, 0, 0, 255) ./ 255,
    "BlanchedAlmond" => (255, 235, 205, 255) ./ 255,
    "Blue" => (0, 0, 255, 255) ./ 255,
    "BlueViolet" => (138, 43, 226, 255) ./ 255,
    "Brown" => (165, 42, 42, 255) ./ 255,
    "BurlyWood" => (222, 184, 135, 255) ./ 255,
    "CadetBlue" => (95, 158, 160, 255) ./ 255,
    "Chartreuse" => (127, 255, 0, 255) ./ 255,
    "Chocolate" => (210, 105, 30, 255) ./ 255,
    "Coral" => (255, 127, 80, 255) ./ 255,
    "CornflowerBlue" => (100, 149, 237, 255) ./ 255,
    "Cornsilk" => (255, 248, 220, 255) ./ 255,
    "Crimson" => (220, 20, 60, 255) ./ 255,
    "Cyan" => (0, 255, 255, 255) ./ 255,
    "DarkBlue" => (0, 0, 139, 255) ./ 255,
    "DarkCyan" => (0, 139, 139, 255) ./ 255,
    "DarkGoldenrod" => (184, 134, 11, 255) ./ 255,
    "DarkGray" => (169, 169, 169, 255) ./ 255,
    "DarkGreen" => (0, 100, 0, 255) ./ 255,
    "DarkKhaki" => (189, 183, 107, 255) ./ 255,
    "DarkMagenta" => (139, 0, 139, 255) ./ 255,
    "DarkOliveGreen" => (85, 107, 47, 255) ./ 255,
    "DarkOrange" => (255, 140, 0, 255) ./ 255,
    "DarkOrchid" => (153, 50, 204, 255) ./ 255,
    "DarkRed" => (139, 0, 0, 255) ./ 255,
    "DarkSalmon" => (233, 150, 122, 255) ./ 255,
    "DarkSeaGreen" => (143, 188, 139, 255) ./ 255,
    "DarkSlateBlue" => (72, 61, 139, 255) ./ 255,
    "DarkSlateGray" => (47, 79, 79, 255) ./ 255,
    "DarkTurquoise" => (0, 206, 209, 255) ./ 255,
    "DarkViolet" => (148, 0, 211, 255) ./ 255,
    "DeepPink" => (255, 20, 147, 255) ./ 255,
    "DeepSkyBlue" => (0, 191, 255, 255) ./ 255,
    "DimGray" => (105, 105, 105, 255) ./ 255,
    "DodgerBlue" => (30, 144, 255, 255) ./ 255,
    "Firebrick" => (178, 34, 34, 255) ./ 255,
    "FloralWhite" => (255, 250, 240, 255) ./ 255,
    "ForestGreen" => (34, 139, 34, 255) ./ 255,
    "Fuchsia" => (255, 0, 255, 255) ./ 255,
    "Gainsboro" => (220, 220, 220, 255) ./ 255,
    "GhostWhite" => (248, 248, 255, 255) ./ 255,
    "Gold" => (255, 215, 0, 255) ./ 255,
    "Goldenrod" => (218, 165, 32, 255) ./ 255,
    "Gray" => (128, 128, 128, 255) ./ 255,
    "Green" => (0, 128, 0, 255) ./ 255,
    "GreenYellow" => (173, 255, 47, 255) ./ 255,
    "Honeydew" => (240, 255, 240, 255) ./ 255,
    "HotPink" => (255, 105, 180, 255) ./ 255,
    "IndianRed" => (205, 92, 92, 255) ./ 255,
    "Indigo" => (75, 0, 130, 255) ./ 255,
    "Ivory" => (255, 255, 240, 255) ./ 255,
    "Khaki" => (240, 230, 140, 255) ./ 255,
    "Lavender" => (230, 230, 250, 255) ./ 255,
    "LavenderBlush" => (255, 240, 245, 255) ./ 255,
    "LawnGreen" => (124, 252, 0, 255) ./ 255,
    "LemonChiffon" => (255, 250, 205, 255) ./ 255,
    "LightBlue" => (173, 216, 230, 255) ./ 255,
    "LightCoral" => (240, 128, 128, 255) ./ 255,
    "LightCyan" => (224, 255, 255, 255) ./ 255,
    "LightGoldenrodYellow" => (250, 250, 210, 255) ./ 255,
    "LightGray" => (211, 211, 211, 255) ./ 255,
    "LightGreen" => (144, 238, 144, 255) ./ 255,
    "LightPink" => (255, 182, 193, 255) ./ 255,
    "LightSalmon" => (255, 160, 122, 255) ./ 255,
    "LightSeaGreen" => (32, 178, 170, 255) ./ 255,
    "LightSkyBlue" => (135, 206, 250, 255) ./ 255,
    "LightSlateGray" => (119, 136, 153, 255) ./ 255,
    "LightSteelBlue" => (176, 196, 222, 255) ./ 255,
    "LightYellow" => (255, 255, 224, 255) ./ 255,
    "Lime" => (0, 255, 0, 255) ./ 255,
    "LimeGreen" => (50, 205, 50, 255) ./ 255,
    "Linen" => (250, 240, 230, 255) ./ 255,
    "Magenta" => (255, 0, 255, 255) ./ 255,
    "Maroon" => (128, 0, 0, 255) ./ 255,
    "MediumAquamarine" => (102, 205, 170, 255) ./ 255,
    "MediumBlue" => (0, 0, 205, 255) ./ 255,
    "MediumOrchid" => (186, 85, 211, 255) ./ 255,
    "MediumPurple" => (147, 112, 219, 255) ./ 255,
    "MediumSeaGreen" => (60, 179, 113, 255) ./ 255,
    "MediumSlateBlue" => (123, 104, 238, 255) ./ 255,
    "MediumSpringGreen" => (0, 250, 154, 255) ./ 255,
    "MediumTurquoise" => (72, 209, 204, 255) ./ 255,
    "MediumVioletRed" => (199, 21, 133, 255) ./ 255,
    "MidnightBlue" => (25, 25, 112, 255) ./ 255,
    "MintCream" => (245, 255, 250, 255) ./ 255,
    "MistyRose" => (255, 228, 225, 255) ./ 255,
    "Moccasin" => (255, 228, 181, 255) ./ 255,
    "NavajoWhite" => (255, 222, 173, 255) ./ 255,
    "Navy" => (0, 0, 128, 255) ./ 255,
    "OldLace" => (253, 245, 230, 255) ./ 255,
    "Olive" => (128, 128, 0, 255) ./ 255,
    "OliveDrab" => (107, 142, 35, 255) ./ 255,
    "Orange" => (255, 165, 0, 255) ./ 255,
    "OrangeRed" => (255, 69, 0, 255) ./ 255,
    "Orchid" => (218, 112, 214, 255) ./ 255,
    "PaleGoldenrod" => (238, 232, 170, 255) ./ 255,
    "PaleGreen" => (152, 251, 152, 255) ./ 255,
    "PaleTurquoise" => (175, 238, 238, 255) ./ 255,
    "PaleVioletRed" => (219, 112, 147, 255) ./ 255,
    "PapayaWhip" => (255, 239, 213, 255) ./ 255,
    "PeachPuff" => (255, 218, 185, 255) ./ 255,
    "Peru" => (205, 133, 63, 255) ./ 255,
    "Pink" => (255, 192, 203, 255) ./ 255,
    "Plum" => (221, 160, 221, 255) ./ 255,
    "PowderBlue" => (176, 224, 230, 255) ./ 255,
    "Purple" => (128, 0, 128, 255) ./ 255,
    "Red" => (255, 0, 0, 255) ./ 255,
    "RosyBrown" => (188, 143, 143, 255) ./ 255,
    "RoyalBlue" => (65, 105, 225, 255) ./ 255,
    "SaddleBrown" => (139, 69, 19, 255) ./ 255,
    "Salmon" => (250, 128, 114, 255) ./ 255,
    "SandyBrown" => (244, 164, 96, 255) ./ 255,
    "SeaGreen" => (46, 139, 87, 255) ./ 255,
    "SeaShell" => (255, 245, 238, 255) ./ 255,
    "Sienna" => (160, 82, 45, 255) ./ 255,
    "Silver" => (192, 192, 192, 255) ./ 255,
    "SkyBlue" => (135, 206, 235, 255) ./ 255,
    "SlateBlue" => (106, 90, 205, 255) ./ 255,
    "SlateGray" => (112, 128, 144, 255) ./ 255,
    "Snow" => (255, 250, 250, 255) ./ 255,
    "SpringGreen" => (0, 255, 127, 255) ./ 255,
    "SteelBlue" => (70, 130, 180, 255) ./ 255,
    "Tan" => (210, 180, 140, 255) ./ 255,
    "Teal" => (0, 128, 128, 255) ./ 255,
    "Thistle" => (216, 191, 216, 255) ./ 255,
    "Tomato" => (255, 99, 71, 255) ./ 255,
    "Turquoise" => (64, 224, 208, 255) ./ 255,
    "Violet" => (238, 130, 238, 255) ./ 255,
    "Wheat" => (245, 222, 179, 255) ./ 255,
    "White" => (255, 255, 255, 255) ./ 255,
    "WhiteSmoke" => (245, 245, 245, 255) ./ 255,
    "Yellow" => (255, 255, 0, 255) ./ 255,
    "YellowGreen" => (154, 205, 50, 255) ./ 255,
))

end

ext: jl


text:
## MultilevelEstimators.jl : module file
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

module MultilevelEstimators

## dependencies ##

using Dates, DelimitedFiles, Distributed, JLD2, LinearAlgebra, Printf, Random, SpecialFunctions, Statistics

## import statements ##

import Base: filter, ndims, eltype, show, diff, sum, :, print, length, getindex, run, push!, isless, >, <, ≥, ≤, +, zero, one, keys, haskey

import Random: GLOBAL_RNG

import Statistics: mean, var

## export statements ##

export Level, Index

export AbstractIndexSet, SL, ML, FT, TD, HC, ZC, AD, U, get_index_set

export AbstractDistribution, Uniform, Normal, TruncatedNormal, Weibull , transform

export AbstractSampleMethod, MC, QMC

export LatticeRule32, ShiftedLatticeRule, get_point

export Estimator

export History, haskey

export run

## include statements ##

# core

include("check_input.jl")

include("index.jl")

include("index_set.jl")

include("distribution.jl")

include("sample_method.jl")

include("lattice.jl")

include("estimator.jl")

include("options.jl")

include("parse.jl")

include("internals.jl")

include("print.jl")

include("sample.jl")

include("history.jl")

include("run.jl")

include("methods.jl")

end # module


text:
using Test, Batched

@testset "checking BatchedMatrix matmul" begin

    A = BatchedMatrix(rand(3, 4, 2, 3))
    B = BatchedMatrix(rand(4, 3, 2, 3))

    C = A * B

    vA = A.parent
    vB = B.parent
    vC = zeros(3, 3, 2, 3)

    for i in 1:2
        for j in 1:3
            vC[:, :, i, j] = vA[:, :, i, j] * vB[:, :, i, j]
        end
    end

    @test vC ≈ C.parent

end

@testset "checking BatchedMatrix matmul (transposed)" begin

    A = BatchedMatrix(rand(3, 4, 2, 3))
    B = BatchedMatrix(rand(3, 4, 2, 3))

    C = A * transpose(B)

    vA = A.parent
    vB = B.parent
    vC = zeros(3, 3, 2, 3)

    for i in 1:2
        for j in 1:3
            vC[:, :, i, j] = vA[:, :, i, j] * transpose(vB[:, :, i, j])
        end
    end

    @test vC ≈ C.parent

end

@testset "checking BatchedMatrix matmul (scale)" begin

    A = BatchedMatrix(rand(3, 4, 2, 3))
    B = BatchedUniformScaling(rand(2, 3))

    C = B * A

    vA = A.parent
    vB = B.scalars
    vC = zeros(3, 4, 2, 3)

    for i in 1:2
        for j in 1:3
            vC[:, :, i, j] = vB[i, j] * vA[:, :, i, j]
        end
    end

    @test vC ≈ C.parent

end

@testset "checking BatchedMatrix matmul (scale & transpose)" begin

    A = BatchedMatrix(rand(3, 4, 2, 3))
    B = BatchedUniformScaling(rand(2, 3))

    C = B * transpose(A)

    vA = A.parent
    vB = B.scalars
    vC = zeros(4, 3, 2, 3)

    for i in 1:2
        for j in 1:3
            vC[:, :, i, j] = vB[i, j] * transpose(vA[:, :, i, j])
        end
    end

    @test vC ≈ C.parent
end


text:
function simpson_integrate(u, stepsize)
    numpts = length(u)
    @assert iseven(numpts - 1)
    mid = div(numpts,2)
    integral = 0.0
    for i = 1:mid
        integral += stepsize / 3.0 * (u[2i-1] + 4.0 * u[2i] + u[2i+1])
    end
    return integral
end

max_stars_repo_name: ArjunNarayanan/ExploreJacobi

max_issues_repo_name: ArjunNarayanan/ExploreJacobi

max_forks_repo_name: ArjunNarayanan/ExploreJacobi


text:
push!(LOAD_PATH,"../src/")
using Documenter, HarmonicPowerModels

# A flag to check if we are running in a GitHub action.
const _IS_GITHUB_ACTIONS = get(ENV, "GITHUB_ACTIONS", "false") == "true"

# Pass --pdf to build the PDF. On GitHub actions, we always build the PDF.
const _PDF = findfirst(isequal("--pdf"), ARGS) !== nothing || _IS_GITHUB_ACTIONS

const _PAGES = [
    "Home" => "index.md",
];

@time Documenter.makedocs(
    modules = [HarmonicPowerModels],
    format = Documenter.HTML(mathengine = Documenter.MathJax()),
    sitename = "HarmonicPowerModels",
    authors = "Tom Van Acker, Frederik Geth and contributors.",
    pages = _PAGES,
    )

if _PDF
    # latex_platform = _IS_GITHUB_ACTIONS ? "docker" : "native"
    latex_platform = "docker"
    @time Documenter.makedocs(
        sitename = "HarmonicPowerModels",
        authors = "The HarmonicPowerModels core developers and contributors",
        format = Documenter.LaTeX(platform = latex_platform),
        build = "latex_build",
    )
    # Hack for deploying: copy the pdf (and only the PDF) into the HTML build
    # directory! We don't want to copy everything in `latex_build` because it
    # includes lots of extraneous LaTeX files.
    cp(
        joinpath(@__DIR__, "latex_build", "HarmonicPowerModels.pdf"),
        joinpath(@__DIR__, "build", "HarmonicPowerModels.pdf"),
    )
end

Documenter.deploydocs(
    repo = "github.com/timmyfaraday/HarmonicPowerModels.jl.git",
    push_preview = true,
)


text:
"""
    NexPicture(nexSerial, name; pid=pid, cid=cid)

A Nextion Picture UI control.

On Nextion side, it have the following properties:
- pic
"""
struct NexPicture <: AbstractNexObject
    _nid::NexID

    viewable::IViewable
    picturable::IPicturable

    function NexPicture(nexSerial::T, name::Name; pid=PageID(), cid=ComponentID()) where {T <: AbstractNexSerial}
        nid = NexID(nexSerial, name, pid, cid)
        new(nid, IViewable(nid), IPicturable(nid))
    end
end

function setproperty!(obj::NexPicture, property::Symbol, new_val)
    _setcommonproperty!(obj, property, new_val)
end

function getproperty(obj::NexPicture, property::Symbol)
    _getcommonproperty(obj, property)
end


text:
using JLD2
using Printf
using GLMakie

Nx = 128
Ny = 60

bathymetry_path = "bathy_128x60var4.bin"
east_west_stress_path = "off_TAUXvar1.bin"
north_south_stress_path = "off_TAUY.bin"
sea_surface_temperature_path="sst25_128x60x12.bin"

Nbytes = sizeof(Float32) * Nx * Ny
bathymetry = reshape(bswap.(reinterpret(Float32, read(bathymetry_path, Nbytes))), (Nx, Ny))

τˣ = reshape(bswap.(reinterpret(Float32, read(east_west_stress_path,  12Nbytes))), (Nx, Ny, 12))
τʸ = reshape(bswap.(reinterpret(Float32, read(north_south_stress_path, 12Nbytes))), (Nx, Ny, 12))
target_sea_surface_temperature = reshape(bswap.(reinterpret(Float32, read(sea_surface_temperature_path, 12Nbytes))), (Nx, Ny, 12))

τˣ = τˣ[:, :, 1]
τʸ = τʸ[:, :, 1]
target_sea_surface_temperature = target_sea_surface_temperature[:, :, 1]

# bathymetry = Array{Float64, 2}(bathymetry)
# τˣ = Array{Float64, 2}(τˣ)
# τʸ = Array{Float64, 2}(τʸ)
# target_sea_surface_temperature = Array{Float64, 2}(target_sea_surface_temperature)

bathymetry_file = jldopen("earth_bathymetry_128_60.jld2", "a+")
bathymetry_file["bathymetry"] = bathymetry
close(bathymetry_file)

println("Bathymetry: min= ", minimum(bathymetry)," , max= ", maximum(bathymetry))
println("τˣ: min= ", minimum(τˣ)," , max= ", maximum(τˣ))
println("τʸ: min= ", minimum(τʸ)," , max= ", maximum(τʸ))
println("SST★: min= ", minimum(target_sea_surface_temperature)," , max= ",maximum(target_sea_surface_temperature))

fig = Figure(resolution = (1920, 1080))
ax = fig[1, 1] = LScene(fig, title = "Bathymetry")
heatmap!(ax, bathymetry)

ax = fig[2, 1] = LScene(fig, title = "Bathymetry")
heatmap!(ax, target_sea_surface_temperature)

ax = fig[1, 2] = LScene(fig, title = "East-west wind stress")
heatmap!(ax, τˣ)

ax = fig[2, 2] = LScene(fig, title = "North-south wind stress")
heatmap!(ax, τʸ)

display(fig)


text:
"""
    proteinCross(genotype1::GEGenotype, genotype2::GEGenotype,
                            gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                            integerCross::Function = singlePointCross, varargs...)::Array{GEGenotype}

Perform codons based cross for two **Grammatical Evolution** Individuals, by means of
an integer cross.

# Self-provided Arguments
- `genotype1::GEGenotype`: first parent genotype.
- `genotype2::GEGenotype`: second parent genotype.
- `gpExperimentInfo::GEInfo`: information about grammatical evolution.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

# User Arguments
- `integerCross::Function`: integer crossover function.
- `varargs`: arguments for integerCross.

!!! note
    `Selfprovided Arguments` are already provided by the library, so only `User Arguments` must be provided.
"""
function proteinCross(genotype1::GEGenotype, genotype2::GEGenotype,
                            gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                            integerCross::Function = singlePointCross, varargs...)::Array{GEGenotype}

    protein1, introns1, protein2, introns2 = pruneProtein(genotype1, genotype2)

    children = integerCross(protein1, protein2, rng, varargs...)

    child1 = generateTree(gpExperimentInfo, rng::Random.AbstractRNG, vcat(children[1]._representation, introns1))
    child2 = generateTree(gpExperimentInfo, rng::Random.AbstractRNG, vcat(children[2]._representation, introns2))
    return [child1, child2]
end

"""
    proteinCrossGrow(genotype1::GEGenotype, genotype2::GEGenotype,
                            gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                            integerCross::Function = singlePointCross, varargs...)::Array{GEGenotype}

As proteinCross, but instead, a growGenerator is used to ensure that all individuals are valids.
It is advisable to used along proteinMutationGrow.

# Self-provided Arguments
- `genotype1::GEGenotype`: first parent genotype.
- `genotype2::GEGenotype`: second parent genotype.
- `gpExperimentInfo::GEInfo`: information about grammatical evolution.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

# User Arguments
- `integerCross::Function`: integer crossover function.
- `varargs`: arguments for integerCross.

!!! note
    `Selfprovided Arguments` are already provided by the library, so only `User Arguments` must be provided.
"""
function proteinCrossGrow(genotype1::GEGenotype, genotype2::GEGenotype,
                            gpExperimentInfo::GEInfo, rng::Random.AbstractRNG,
                            integerCross::Function = singlePointCross, varargs...)

    protein1, introns1, protein2, introns2 = pruneProtein(genotype1, genotype2)

    children = integerCross(protein1, protein2, rng, varargs...)

    child1 = growGenerator(gpExperimentInfo, rng, vcat(children[1]._representation, introns1))
    child2 = growGenerator(gpExperimentInfo, rng, vcat(children[2]._representation, introns2))

    return [child1, child2]
end

# TODO Fails
"""
function onePointCrossDepth(parent1::GEGenotype, parent2::GEGenotype,
                       gpExperimentInfo::GEInfo, rng::Random.AbstractRNG)

    maxDepth = gpExperimentInfo._maxDepth
    child1 = copyGenotype(parent1)
    child2 = copyGenotype(parent2)

    preOrder1 = getPreOrder(child1._tree)
    preOrder2 = getPreOrder(child2._tree)

    posibleValues = collect(1:length(preOrder1))
    posibleValues = filter(x->!isLeaf(preOrder1[x]), posibleValues)

    heights = zeros(UInt16, length(preOrder2))

    heights1 = getHeights(child1._tree)
    heights2 = getHeights(child2._tree)

    found = false

    while !isempty(posibleValues) && !found
        selected = rand(rng, posibleValues)
        posibleValues[selected], posibleValues[end] = posibleValues[end], posibleValues[selected]

        nodeIndex = pop!(posibleValues)
        symbol1 = getSymbol(preOrder1[nodeIndex])
        depth1 = getDepth(preOrder1[nodeIndex])
        height1 = heights1[nodeIndex]

        posiblesNodes = Array{GERep}(undef,0)

        for i=1:length(preOrder2)
            node = preOrder2[i]
            if getSymbol(node) == symbol1
                if height1 == heights2[i]
                    found = true
                    changeNodes!(preOrder1[nodeIndex], node)
                elseif height1 > heights2[i]
                    if height1 + getDepth(node) <= maxDepth
                        found = true
                        changeNodes!(preOrder1[nodeIndex], node)
                    end
                else
                    if height2 + depth1 <= maxDepth
                        found = true
                        changeNodes!(preOrder1[nodeIndex], node)
                    end
                end
            end
        end
    end

    return [child1, child2]
end
"""

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
cd(@__DIR__)
open("cwrapper.cpp", "w") do io
    println(io, """
    #include "tetgen.h" // Defined tetgenio, tetrahedralize().
    """)

    for (REAL, ext) in ["double" => "f64"]
        println(io, """
        typedef struct {
          int firstnumber; // 0 or 1, default 0.
          int mesh_dim;    // must be 3.

          $REAL *pointlist;
          $REAL *pointattributelist;
          $REAL *pointmtrlist;
          int *pointmarkerlist;
          int numberofpoints;
          int numberofpointattributes;
          int numberofpointmtrs;

          int *tetrahedronlist;
          $REAL *tetrahedronattributelist;
          $REAL *tetrahedronvolumelist;
          int *neighborlist;
          int numberoftetrahedra;
          int numberofcorners;
          int numberoftetrahedronattributes;

          tetgenio::facet *facetlist;
          int *facetmarkerlist;
          int numberoffacets;

          $REAL *holelist;
          int numberofholes;

          $REAL *regionlist;
          int numberofregions;

          $REAL *facetconstraintlist;
          int numberoffacetconstraints;

          $REAL *segmentconstraintlist;
          int numberofsegmentconstraints;

          int *trifacelist;
          int *trifacemarkerlist;
          int numberoftrifaces;

          int *edgelist;
          int *edgemarkerlist;
          int numberofedges;
        } TetGenIO$ext;

        void copy_tetio(tetgenio* in, TetGenIO$ext* out){
            out->firstnumber = in->firstnumber;
            out->mesh_dim = in->mesh_dim;

            out->pointlist = in->pointlist;
            out->pointattributelist = in->pointattributelist;
            out->pointmtrlist = in->pointmtrlist;
            out->pointmarkerlist = in->pointmarkerlist;
            out->numberofpoints = in->numberofpoints;
            out->numberofpointattributes = in->numberofpointattributes;
            out->numberofpointmtrs = in->numberofpointmtrs;

            out->tetrahedronlist = in->tetrahedronlist;
            out->tetrahedronattributelist = in->tetrahedronattributelist;
            out->tetrahedronvolumelist = in->tetrahedronvolumelist;
            out->neighborlist = in->neighborlist;
            out->numberoftetrahedra = in->numberoftetrahedra;
            out->numberofcorners = in->numberofcorners;
            out->numberoftetrahedronattributes = in->numberoftetrahedronattributes;

            out->facetlist = in->facetlist;
            out->facetmarkerlist = in->facetmarkerlist;
            out->numberoffacets = in->numberoffacets;

            out->holelist = in->holelist;
            out->numberofholes = in->numberofholes;

            out->regionlist = in->regionlist;
            out->numberofregions = in->numberofregions;

            out->facetconstraintlist = in->facetconstraintlist;
            out->numberoffacetconstraints = in->numberoffacetconstraints;

            out->segmentconstraintlist = in->segmentconstraintlist;
            out->numberofsegmentconstraints = in->numberofsegmentconstraints;

            out->trifacelist = in->trifacelist;
            out->trifacemarkerlist = in->trifacemarkerlist;
            out->numberoftrifaces = in->numberoftrifaces;

            out->edgelist = in->edgelist;
            out->edgemarkerlist = in->edgemarkerlist;
            out->numberofedges = in->numberofedges;

        }
        void copy_tetio(TetGenIO$ext* in, tetgenio* out){
            out->firstnumber = in->firstnumber;
            out->mesh_dim = in->mesh_dim;

            out->pointlist = in->pointlist;
            out->pointattributelist = in->pointattributelist;
            out->pointmtrlist = in->pointmtrlist;
            out->pointmarkerlist = in->pointmarkerlist;
            out->numberofpoints = in->numberofpoints;
            out->numberofpointattributes = in->numberofpointattributes;
            out->numberofpointmtrs = in->numberofpointmtrs;

            out->tetrahedronlist = in->tetrahedronlist;
            out->tetrahedronattributelist = in->tetrahedronattributelist;
            out->tetrahedronvolumelist = in->tetrahedronvolumelist;
            out->neighborlist = in->neighborlist;
            out->numberoftetrahedra = in->numberoftetrahedra;
            out->numberofcorners = in->numberofcorners;
            out->numberoftetrahedronattributes = in->numberoftetrahedronattributes;

            out->facetlist = in->facetlist;
            out->facetmarkerlist = in->facetmarkerlist;
            out->numberoffacets = in->numberoffacets;

            out->holelist = in->holelist;
            out->numberofholes = in->numberofholes;

            out->regionlist = in->regionlist;
            out->numberofregions = in->numberofregions;

            out->facetconstraintlist = in->facetconstraintlist;
            out->numberoffacetconstraints = in->numberoffacetconstraints;

            out->segmentconstraintlist = in->segmentconstraintlist;
            out->numberofsegmentconstraints = in->numberofsegmentconstraints;

            out->trifacelist = in->trifacelist;
            out->trifacemarkerlist = in->trifacemarkerlist;
            out->numberoftrifaces = in->numberoftrifaces;

            out->edgelist = in->edgelist;
            out->edgemarkerlist = in->edgemarkerlist;
            out->numberofedges = in->numberofedges;

        }
        extern "C" {

          TetGenIO$ext tetrahedralize$ext(TetGenIO$ext jl_in, char* command){
            tetgenio in, out;
            copy_tetio(&jl_in, &in);
            tetrahedralize(command, &in, &out);
            TetGenIO$ext jl_out;
            copy_tetio(&out, &jl_out);
            in.initialize(); // don't free pointers from julia!
            out.initialize(); // don't free pointers for julia!
            return jl_out;
          }

        }

        """)
    end
end


text:
abstract type AbstractSDEKernel
end

struct SDEKernel{fType,gType,tType,pType,NRPType} <: AbstractSDEKernel
  f::fType
  g::gType
  trange::tType
  p::pType
  noise_rate_prototype::NRPType
  constant_diffusity::Bool
end

struct SDEKernel!{fType,gType,gstep!Type,tType,pType,NRPType,wsType} <: AbstractSDEKernel
  f::fType
  g::gType
  gstep!::gstep!Type
  trange::tType
  p::pType
  noise_rate_prototype::NRPType
  ws::wsType
  constant_diffusity::Bool
end
function make_gstep(g) #g has signature  g(du,u,p,t)
  function gstep!(dx, ws, x, p, t, dw, noise_rate_prototype)
    if noise_rate_prototype===nothing
      g(ws, x, p, t)
      dx .+= ws .* dw
    else
      g(ws, x, p,t)
      dx .+= ws * dw
    end
  end
end

function SDEKernel(f,g,trange,p=nothing,noise_rate_prototype=nothing,constant_diffusity=false)
  SDEKernel{typeof(f),typeof(g),typeof(trange),typeof(p),typeof(noise_rate_prototype)}(f,g,trange,p,noise_rate_prototype,constant_diffusity)
end

function SDEKernel!(f,g,gstep!,trange,p=nothing,noise_rate_prototype=nothing,constant_diffusity=false; ws = copy(noise_rate_prototype))
  SDEKernel!{typeof(f),typeof(g),typeof(gstep!),typeof(trange),typeof(p),typeof(noise_rate_prototype),typeof(ws)}(f,g,gstep!, trange,p,noise_rate_prototype, ws, constant_diffusity)
end
abstract type AbstractFilteringAlgorithm end

struct CovarianceFilter <: AbstractFilteringAlgorithm end
struct InformationFilter <: AbstractFilteringAlgorithm end
struct LyapunovFilter <: AbstractFilteringAlgorithm end

struct solLyapunov{kernelType,ΣType,νTType,PTType,CTType}
  ktilde::kernelType
  Σ::ΣType
  νT::νTType
  PT::PTType
  cT::CTType
end

struct Message{kernelType,solType,sol2Type,tType,filterType}
  ktilde::kernelType
  sol::solType
  soldis::sol2Type
  ts::tType
  filter::filterType
end

function Message(sol, sdekernel::SDEKernel, filter, apply_timechange=false)
  soldis = construct_discrete_sol(sol)
  if OrdinaryDiffEq.isadaptive(sol.alg) || apply_timechange
    ts = reverse(sol.t)
  else
    ts = sdekernel.trange
  end
  Message{typeof(sdekernel),typeof(sol),typeof(soldis),typeof(ts),typeof(filter)}(sdekernel,
    sol,soldis,ts,filter)
end

struct GuidingDriftCache{kernelType,messageType}
  k::kernelType
  message::messageType
end

struct GuidingDiffusionVectorCache{gType}
  g::gType
end

struct GuidingDiffusionCache{gType,sizeType}
  g::gType
  padded_size::sizeType
end

abstract type AbstractRegression{inplace} end

struct Regression{kernelType,pJType,ifuncType,pfType,θType,duType,inplace} <: AbstractRegression{inplace}
  k::kernelType
  fjac::pJType
  ϕ0func::ifuncType
  pf::pfType
  θ::θType
  dy::duType
  isscalar::Bool
end

function Regression{iip}(sdekernel::Union{SDEKernel,SDEProblem}; paramjac=nothing,
    intercept=nothing, θ=sdekernel.p, yprototype=nothing,
    dyprototype=nothing, m=nothing) where iip

  if paramjac === nothing
    if iip
      pf = ParamJacobianWrapper2(sdekernel.f,first(sdekernel.trange),yprototype)
    else
      pf = ParamJacobianWrapper(sdekernel.f,first(sdekernel.trange),yprototype)
    end
  else
    pf = nothing
  end

  if !iip
    dy = nothing
  else
    dyprototype === nothing && error("dyprototype needs to be known for using inplace functions.")
    if m===nothing
      # default, diagonal noise
      dy = similar(dyprototype)
      isscalar = false
    else
      if m==1
        # scalar noise case
        dy = similar(dyprototype, length(dyprototype))
        isscalar = true
      else
        # non-diagonal noise
        dy = similar(dyprototype, length(dyprototype,m))
        isscalar = false
      end
    end
  end

  if m===nothing
    isscalar = false
  else
    isscalar = (m==1)
  end

  Regression{typeof(sdekernel),typeof(paramjac),typeof(intercept),typeof(pf),
    typeof(θ),typeof(dy),iip}(sdekernel,
    paramjac,intercept,pf,θ,dy,isscalar)
end

function Regression(sdekernel;kwargs...)
  iip=DiffEqBase.isinplace(sdekernel.g,4)
  Regression{iip}(sdekernel;kwargs...)
end

mutable struct Regression!{kernelType,pJType,ifuncType,phiType,uType,pfType,θType,duType,inplace} <: AbstractRegression{inplace}
  k::kernelType
  fjac!::pJType
  ϕ0func!::ifuncType
  ϕ::phiType
  ϕ0::uType
  y::uType
  y2::uType
  pf::pfType
  θ::θType
  dy::duType
  isscalar::Bool
end

function Regression!{iip}(sdekernel::Union{SDEKernel,SDEProblem}, yprototype;
      paramjac_prototype=nothing, paramjac=nothing, intercept=nothing, isdiagonal=true,
      θ=sdekernel.p,
      dyprototype=yprototype, m=length(yprototype)) where iip

  y = similar(yprototype)
  y2 = similar(y)
  ϕ0 = similar(y)

  if paramjac_prototype !== nothing
    ϕ = similar(paramjac_prototype)
  else
    ϕ = zeros((length(yprototype),length(θ)))
  end

  if paramjac === nothing
    if iip
      pf = ParamJacobianWrapper2(sdekernel.f,first(sdekernel.trange),y)
    else
      pf = ParamJacobianWrapper(sdekernel.f,first(sdekernel.trange),y)
    end
  else
    pf = nothing
  end

  if !iip
    dy = nothing
  else
    dyprototype === nothing && error("dyprototype needs to be known for using inplace functions.")
    if m==length(yprototype)
      # default, diagonal noise
      dy = similar(dyprototype)
    else
      if m==1
        # scalar noise case
        dy = similar(dyprototype, length(dyprototype))
      else
        # non-diagonal noise
        dy = similar(dyprototype, length(dyprototype,m))
      end
    end
  end

  isscalar = (m==1 && length(yprototype)!=1)

  Regression!{typeof(sdekernel),typeof(paramjac),typeof(intercept),typeof(ϕ),
    typeof(y),typeof(pf),typeof(θ),typeof(dy),iip}(sdekernel,paramjac,intercept,
    ϕ,ϕ0,y,y2,pf,θ,dy,isscalar)
end

function Regression!(sdekernel,yprototype;kwargs...)
  iip=DiffEqBase.isinplace(sdekernel.g,4)
  Regression!{iip}(sdekernel,yprototype;kwargs...)
end

# internal solvers
abstract type AbstractInternalSolver end

struct EulerMaruyama! <: AbstractInternalSolver end

struct DefaultForwardGuidingP end

"""
    tangent!(du, u, dz, P)

!! May change `du`, but has to return it. !!

"""
function tangent!
end

struct GuidedSDE{kType, mType}
  kernel::kType
  message::mType
end

struct GuidedSDE!{kType, mType}
  kernel::kType
  message::mType
end


text:
"""
    WeibullExpert(k, θ)

PDF:

```math
f(x; k, \\theta) = \\frac{k}{\\theta} \\left( \\frac{x}{\\theta} \\right)^{k-1} e^{-(x/\\theta)^k},
\\quad x \\geq 0
```

See also: [Weibull Distribution](https://en.wikipedia.org/wiki/Weibull_distribution) (Wikipedia)

"""
struct WeibullExpert{T<:Real} <: NonZIContinuousExpert
    k::T
    θ::T
    WeibullExpert{T}(k::T, θ::T) where {T<:Real} = new{T}(k, θ)
end

function WeibullExpert(k::T, θ::T; check_args=true) where {T <: Real}
    check_args && @check_args(WeibullExpert, k >= one(k) && θ > zero(θ))
    return WeibullExpert{T}(k, θ)
end

## Outer constructors
WeibullExpert(k::Real, θ::Real) = WeibullExpert(promote(k, θ)...)
WeibullExpert(k::Integer, θ::Integer) = WeibullExpert(float(k), float(θ))
WeibullExpert() = WeibullExpert(2.0, 1.0)

## Conversion
function convert(::Type{WeibullExpert{T}}, k::S, θ::S) where {T <: Real, S <: Real}
    WeibullExpert(T(k), T(θ))
end
function convert(::Type{WeibullExpert{T}}, d::WeibullExpert{S}) where {T <: Real, S <: Real}
    WeibullExpert(T(d.k), T(d.θ), check_args=false)
end
copy(d::WeibullExpert) = WeibullExpert(d.k, d.θ, check_args=false)

## Loglikelihood of Expert
logpdf(d::WeibullExpert, x...) = Distributions.logpdf.(Distributions.Weibull(d.k, d.θ), x...)
pdf(d::WeibullExpert, x...) = Distributions.pdf.(Distributions.Weibull(d.k, d.θ), x...)
logcdf(d::WeibullExpert, x...) = Distributions.logcdf.(Distributions.Weibull(d.k, d.θ), x...)
cdf(d::WeibullExpert, x...) = Distributions.cdf.(Distributions.Weibull(d.k, d.θ), x...)

## expert_ll, etc
expert_ll_exact(d::WeibullExpert, x::Real) = LRMoE.logpdf(d, x)
function expert_ll(d::WeibullExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_ll = (yl == yu) ? logpdf.(d, yl) : logcdf.(d, yu) + log1mexp.(logcdf.(d, yl) - logcdf.(d, yu))
    expert_ll = (tu == 0.) ? -Inf : expert_ll
    return expert_ll
end
function expert_tn(d::WeibullExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn = (tl == tu) ? logpdf.(d, tl) : logcdf.(d, tu) + log1mexp.(logcdf.(d, tl) - logcdf.(d, tu))
    expert_tn = (tu == 0.) ? -Inf : expert_tn
    return expert_tn
end
function expert_tn_bar(d::WeibullExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn_bar = (tl == tu) ? 0.0 : log1mexp.(logcdf.(d, tu) + log1mexp.(logcdf.(d, tl) - logcdf.(d, tu)))
    return expert_tn_bar
end

exposurize_expert(d::WeibullExpert; exposure = 1) = d

## Parameters
params(d::WeibullExpert) = (d.k, d.θ)
function params_init(y, d::WeibullExpert)
    pos_idx = (y .> 0.0)

    function init_obj(logk, y)
        n = length(y)
        k_tmp = exp(logk)
        θ_tmp = ( sum(y.^k_tmp) / n)^(1/k_tmp)
        return -1 * (n*logk - n*log(θ_tmp) + k_tmp*sum(log.(y)) - n*(k_tmp-1)*log(θ_tmp) - 1/(θ_tmp^k_tmp)*sum(y.^k_tmp) )
    end

    logk_init = Optim.minimizer( Optim.optimize(x -> init_obj(x, y[pos_idx]),  0.0, 2.0 ))

    k_init = exp(logk_init)
    θ_init = ( sum(y[pos_idx].^k_init) / length(y[pos_idx]))^(1/k_init)

    try
        return WeibullExpert(k_init, θ_init)
    catch;
        WeibullExpert()
    end
end

## KS stats for parameter initialization
function ks_distance(y, d::WeibullExpert)
    p_zero = sum(y .== 0.0) / sum(y .>= 0.0)
    return max(abs(p_zero-0.0), (1-0.0)*HypothesisTests.ksstats(y[y .> 0.0], Distributions.Weibull(d.k, d.θ))[2])
end

## Simululation
sim_expert(d::WeibullExpert) = Distributions.rand(Distributions.Weibull(d.k, d.θ), 1)[1]

## penalty
penalty_init(d::WeibullExpert) = [2.0 10.0 1.0 Inf]
no_penalty_init(d::WeibullExpert) = [1.0 Inf 1.0 Inf]
penalize(d::WeibullExpert, p) = (p[1]-1)*log(d.k) - d.k/p[2] + (p[3]-1)*log(d.θ) - d.θ/p[4]

## statistics
mean(d::WeibullExpert) = mean(Distributions.Weibull(d.k, d.θ))
var(d::WeibullExpert) = var(Distributions.Weibull(d.k, d.θ))
quantile(d::WeibullExpert, p) = quantile(Distributions.Weibull(d.k, d.θ), p)
lev(d::WeibullExpert, u) = d.θ*gamma_inc(float(1/d.k+1), (u/d.θ)^d.k, 0)[1]*gamma(float(1/d.k+1)) + u*exp(-(u/d.θ)^d.k)
excess(d::WeibullExpert, u) = mean(d) - lev(d, u)

## Misc functions for E-Step

function _int_logy_func(d::WeibullExpert, x)
    return (iszero(x) || isinf(x)) ? 0.0 : log(x) * pdf.(d, x)
end

function _int_obs_logY_raw(d::WeibullExpert, yl, yu)
    if yl == yu
        return log(yl) * pdf.(d, yl)
    else
        return quadgk.(x -> _int_logy_func(d, x), yl, yu, rtol = 1e-8)[1]
    end
end

function _int_lat_logY_raw(d::WeibullExpert, tl, tu)
    return (tl==0 ? 0.0 : quadgk.(x -> _int_logy_func(d, x), 0.0, tl, rtol = 1e-8)[1]) + (isinf(tu) ? 0.0 : quadgk.(x -> _int_logy_func(d, x), tu, Inf, rtol = 1e-8)[1])
end

function _int_powy_func(d_old::WeibullExpert, k_new, l, u)
    term1 = isinf(l) ? 0.0 : gamma((d_old.k+k_new)/d_old.k) * gamma_inc((d_old.k+k_new)/d_old.k, (l/d_old.θ)^(d_old.k), 1)[2]
    term2 = isinf(u) ? 0.0 : gamma((d_old.k+k_new)/d_old.k) * gamma_inc((d_old.k+k_new)/d_old.k, (u/d_old.θ)^(d_old.k), 1)[2]
    return (d_old.θ)^(k_new) * (term1 - term2)
end

function _int_obs_powY(d_old::WeibullExpert, k_new, yl, yu, expert_ll_pos)
    if yl == yu
        return yl^k_new
    else
        return exp(-expert_ll_pos) * _int_powy_func(d_old, k_new, yl, yu)
    end
end

function _int_lat_powY(d_old::WeibullExpert, k_new, tl, tu, expert_tn_bar_pos)
    return exp(-expert_tn_bar_pos) * (_int_powy_func(d_old, k_new, 0.0, tl) + _int_powy_func(d_old, k_new, tu, Inf))
end

function _weibull_k_to_λ(k, sum_term_zkz, sum_term_zkzpowy;
                        penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf])

    return (sum_term_zkzpowy / sum_term_zkz)^(1/k)
end

function _weibull_optim_k(logk,
                        d_old,
                        tl, yl, yu, tu,
                        expert_ll_pos, expert_tn_bar_pos,
                        z_e_obs, z_e_lat, k_e,
                        logY_e_obs, logY_e_lat;
                        penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf])
    # Optimization in log scale for unconstrained computation
    k_tmp = exp(logk)

    # Further E-step
    powY_e_obs = _int_obs_powY.(d_old, k_tmp, yl, yu, expert_ll_pos)
    powY_e_lat = _int_lat_powY.(d_old, k_tmp, tl, tu, expert_tn_bar_pos)
    nan2num(powY_e_obs, 0.0) # get rid of NaN
    nan2num(powY_e_lat, 0.0) # get rid of NaN

    term_zkz = z_e_obs .+ (z_e_lat .* k_e)
    term_zkz_logY = (z_e_obs .* logY_e_obs) .+ (z_e_lat .* k_e .* logY_e_lat)
    term_zkz_powY = (z_e_obs .* powY_e_obs) .+ (z_e_lat .* k_e .* powY_e_lat)

    sum_term_zkz = sum(term_zkz)[1]
    sum_term_zkz_logY = sum(term_zkz_logY)[1]
    sum_term_zkz_powY = sum(term_zkz_powY)[1]

    θ_tmp = _weibull_k_to_λ(k_tmp, sum_term_zkz, sum_term_zkz_powY, penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    obj = log(k_tmp)*sum_term_zkz - k_tmp*log(θ_tmp)*sum_term_zkz + (k_tmp-1)*sum_term_zkz_logY - (θ_tmp^(-k_tmp))*sum_term_zkz_powY
    p = penalty ? (pen_pararms_jk[1]-1)*log(k_tmp) - k_tmp/pen_pararms_jk[2] + (pen_pararms_jk[3]-1)*log(θ_tmp) - θ_tmp/pen_pararms_jk[4] : 0.0
    return (obj + p) * (-1.0)

end

## EM: M-Step
function EM_M_expert(d::WeibullExpert,
                     tl, yl, yu, tu,
                     exposure,
                     z_e_obs, z_e_lat, k_e;
                     penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf])

    expert_ll_pos = expert_ll.(d, tl, yl, yu, tu)
    expert_tn_bar_pos = expert_tn_bar.(d, tl, yl, yu, tu)

    # Further E-Step
    yl_yu_unique = unique_bounds(yl, yu)
    int_obs_logY_tmp = _int_obs_logY_raw.(d, yl_yu_unique[:,1], yl_yu_unique[:,2])
    logY_e_obs = exp.(-expert_ll_pos) .* int_obs_logY_tmp[match_unique_bounds(hcat(vec(yl), vec(yu)), yl_yu_unique)]
    nan2num(logY_e_obs, 0.0) # get rid of NaN

    tl_tu_unique = unique_bounds(tl, tu)
    int_lat_logY_tmp = _int_lat_logY_raw.(d, tl_tu_unique[:,1], tl_tu_unique[:,2])
    logY_e_lat = exp.(-expert_tn_bar_pos) .* int_lat_logY_tmp[match_unique_bounds(hcat(vec(tl), vec(tu)), tl_tu_unique)]
    nan2num(logY_e_lat, 0.0) # get rid of NaN

    # Update parameters
    pos_idx = (yu .!= 0.0)
    logk_new = Optim.minimizer( Optim.optimize(x -> _weibull_optim_k(x, d,
                                                            tl[pos_idx], yl[pos_idx], yu[pos_idx], tu[pos_idx],
                                                            expert_ll_pos[pos_idx],  expert_tn_bar_pos[pos_idx],
                                                            z_e_obs[pos_idx], z_e_lat[pos_idx], k_e[pos_idx],
                                                            logY_e_obs[pos_idx], logY_e_lat[pos_idx],
                                                            penalty = penalty, pen_pararms_jk = pen_pararms_jk),
                                               max(log(d.k)-2.0, 0.0), log(d.k)+2.0 )) # ,
                                               # log(d.k)-2.0, log(d.k)+2.0 )) # ,
                                               # GoldenSection() ))
                                               # , rel_tol = 1e-8) )
    k_new = exp(logk_new)

    # Further E-step
    powY_e_obs = _int_obs_powY.(d, k_new, yl, yu, expert_ll_pos)
    powY_e_lat = _int_lat_powY.(d, k_new, tl, tu, expert_tn_bar_pos)
    nan2num(powY_e_obs, 0.0) # get rid of NaN
    nan2num(powY_e_lat, 0.0) # get rid of NaN
    inf2num(powY_e_obs, 0.0) # get rid of Inf
    inf2num(powY_e_lat, 0.0) # get rid of Inf

    term_zkz = z_e_obs .+ (z_e_lat .* k_e)
    # term_zkz_logY = (z_e_obs .* logY_e_obs) .+ (z_e_lat .* k_e .* logY_e_lat)
    term_zkz_powY = (z_e_obs .* powY_e_obs) .+ (z_e_lat .* k_e .* powY_e_lat)

    sum_term_zkz = sum(term_zkz)[1]
    # sum_term_zkz_logY = sum(term_zkz_logY)[1]
    sum_term_zkz_powY = sum(term_zkz_powY)[1]

    θ_new = _weibull_k_to_λ(k_new, sum_term_zkz, sum(sum_term_zkz_powY)[1], penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    return WeibullExpert(k_new, θ_new)

end

## EM: M-Step, exact observations
function _weibull_optim_k_exact(logk,
                        ye,
                        z_e_obs,
                        logY_e_obs;
                        penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf])
    # Optimization in log scale for unconstrained computation
    k_tmp = exp(logk)

    # Further E-step
    powY_e_obs = ye.^(k_tmp)
    nan2num(powY_e_obs, 0.0)

    term_zkz = z_e_obs
    term_zkz_logY = z_e_obs .* logY_e_obs
    term_zkz_powY = z_e_obs .* powY_e_obs

    sum_term_zkz = sum(term_zkz)[1]
    sum_term_zkz_logY = sum(term_zkz_logY)[1]
    sum_term_zkz_powY = sum(term_zkz_powY)[1]

    θ_tmp = _weibull_k_to_λ(k_tmp, sum_term_zkz, sum_term_zkz_powY, penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    obj = log(k_tmp)*sum_term_zkz - k_tmp*log(θ_tmp)*sum_term_zkz + (k_tmp-1)*sum_term_zkz_logY - (θ_tmp^(-k_tmp))*sum_term_zkz_powY
    p = penalty ? (pen_pararms_jk[1]-1)*log(k_tmp) - k_tmp/pen_pararms_jk[2] + (pen_pararms_jk[3]-1)*log(θ_tmp) - θ_tmp/pen_pararms_jk[4] : 0.0
    return (obj + p) * (-1.0)

end
function EM_M_expert_exact(d::WeibullExpert,
                            ye, exposure,
                            z_e_obs;
                            penalty = true, pen_pararms_jk = [Inf 1.0 Inf])

    # Further E-Step
    logY_e_obs = log.(ye)
    nan2num(logY_e_obs, 0.0) # get rid of NaN

    # Update parameters
    pos_idx = (ye .!= 0.0)
    logk_new = Optim.minimizer( Optim.optimize(x -> _weibull_optim_k_exact(x,
                                                            ye[pos_idx],
                                                            z_e_obs[pos_idx],
                                                            logY_e_obs[pos_idx],
                                                            penalty = penalty, pen_pararms_jk = pen_pararms_jk),
                                               max(log(d.k)-2.0, 0.0), log(d.k)+2.0 ))

    k_new = exp(logk_new)

    # Further E-step
    powY_e_obs = ye.^(k_new)
    nan2num(powY_e_obs, 0.0)

    term_zkz = z_e_obs
    term_zkz_powY = z_e_obs .* powY_e_obs

    sum_term_zkz = sum(term_zkz)[1]
    sum_term_zkz_powY = sum(term_zkz_powY)[1]

    θ_new = _weibull_k_to_λ(k_new, sum_term_zkz, sum(sum_term_zkz_powY)[1], penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    return WeibullExpert(k_new, θ_new)

end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
const RL = CommonRLInterface

@testset "MDP to Env" begin
    struct RLTestMDP <: MDP{Int, Int} end

    POMDPs.actions(m::RLTestMDP) = [-1, 1]
    POMDPs.transition(m::RLTestMDP, s, a) = Deterministic(clamp(s + a, 1, 3))
    POMDPs.initialstate(m::RLTestMDP) = Deterministic(1)
    POMDPs.isterminal(m::RLTestMDP, s) = s == 3
    POMDPs.reward(m::RLTestMDP, s, a, sp) = sp

    env = convert(RL.AbstractEnv, RLTestMDP())

    @test !RL.provided(RL.observations, env)
    POMDPs.states(m::RLTestMDP) = 1:3
    @test RL.provided(RL.observations, env)
    @test all(only.(RL.observations(env)) .== states(RLTestMDP()))

    @test RL.actions(env) == [-1, 1]
    @test RL.valid_actions(env) == [-1, 1]
    @test RL.observe(env) == [1]
    @test RL.state(env) == [1]
    @test RL.act!(env, 1) == 2
    @test !RL.terminated(env)

    RL.reset!(env)
    @test RL.observe(env) == [1]
    @test RL.act!(env, 1) == 2
    @test RL.act!(env, 1) == 3
    @test RL.observe(env) == [3]
    @test RL.terminated(env)

    RL.setstate!(env, [2])
    @test RL.observe(env) == [2]

    env2 = RL.clone(env)
    @test RL.act!(env2, 1) == 3
    @test RL.observe(env2) == [3]
    @test RL.observe(env) == [2]

    @test convert(MDP, env) isa RLTestMDP
end

@testset "POMDP to Env" begin
    struct RLTestPOMDP <: POMDP{Int, Int, Int} end

    POMDPs.actions(m::RLTestPOMDP) = [-1, 1]
    POMDPs.states(m::RLTestPOMDP) = 1:3
    POMDPs.transition(m::RLTestPOMDP, s, a) = Deterministic(clamp(s + a, 1, 3))
    POMDPs.observation(m::RLTestPOMDP, s, a, sp) = Deterministic(sp + 1)
    POMDPs.initialstate(m::RLTestPOMDP) = Deterministic(1)
    POMDPs.initialobs(m::RLTestPOMDP, s) = Deterministic(s + 1)
    POMDPs.isterminal(m::RLTestPOMDP, s) = s == 3
    POMDPs.reward(m::RLTestPOMDP, s, a, sp) = sp

    env = convert(RL.AbstractEnv, RLTestPOMDP())

    @test !RL.provided(RL.observations, env)
    POMDPs.observations(m::RLTestPOMDP) = 2:4
    @test RL.provided(RL.observations, env)
    @test all(only.(RL.observations(env)) .== observations(RLTestPOMDP()))

    @test RL.actions(env) == [-1, 1]
    @test RL.valid_actions(env) == [-1, 1]
    @test RL.observe(env) == [2]
    @test RL.state(env) == (1,2)
    RL.setstate!(env, RL.state(env))
    @test RL.act!(env, 1) == 2
    @test !RL.terminated(env)

    RL.reset!(env)
    @test RL.observe(env) == [2]
    @test RL.act!(env, 1) == 2
    @test RL.act!(env, 1) == 3
    @test RL.observe(env) == [4]
    @test RL.terminated(env)

    RL.setstate!(env, (2,3))
    @test RL.observe(env) == [3]

    env2 = RL.clone(env)
    @test RL.act!(env2, 1) == 3
    @test RL.observe(env2) == [4]
    @test RL.observe(env) == [3]

    @test convert(POMDP, env) isa RLTestPOMDP
end

@testset "Env to MDP" begin
    mutable struct MDPEnv <: RL.AbstractEnv
        s::Int
    end

    RL.reset!(env::MDPEnv) = env.s = 1
    RL.actions(env::MDPEnv) = [-1, 1]
    RL.observe(env::MDPEnv) = [env.s]
    RL.terminated(env::MDPEnv) = env.s >= 3
    function RL.act!(env::MDPEnv, a)
        r = env.s
        env.s = max(1, env.s + a)
        return r
    end

    m1 = convert(MDP, MDPEnv(1))
    @test m1 isa OpaqueRLEnvMDP
    @test simulate(RolloutSimulator(), m1, FunctionPolicy(s->1)) == 3.0

    RL.@provide RL.state(env::MDPEnv) = env.s
    RL.@provide RL.setstate!(env::MDPEnv, s) = env.s = s

    env = MDPEnv(1)
    m2 = convert(MDP, env)
    @test m2 isa RLEnvMDP
    @test simulate(RolloutSimulator(), m2, FunctionPolicy(s->1)) == 3.0

    RL.setstate!(env, 1)
    @test isterminal(m2, 3)
    @test RL.state(env) == 1

    @test convert(RL.AbstractEnv, m1) isa MDPEnv
end

@testset "Env to POMDP" begin
    mutable struct POMDPEnv <: RL.AbstractEnv
        s::Int
    end

    RL.reset!(env::POMDPEnv) = env.s = 1
    RL.actions(env::POMDPEnv) = [-1, 1]
    RL.observe(env::POMDPEnv) = [env.s > 0]
    RL.terminated(env::POMDPEnv) = env.s >= 3
    function RL.act!(env::POMDPEnv, a)
        r = env.s
        env.s = max(1, env.s + a)
        return r
    end

    m1 = convert(POMDP, POMDPEnv(1))
    @test m1 isa OpaqueRLEnvPOMDP
    @test simulate(RolloutSimulator(), m1, FunctionPolicy(s->1)) == 3.0

    RL.@provide RL.state(env::POMDPEnv) = env.s
    RL.@provide RL.setstate!(env::POMDPEnv, s) = env.s = s

    env = POMDPEnv(1)
    m2 = convert(POMDP, env)
    @test m2 isa RLEnvPOMDP
    @test simulate(RolloutSimulator(), m2, FunctionPolicy(s->1)) == 3.0

    RL.setstate!(env, 1)
    @test isterminal(m2, 3)
    @test RL.state(env) == 1

    @test convert(RL.AbstractEnv, m1) isa POMDPEnv
end


text:
for ( julia_fun, ippf_prefix )  in  [ (   :"insert julia function name",  "ippsFunctionBaseName"  ) ]
    for ( "TypeSignatures" ) in "AnArrayOfTuples"

        julia_fun! = symbol(string(julia_fun, '!')) # in-place version if it makes sense

        @eval begin
            function $(julia_fun!)( buffer, "julia function arguments" )
                sigLen = length( signal )
                outLen = length( buffer )
                @ippscall( $ippfsr,  (  "the c function argument types" ),
                                        "julia varables to pass to the c function")
                buffer
            end # function

            $(julia_fun)( "julia function arguments" ) = $(julia_fun!)( "create a buffer", "julia function arguments" )

        end # eval
    end # type loop
end # function loop


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

mutable struct AccessUriRaw <: SwaggerModel
    accessSAS::Any # spec type: Union{ Nothing, String } # spec name: accessSAS

    function AccessUriRaw(;accessSAS=nothing)
        o = new()
        validate_property(AccessUriRaw, Symbol("accessSAS"), accessSAS)
        setfield!(o, Symbol("accessSAS"), accessSAS)
        o
    end
end # type AccessUriRaw

const _property_map_AccessUriRaw = Dict{Symbol,Symbol}(Symbol("accessSAS")=>Symbol("accessSAS"))
const _property_types_AccessUriRaw = Dict{Symbol,String}(Symbol("accessSAS")=>"String")
Base.propertynames(::Type{ AccessUriRaw }) = collect(keys(_property_map_AccessUriRaw))
Swagger.property_type(::Type{ AccessUriRaw }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_AccessUriRaw[name]))}
Swagger.field_name(::Type{ AccessUriRaw }, property_name::Symbol) =  _property_map_AccessUriRaw[property_name]

function check_required(o::AccessUriRaw)
    true
end

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


text:
#get dictionary for storing log information
function getLog(m0,nc,maxit)
	 data=Dict{Symbol, Any}()
	 data[:residuals]=zeros(Float64,maxit,m0)
	 data[:feastResidual]=zeros(Float64,maxit)
	 data[:eigenvalues]=zeros(ComplexF64,maxit,m0)
	 data[:ninside]=zeros(Int64,maxit)
	 data[:iterations]=0
	 data[:insideIndices]=Dict{Int64, Any}()
	 data[:shiftIndex]=Dict{ComplexF64,Int64}()

	 data[:linResiduals]=zeros(Float64,m0,nc,maxit)
	 data[:linIts]=zeros(Int64,m0,nc,maxit)

	 data[:hlinResiduals]=zeros(Float64,m0,nc,maxit)
	 data[:hlinIts]=zeros(Int64,m0,nc,maxit)

	 return data
end

########################
# CORE FEAST FUNCTIONS #
########################

#These functions are one- and two-sided implementations of the FEAST algorithm. They are essentially outlines that encompass all possible variations of the FEAST algorithm.

#They are used by passing functions to them to perform the essential operations of the FEAST algorithm, such as evaluating the eigenvector residual and calculating the contour integrand.

#core one-sided FEAST algorithm; this function is used for all variations of FEAST
function feast_core(Tf,integrand,rrsolve,x0,nc,emid,ra,rb,eps::Float64,maxit; log=false, insideEps=1e-2, verbose=1, set_nin=false)
	 #Tf: eigenvector residual function. E.g. for linear problems Tf(l,x) = B*x*diagm(l)-A*x
	 #integrand: function that returns integral integrand evaluated at point z in complex plane, applied to a vector x. E.g. for linear FEAST integrand(z,x,l)=\(z*B-A,x)
	 #rrsolve: function for returning eigenvalue and eigenvector approximations from rayleigh ritz. E.g. for linear FEAST rrsolve(Q)=eig(Q'*A*Q,Q'*B*Q)
	 #x0: initial guess for eigenvectors
	 #nc: number of quadrature points for numerical integration
	 #emid: center of integration ellipse
	 # ra,rb: radii of integration ellipse
	 # eps: convergence tolerance
	 # maxit: maximum number of FEAST iterations
	 # insideEps: residual threshold for determining whether or not an eigenvalue inside the contour is spurious
	 # log: if true, returns dictionary object with logged convergence information
	 # verbose(0-1): how much information to print out
	 # set_nin: if equal to an integer, sets the number of eigenvalues inside the contour to 'set_nin'; used for diagnostic purposes

	 #get shapes of arrays:
	(n,m0)=size(x0)

	 #storing convergence info:
	 data=getLog(m0,nc,maxit)

	 #data[:linsysresidual]=zeros(Float64,1)
	 #data[:linsysits]=zeros(Int64,1)

	 #start with initial guess
	 x=copy(x0)

	 #initialize eigenvalue estimates:
	lest=zeros(ComplexF64,m0,1)

	#initialize iteration number:
	it=0

	 #get quadrature points from trapezoidal rule:
	 (gk,wk)=trapezoidal(nc)
	 offset=pi/nc #offset angle of first quadrature point; makes sure it isn't on real axis for hermitian problems

	 #save indices for contour points
	 for k in 1:nc
		  #curve parametrization angle:
		  theta=gk[k]*pi+pi+offset

		  #quadrature point:
		  z=emid+((ra+rb)/2)*exp(im*theta)+((ra-rb)/2)*exp(-1.0*im*theta)
		  data[:shiftIndex][z]=k
	 end

	 #Initialize FEAST subspace:
	 Q=copy(x)

	 res=1.0 #initial residual
	while res>eps && it<maxit
		it=it+1 #update number of iterations

		  data[:iterations]=it

		  #rayleigh ritz
		  (lest,x)=rrsolve(Q)

		  #sort everything by real part of eigenvalue
		  p=sortperm(real(lest))
		  lest=lest[p]
		  x=x[:,p]

		  #store convergence data
		  #indicate which eigenvalues are inside the contour
		  data[:insideIndices][it]=getInsideIndex(lest,emid,ra,rb)
		  #store all the eigenvalues
		  data[:eigenvalues][it,:]=lest

		  #calculate eigenvector residuals
		  resvecs=Tf(lest,x)
		  for i in 1:m0
				data[:residuals][it,i]=norm(resvecs[:,i])/norm(x[:,i])
		  end

		  #find the largest residual inside the contour
		  res,ninside=getMaxResInside(lest,x,resvecs,emid,ra,rb;inEps=insideEps, set_nin=set_nin)

		  data[:feastResidual][it]=res
		  data[:ninside][it]=ninside

		  verbose>0 && println("  $it  res=$res  minres=$(minimum(data[:residuals][it,:]))  nin=$ninside")

		  if(res<eps)
				if(log)
					return (lest,x,data)
			  else
					return (lest,x)
			  end
		  end

		  #apply contour integral to get FEAST subspace:
		 Q=zero(x)
		 for k in 1:nc
		 #Q=@parallel (+) for k in 1:nc
			  #integration curve is an ellipse centered at emid, with radii ra and rb

			  #curve parametrization angle:
			  theta=gk[k]*pi+pi+offset

			  #quadrature point:
			  z=emid+((ra+rb)/2)*exp(im*theta)+((ra-rb)/2)*exp(-1.0*im*theta)

				#integrand evaluated at quadrature point z

				#println("     linear system $k")

				Qk=integrand(z,x,lest,data,resvecs)

				#add integrand contribution to quadrature:
			  Q=Q+wk[k]*(((ra+rb)/2)*exp(im*theta)-((ra-rb)/2)*exp(-1.0*im*theta))*Qk
			  #    wk[k]*(((ra+rb)/2)*exp(im*theta)-((ra-rb)/2)*exp(-1.0*im*theta))*Qk
		  end

		  #Orthonormalize FEAST subspace to avoid spurious eigenpairs:
		  #0.63
		  #(Qq,Rq)=qr(Q)
		  #Q[:]=Qq

		  #0.7
		  Q=Array(qr(Q).Q)
	end

	 if(log)
		  return (lest,x,data)
	 else
		  return (lest,x)
	 end

end

#core two-sided FEAST algorithm
#Difference from symmetric FEAST: have to solve for left and right eigenvectors simultaneously and biorthogonalize them
function feastNS_core(Tf,hTf,integrand,hintegrand,factorizer,hfactorizer,FactoType,rrsolvens,biortho,x0,y0,nc,emid,ra,rb,eps,maxit; insideEps=1e-2, log=false)
	 #Tf: eigenvector residual function. E.g. for linear problems Tf(l,x) = B*x*diagm(l)-A*x
	 #hTf: hermitian conjugate transpose of residual function
	 #integrand: function that returns integral integrand evaluated at point z in complex plane, applied to a vector x. E.g. for linear FEAST integrand(z,x,l)=\(z*B-A,x)
	 #rrsolvens: function for returning eigenvalue and eigenvector approximations from rayleigh ritz. E.g. for linear FEAST rrsolve(Q)=eig(Q'*A*Q,Q'*B*Q)
	 #x0: initial guess for right eigenvectors
	 #y0: initial guess for left eigenvectors
	 #nc: number of quadrature points for numerical integration
	 #emid: center of integration ellipse
	 # ra,rb: radii of integration ellipse
	 # eps: convergence tolerance
	 # maxit: maximum number of FEAST iterations
	 # insideEps: residual threshold for determining whether or not an eigenvalue inside the contour is spurious
	 # log: if true, returns dictionary object with logged convergence information

	 #get shapes of arrays:
	(n,m0)=size(x0)

	#storing convergence info:
	 data=getLog(m0,nc,maxit)

	 #start with initial guess
	 x=copy(x0)
	 y=copy(y0)

	 #initialize eigenvalue estimates:
	lest=zeros(ComplexF64,m0,1)

	#initialize iteration number:
	it=0

	 #get quadrature points from trapezoidal rule:
	 (gk,wk)=trapezoidal(nc)
	 offset=pi/nc #offset angle of first quadrature point; makes sure it isn't on real axis for hermitian problems

	 #save indices for contour points
	 data[:storedFacto] = Vector{FactoType}(undef, nc)
	 data[:storedFactoH] = Vector{FactoType}(undef, nc)
	 print("Factorizing: ")
	 for k in 1:nc
		  #curve parametrization angle:
		  theta=gk[k]*pi+pi+offset

		  #quadrature point:
		  z=emid+((ra+rb)/2)*exp(im*theta)+((ra-rb)/2)*exp(-1.0*im*theta)
		  data[:shiftIndex][z] = k
		  data[:storedFacto][k] = factorizer(z)
		  data[:storedFactoH][k] = hfactorizer(z)
		  print(k, " ")
	 end
	 println()

	 #Initialize FEAST subspaces:
	 Q=copy(x)
	 R=copy(y)

	 res=1.0 #initial residual
	while res>eps && it<maxit
		it=it+1 #update number of iterations
		data[:iterations]=it

		  #Biorthogonalize Q and R subspaces
		  println("Biorthogonalization")
		  (Q,R)=biortho(Q,R)

		  #rayleigh ritz
		  println("Rayleigh-Ritz")
		  (y,lest,x)=rrsolvens(Q,R)

		  #calculate eigenvector residuals
		  println("Left & Right Residuals")
		  resvecs=Tf(lest,x)
		  hresvecs=hTf(lest,y)

		  println("Convergence data")
		  #store convergence data
		  #indicate which eigenvalues are inside the contour
		  data[:insideIndices][it]=getInsideIndex(lest,emid,ra,rb)
		  #store all the eigenvalues
		  data[:eigenvalues][it,:]=lest

		  println("Right Residuals (duplicate)")
		  #calculate eigenvector residuals
		  resvecs=Tf(lest,x)
		  for i in 1:m0
				data[:residuals][it,i]=norm(resvecs[:,i])/norm(x[:,i])
		  end

		  println("Inside/outside")
		  #find the largest residual inside the contour
		  res,ninside=getMaxResInside(lest,x,resvecs,emid,ra,rb; inEps=insideEps)
		  data[:feastResidual][it]=res
		  data[:ninside][it]=ninside

		  println("  $it  res=$res  minres=$(minimum(data[:residuals][it,:]))  nin=$ninside")
		  if(res<eps)
				if(log)
						  return (y,lest,x,data)
					 else
						  return (y,lest,x)
				end
		  end

		  print("Contour integral: ")
		  #apply contour integral to get FEAST subspace:
		 Q=zero(x)
		 R=zero(y)
		 for k in 1:nc
			  #integration curve is an ellipse centered at emid, with radii ra and rb
			  print(k," ")
			  #curve parametrization angle:
			  theta=gk[k]*pi+pi+offset

			  #quadrature point:
			  z=emid+((ra+rb)/2)*exp(im*theta)+((ra-rb)/2)*exp(-1.0*im*theta)

				#integrand evaluated at quadrature point z
				Qk=integrand(data[:storedFacto][k],x,lest,data,resvecs)
				Rk=hintegrand(data[:storedFactoH][k],y,lest,data,hresvecs)

				#add integrand contribution to quadrature:
			  Q=Q+wk[k]*(((ra+rb)/2)*exp(im*theta)-((ra-rb)/2)*exp(-1.0*im*theta))*Qk
			  R=R+(wk[k]*(((ra+rb)/2)*exp(im*theta)-((ra-rb)/2)*exp(-1.0*im*theta)))'*Rk
		  end
		  println()
	end

		  if(log)
				return (y,lest,x,data)
		  else
				return (y,lest,x)
		  end

end

max_stars_repo_name: rleegates/feastjl

max_issues_repo_name: rleegates/feastjl

max_forks_repo_name: rleegates/feastjl


text:
# Empty Stream
# ============

"""
`EmptyStream` is a dummy stream to allow a buffered stream to wrap an
array without additional buffering.
"""
struct EmptyStream end

Base.close(::EmptyStream) = nothing

Base.position(stream::BufferedInputStream{EmptyStream}) = stream.position - 1
Base.position(stream::BufferedOutputStream{EmptyStream}) = stream.position - 1

# Buffered input stream
# ---------------------

function BufferedInputStream(data::Vector{UInt8}, len::Integer = lastindex(data))
    return BufferedInputStream{EmptyStream}(EmptyStream(), data, 1, len, 0, false)
end

function fillbuffer!(stream::BufferedInputStream{EmptyStream})
    return 0
end

function Base.eof(stream::BufferedInputStream{EmptyStream})
    return stream.position > stream.available
end

function Base.read(stream::BufferedInputStream{EmptyStream}, ::Type{UInt8})
    checkopen(stream)
    if eof(stream)
        throw(EOFError())
    end
    byte = stream.buffer[stream.position]
    stream.position += 1
    return byte
end

function Base.seek(stream::BufferedInputStream{EmptyStream}, pos::Integer)
    checkopen(stream)
    upanchor!(stream)
    if 1 <= pos + 1 <= stream.available
        stream.position = pos + 1
    else
        throw(BoundsError)
    end
end

# Buffered output stream
# ----------------------

function BufferedOutputStream()
    return BufferedOutputStream(Vector{UInt8}(undef, 1024))
end

function BufferedOutputStream(data::Vector{UInt8})
    return BufferedOutputStream{EmptyStream}(EmptyStream(), data, 1)
end

function flushbuffer!(stream::BufferedOutputStream{EmptyStream}, eof::Bool=false)
    if available_bytes(stream) == 0
        sz = 2 * length(stream.buffer)
        sz = sz > 0 ? nextpow(2, sz) : 0
        resize!(stream.buffer, max(sz, 16))
    end
    return
end

function Base.write(stream::BufferedOutputStream{EmptyStream}, byte::UInt8)
    checkopen(stream)
    flushbuffer!(stream)
    stream.buffer[stream.position] = byte
    stream.position += 1
    return 1
end

function Base.write(stream::BufferedOutputStream{EmptyStream}, data::Vector{UInt8})
    checkopen(stream)
    n = length(data)
    if n > available_bytes(stream)
        sz =
        resize!(stream.buffer, nextpow(2, n + stream.position - 1))
    end
    copyto!(stream.buffer, stream.position, data, 1)
    stream.position += n
    return n
end

# Faster append for vector-backed output streams.
@inline function Base.append!(stream::BufferedOutputStream{EmptyStream},
                              data::Vector{UInt8}, start::Int, stop::Int)
    n = stop - start + 1
    if stream.position + n > length(stream.buffer)
        resize!(stream.buffer, max(n, max(1024, 2 * length(stream.buffer))))
    end
    copyto!(stream.buffer, stream.position, data, start, n)
    stream.position += n
end

function Base.take!(stream::BufferedOutputStream{EmptyStream})
    # TODO: benchmark resizing stream.buffer, returning it, and replacing it
    # with an zero-length array, which might be fast in the common case of
    # building just one array/string.
    chunk = stream.buffer[1:stream.position-1]
    stream.position = 1
    return chunk
end

function Base.empty!(stream::BufferedOutputStream{EmptyStream})
    return stream.position = 1
end

function Base.isempty(stream::BufferedOutputStream{EmptyStream})
    return stream.position == 1
end

function Base.length(stream::BufferedOutputStream{EmptyStream})
    return stream.position - 1
end

function Base.:(==)(a::BufferedOutputStream{EmptyStream}, b::BufferedOutputStream{EmptyStream})
    if a.position == b.position
        return ccall(:memcmp, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
                     a.buffer, b.buffer, a.position - 1) == 0
    else
        return false
    end
end


text:
using Plots
U(r,eps,sig) = 4*eps*((sig/r)^12 - (sig/r)^6)
# r vector between 0.01 and 15.
r = collect(0.4:0.01:2.5)
# compute U for every r
eps = 2. ; sig = 0.5 ;
u = [ U(x,eps,sig) for x in r ]
plot(r,u,label="",xlabel="r",ylabel="U(r)",
     ylim=[-2.5,10],framestyle=:box,linewidth=2,
     size=(400,300))
savefig("./cutoff1.pdf")


text:
using Test
using ModelObjectsLH
using CollegeEntry

ce = CollegeEntry;

function adm_prob_test(switches)
    @testset "$switches" begin
        af = init_admprob_fct(switches);
        @test ce.validate_admprob_fct(af);
        @test get_object_id(af) isa ObjectId;
        @test get_object_id(switches) isa ObjectId;

        if !isnothing(ce.param_names(af))
            pNameV = ce.param_names(af);
            if !isnothing(ce.param_names(af))
                for pName in pNameV
                    ce.by_college!(switches, pName);
                    @test ce.by_college(switches, pName);
                    ce.not_by_college!(switches, pName);
                    @test !ce.by_college(switches, pName);
                end
            end
        end

        xV = 0.01 : 0.1 : 0.99;
        for ic = 1 : n_colleges(af)
            f = make_admprob_function(af, ic);
            probV = f.(xV);
            @test size(probV) == size(xV);
            @test all(0.0 .<= probV .<= 1.0);
            @test all(diff(probV) .>= 0.0);

            prob2V = prob_admit(af, ic, xV);
            @test isapprox(probV, prob2V);
        end
    end
end

@testset "Admission probs" begin
    nc = 4;
    for switches in (
        ce.make_test_admprob_fct_step_switches(nc),
        ce.make_test_admprob_fct_logistic_switches(nc),
        ce.make_test_admprob_fct_linear_switches(nc),
        AdmProbFctOpenSwitches{Float64}(ObjectId(:test), nc)
        )
        adm_prob_test(switches);
    end
end

# ----------------


text:
#= 𝐵-effective =#
# Methods for getting 𝐵-effective
export rfgr2B
"""
    B = rfgr2B(rf, gr, loc=[0 0 0]u"cm"; Δf=0u"Hz", b1Map=1, γ=γ¹H)
Turn rf, `rf`, and gradient, `gr`, into 𝐵-effective magnetic field.

*INPUTS*:
- `rf::TypeND(RF0D, [1,2])` (nT, (nCoil))
- `gr::TypeND(GR0D, [2])`   (nT, 3)
- `loc::TypeND(L0D, [2])`   (1,3) or (nM, 3), locations.
*KEYWORDS*:
- `Δf::TypeND(F0D, [0,1,2])` (1,)  or (nM,), off-resonance.
- `b1Map::TypeND(Union{Real,Complex},[0,1,2])` (1,) or (nM,(nCoils)),
   transmit sensitivity.
- `γ::TypeND(Γ0D, [0,1])` (1,)  or (nM,), gyro-ratio
*OUTPUS*:
- `B`, generator of `TypeND(B0D, [2])` (1,1,nT), 𝐵 field.

See also: [`Pulse2B`](@ref), [`blochSim`](@ref).

# TODO:
Support `loc`, `Δf`, and `b1Map` being `Base.Generators`.
"""
function rfgr2B(rf   ::TypeND(RF0D, [1,2]),
                gr   ::TypeND(GR0D, [2]),
                loc  ::TypeND(L0D,  [2]) = [0 0 0]u"cm";
                Δf   ::TypeND(F0D,  [0,1]) = 0u"Hz",
                b1Map::TypeND(Union{Real,Complex}, [0,1,2,3]) = 1,
                γ    ::TypeND(Γ0D,  [0,1]) = γ¹H)

  nM = maximum(map(x->size(x,1), (loc, Δf, b1Map, γ)))

  Bxy_gen = b1Map==1 ?
    @inbounds(view(rf,t,:)       |> x->repeat([real(x) imag(x)], nM)
              for t in axes(rf,1)) :
    @inbounds(b1Map*view(rf,t,:) |> x->       [real(x) imag(x)]
              for t in axes(rf,1))

  Bz_gen = Δf==0u"Hz" ?
    @inbounds(loc*view(gr,t,:)          for t in axes(gr,1)) :
    @inbounds(loc*view(gr,t,:).+(Δf./γ) for t in axes(gr,1))

  B_gen = @inbounds([bxy bz] for (bxy, bz) in zip(Bxy_gen, Bz_gen))
  return B_gen
end

#= rotation axis/angle, U/Φ =#
export B2UΦ, B2UΦ!
"""
    B2UΦ(B::TypeND(B0D,[2,3]); γ::TypeND(Γ0D,[0,1]), dt::T0D=4e-6u"s")
Given 𝐵-effective, `B`, compute rotation axis/angle, `U`/`Φ`.

*INPUTS*:
- `B::TypeND(B0D, [2,3])` (1,3,nT) or (nM, 3, nT), 𝐵 field.
*KEYWORDS*:
- `γ::TypeND(Γ0D, [0,1])`: Global, (1,); Spin-wise, (nM, 1). gyro ratio
- `dt::T0D` (1,), simulation temporal step size, i.e., dwell time.
*OUTPUTS*:
- `U::TypeND(Real, [2,3])` (1,3,(nT)) or (nM,3,(nT)), axis.
- `Φ::TypeND(Real, [2,3])` (1,1,(nT)) or (nM,1,(nT)), angle.

See also: [`B2UΦ!`](@ref), [`UΦRot!`](@ref).

# Notes:
Somehow, in-place version, `B2UΦ!(B,U,Φ; γ,dt)`, provokes more allocs in julia.
"""
@inline function B2UΦ(B::TypeND(B0D, [2,3]);
                      γ::TypeND(Γ0D, [0,1]), dt::T0D=4e-6u"s")
  X = γ*dt
  U = ustrip.(Float64, unit.(X).^-1, B)
  Φ = sqrt.(sum(U.*U, dims=2))
  U .= U./Φ .|> x->isnan(x) ? 0 : x
  Φ .*= -ustrip.(X).*2π # negate to make: 𝐵×𝑀 → 𝑀×𝐵
  return (U=U, Φ=Φ)
end

"""
    B2UΦ!(B, U; Φ, γ, dt=(4e-6)u"s")
In-place version of `B2UΦ`. Somehow, `B2UΦ!`, provokes more allocs in julia.

See also: [`B2UΦ`](@ref), [`blochSim`](@ref).
"""
@inline function B2UΦ!(B::TypeND(B0D, [2,3]),
                       U::TypeND(Float64, [2,3]),
                       Φ::TypeND(Float64, [2,3]);
                       γ::TypeND(Γ0D, [0,1]), dt::T0D=(4e-6)u"s")
  X = γ*dt
  U .= ustrip.(Float64, unit.(X).^-1, B)
  Φ .= sqrt.(sum(U.*U, dims=2))
  U .= U./Φ .|> x->isnan(x) ? 0 : x # negate to make: 𝐵×𝑀 → 𝑀×𝐵
  Φ .*= -ustrip.(X).*2π
  return (U=U, Φ=Φ)
end

export UΦRot!, UΦRot
"""
    UΦRot!(U, Φ, V, R)
Apply axis-angle, `U`-`Φ` based rotation on `V`. Rotation is broadcasted on `V`
along its 3rd dimension. Results will overwrite into `R`.

*INPUTS*:
- `U::TypeND(AbstractFloat,[2])` (nM, 3), rotation axes in 3D, assumed unitary;
- `Φ::TypeND(AbstractFloat,[1])` (nM,), rotation angles;
- `V::TypeND(AbstractFloat,[2,3])` (nM, 3, (3)), vectors to be rotated;
- `R::TypeND(AbstractFloat,[2,3])` (nM, 3, (3)), vectors rotated, i.e., results;
*OUTPUTS*:
- `R` the input container `R` is also returned for convenience.

See also: [`UΦRot`](@ref), [`B2UΦ!`](@ref).
"""
@inline function UΦRot!(U::TypeND(AbstractFloat,[2]),
                        Φ::TypeND(AbstractFloat,[1]),
                        V::TypeND(AbstractFloat,[2,3]),
                        R::TypeND(AbstractFloat,[2,3]))
  # en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
  # 𝑅 = 𝑐𝑜𝑠𝜃⋅𝐼 - (𝑐𝑜𝑠𝜃-1)⋅(𝐮𝐮ᵀ) + 𝑠𝑖𝑛𝜃⋅[𝐮]ₓ; 𝜃/𝐮, rotation angle/axis
  cΦ, sΦ = cos.(Φ), sin.(Φ)
  (Vx,Vy,Vz) = map(x->size(V,3)==1 ? view(V,:,[x]) : view(V,:,[x],:), (1,2,3))
  (Ux,Uy,Uz) = map(x->view(U,:,x), (1,2,3))

  R .= cΦ.*V .+ (1 .- cΦ).*sum(U.*V,dims=2).*U .+
       sΦ.*hcat(-Uz.*Vy.+Uy.*Vz, Uz.*Vx.-Ux.*Vz, -Uy.*Vx.+Ux.*Vy)
  return R
end

"""
    UΦRot(U, Φ, V)
Same as `UΦRot!(U, Φ, V, R)`, except not in-place.

See also: [`UΦRot!`](@ref).
"""
@inline UΦRot(U, Φ, V) = UΦRot!(U, Φ, V, copy(V))

#= 𝐴, 𝐵 =#
export B2AB
"""
    B2AB(B; T1=(Inf)u"s", T2=(Inf)u"s", γ=γ¹H, dt=(4e-6)u"s")
Turn B-effective into Hargreave's 𝐴/𝐵, mat/vec, see: doi:10.1002/mrm.1170.

*INPUTS*:
- `B::Union{TypeND(B0D, [2,3]), Base.Generator}`:
  Global, (nT,xyz); Spin-wise, (nM,xyz,nT).
*KEYWORDS*:
- `T1 & T2 ::TypeND(T0D, [0,1])`: Global, (1,); Spin-wise, (nM,1).
- `γ::TypeND(Γ0D, [0,1])`: Global, (1,); Spin-wise, (nM, 1). gyro ratio
- `dt::T0D` (1,), simulation temporal step size, i.e., dwell time.
*OUTPUTS*:
- `A::TypeND(AbstractFloat,[3])` (nM, 3,3), `A[iM,:,:]` is the `iM`-th 𝐴.
- `B::TypeND(AbstractFloat,[2])` (nM, 3), `B[iM,:]` is the `iM`-th 𝐵.

See also: [`rfgr2B`](@ref), [`Pulse2B`](@ref).
"""
function B2AB(B ::Base.Generator;
              T1::TypeND(T0D, [0,1])=(Inf)u"s",
              T2::TypeND(T0D, [0,1])=(Inf)u"s",
              γ ::TypeND(Γ0D, [0,1])=γ¹H,
              dt::T0D=(4e-6)u"s")

  nM = maximum([size(x,1) for x in (T1,T2,γ,first(B))])
  AB = reshape([ones(nM) zeros(nM,3) ones(nM) zeros(nM,3) ones(nM) zeros(nM,3)],
               (nM,3,4)) # as if cat(A,B;dims=3), avoid constructing U, Φ twice.
  AB1 = similar(AB)

  # in unit, convert relaxations into losses/recovery per step
  E1, E2 = exp.(-dt./T1), exp.(-dt./T2)
  E1₋₁ = E1 .- 1

  # u, ϕ = Array{Float64}(undef, nM, 3), Array{Float64}(undef, nM, 1)

  for b in B
    u, ϕ = B2UΦ(b; γ=γ, dt=dt)
    # B2UΦ!(b, u, ϕ; γ=γ, dt=dt)
    any(ϕ.!=0) && @inbounds(UΦRot!(u, view(ϕ,:,1), AB, AB1))
    AB1[:,1:2,:] .*= E2
    AB1[:,3,:]   .*= E1
    AB1[:,3,4]   .-= E1₋₁
    AB, AB1 = AB1, AB
  end

  return (A=AB[:,:,1:3], B=AB[:,:,4]) # Can one avoid the array copies here?
end

function B2AB(B ::TypeND(B0D, [2,3]); kw...)
  if size(B,3) == 1 && size(B,1) != size(M,1)
    B = permutedims(B[:,:,:], [3,2,1])  # best practice?
    println("B not being spin-specific, assuming global")
  end
  return B2AB(@inbounds(view(B,:,:,t) for t in axes(B,3)); kw...)
end

#= blochSim =#
export blochSim!, blochSim
"""
    blochSim!(M, B; T1=(Inf)u"s",T2=(Inf)u"s",γ=γ¹H,dt=(4e-6)u"s",doHist=false)
Old school 𝐵-effective magnetic field, `B`, based bloch simulation. Globally or
spin-wisely apply `B` over spins, `M`. `M` will be updated by the results.

*INPUTS*:
- `M::TypeND(Real, [2])` (nM, xyz): input spins' magnetizations.
- `B::Union{TypeND(B0D, [2,3]), Base.Generator}`:
  Global, (nT,xyz); Spin-wise, (nM,xyz,nT).
*KEYWORDS*:
- `T1 & T2 ::TypeND(T0D, [0,1])`: Global, (1,); Spin-wise, (nM,1).
- `γ::TypeND(Γ0D, [0,1])`: Global, (1,); Spin-wise, (nM, 1). gyro ratio
- `dt::T0D` (1,), simulation temporal step size, i.e., dwell time.
- `doHist::Bool`, whether to output spin history through out `B`.
*OUTPUTS*:
- `M::TypeND(Real, [2])` (nM, xyz): spins after applying `B`.
- `Mhst::TypeND(Real, [3])` (nM, xyz, nT): spins history during `B`.

See also: [`applyPulse`](@ref), [`blochSim`](@ref).

# Notes:
1. Not much sanity check inside this function, user is responsible for
   matching up the dimensions.
2. Put relax at the end of each time step may still be inaccurate, since
   physically spins relax continuously, this noise/nuance may worth study
   for applications like fingerprinting simulations, etc.
"""
function blochSim!(M ::TypeND(AbstractFloat, [2]),
                   B ::Base.Generator;
                   T1::TypeND(T0D, [0,1])=(Inf)u"s",
                   T2::TypeND(T0D, [0,1])=(Inf)u"s",
                   γ ::TypeND(Γ0D, [0,1])=γ¹H,
                   dt::T0D=(4e-6)u"s", doHist=false)

  # in unit, convert relaxations into losses/recovery per step
  E1, E2 = exp.(-dt./T1), exp.(-dt./T2)
  E1₋₁ = E1 .- 1

  nM = maximum([size(x,1) for x in (T1,T2,γ,first(B))])
  size(M,1) == 1 && (M = repeat(M, nM))

  Mhst = doHist ? zeros(size(M,1), 3, length(B)) : nothing
  Mi, M1 = M, similar(M) # Mi refers to the original array `M` in memory.

  # u, ϕ = Array{Float64}(undef, nM, 3), Array{Float64}(undef, nM, 1)

  if doHist
    for (t, b) in enumerate(B)
      u, ϕ = B2UΦ(b; γ=γ, dt=dt)
      # B2UΦ!(b, u, ϕ; γ=γ, dt=dt)
      any(ϕ.!=0) && @inbounds(UΦRot!(u, view(ϕ,:,1), M, M1))
      # relaxation
      M1[:,1:2] .*= E2
      M1[:,3]   .*= E1
      M1[:,3]   .-= E1₋₁
      @inbounds(Mhst[:,:,t] = M1)
      M, M1 = M1, M
    end
  else
    for b in B
      u, ϕ = B2UΦ(b; γ=γ, dt=dt)
      # B2UΦ!(b, u, ϕ; γ=γ, dt=dt)
      any(ϕ.!=0) && @inbounds(UΦRot!(u, view(ϕ,:,1), M, M1))
      # relaxation
      M1[:,1:2] .*= E2
      M1[:,3]   .*= E1
      M1[:,3]   .-= E1₋₁
      M, M1 = M1, M
    end
  end
  M === Mi || (Mi .= M) # if `M` doesn't point to the input array, update.

  return (M=M, Mhst=Mhst)
end

function blochSim!(M::TypeND(AbstractFloat, [2]), B::TypeND(B0D, [2,3]); kw...)
  if size(B,3) == 1 && size(B,1) != size(M,1)
    B = permutedims(B[:,:,:], [3,2,1])  # best practice?
    println("B not being spin-specific, assuming global")
  end
  return blochSim!(M, @inbounds(view(B,:,:,t) for t in axes(B,3)); kw...)
end

"""
    blochSim!(M, A, B)
Hargreave's 𝐴/𝐵, mat/vec, based bloch simulation. Globally or spin-wisely apply
matrix `A` and vector `B` over spins, `M`, described in doi:10.1002/mrm.1170

*INPUTS*:
- `M::TypeND(Real, [2])` (nM, xyz): input spins' magnetizations.
- `A::TypeND(AbstractFloat,[3])` (nM, 3,3), `A[iM,:,:]` is the `iM`-th 𝐴.
- `B::TypeND(AbstractFloat,[2])` (nM, 3), `B[iM,:]` is the `iM`-th 𝐵.
*OUTPUTS*:
- `M::TypeND(Real, [2])` (nM, xyz): spins after applying `B`.
"""
blochSim!(M, A, B) = M .= blochSim(M, A, B)

"""
    blochSim(M, B; T1, T2, γ, dt, doHist)
Same as `blochSim!(M, B; T1,T2,γ,dt,doHist)`, `M` will not be updated.

See also: [`blochSim!`](@ref)
"""
blochSim(M, B; kw...) = blochSim!(copy(M), B; kw...)

"""
    blochSim(M, A, B)
Same as `blochSim(M, A, B)`, `M` will not be updated.
"""
blochSim(M::TypeND(AbstractFloat,[2]),
         A::TypeND(AbstractFloat,[3]),
         B::TypeND(AbstractFloat,[2])) =
  sum(A.*permutedims(M[:,:,:],(1,3,2));dims=3) .+ B

# No inplace operation for `M::TypeND(Integer,[2])`.
blochSim(M::TypeND(Integer,[2]), a...; kw...) = blochSim(float(M), a...; kw...)

#= freePrec =#
export freePrec!, freePrec
"""
    freePrec!(M, t; Δf=0u"Hz", T1=(Inf)u"s", T2=(Inf)u"s")
Spins, `M`, free precess by time `t`. `M` will be updated by the results.

*INPUTS*:
- `M::TypeND(Real, [2])` (nM, xyz): input spins' magnetizations.
- `t::T0D` (1,): duration of free precession.
*KEYWORDS*:
- `T1 & T2 ::TypeND(T0D, [0,1])`: Global, (1,); Spin-wise, (nM,1).
*OUTPUTS*:
- `M::TypeND(Real, [2])` (nM, xyz): output spins' magnetizations.

See also: [`freePrec`](@ref).
"""
function freePrec!(M ::TypeND(AbstractFloat,[2]),
                   t ::T0D;
                   Δf::TypeND(F0D,[0,1])=0u"Hz",
                   T1::TypeND(T0D,[0,1])=(Inf)u"s",
                   T2::TypeND(T0D,[0,1])=(Inf)u"s")

  E1, E2 = exp.(-t./T1), exp.(-t./T2)

  M[:,1:2] .*= E2
  M[:,3]   .*= E1
  M[:,3]   .+= 1 .- E1

  eΔθ = exp.(-1im*2π*Δf*t)

  M[:,1:2] .= ((view(M,:,1)+1im*view(M,:,2)).*eΔθ |> x->[real(x) imag(x)])

  return M
end

"""
    freePrec(M, t; Δf, T1, T2)
Same as `freePrec!(M, t; Δf, T1, T2)`, `M` will not be updated.

See also: [`freePrec!`](@ref).
"""
freePrec(M, t; kw...) = freePrec!(copy(M), t; kw...)


text:
"""
# module Potentials

# Examples

```jldoctest
julia>
```
"""
module Potentials

include("InfinitePotentialWells.jl")

end


text:
using Test
import TOMLX: tryparsefile

# These tests need to be updated if the error strings change

tmp = tempname()

@testset "error printing" begin

# Special printing for invalid bare key character
write(tmp, "fooα = 3")
p = tryparsefile(tmp)
err = sprint(Base.showerror, p)
@test contains(err, """
    $tmp:1:4 error: invalid bare key character: 'α'
      fooα = 3
         ^""")

# Error is at EOF
write(tmp, "foo = [1, 2,")
p = tryparsefile(tmp)
err = sprint(Base.showerror, p)
@test contains(err, """
    $tmp:1:12 error: unexpected end of file, expected a value
      foo = [1, 2,
                  ^""")
# A bit of unicode
write(tmp, "\"fαβ\" = [1.2, 1.2.3]")
p = tryparsefile(tmp)
err = sprint(Base.showerror, p)
@test contains(err, """
        $tmp:1:18 error: failed to parse value
          "fαβ" = [1.2, 1.2.3]
                           ^""")

end # testset