text:
module BSeries

@doc read(joinpath(dirname(@__DIR__), "README.md"), String) BSeries

using Reexport: @reexport

@reexport using RootedTrees
using RootedTrees: AbstractRootedTree

@reexport using OrderedCollections: OrderedDict

using Requires: @require

using Latexify: Latexify, LaTeXString

@reexport using Polynomials: Polynomials, Polynomial

export TruncatedBSeries, ExactSolution

export bseries, substitute, compose, evaluate

export modified_equation, modifying_integrator

export elementary_differentials

export MultirateInfinitesimalSplitMethod

# Types used for traits
# These traits may decide between different algorithms based on the
# corresponding complexity etc.
abstract type AbstractEvaluation end
struct LazyEvaluation     <: AbstractEvaluation end # lazy evaluation up to arbitrary order
struct EagerEvaluation    <: AbstractEvaluation end # eager evaluation up to a fixed order
struct MemoizedEvaluation <: AbstractEvaluation end # lazy evaluation memoizing the results

evaluation_type(::Any) = LazyEvaluation() # assume lazy evaluation by default
evaluation_type(::AbstractDict) = EagerEvaluation() # dictionaries store results eagerly

# internal interface
# TODO: `bseries(series::AbstractBSeries) = series` or a `copy`?
# TODO: `bseries(series, order)`` returns a `copy`? up to `order`, if available
# TODO: `truncate!(series, order)`

"""
    TruncatedBSeries

A struct that can describe B-series of both numerical integration methods
(where the coefficient of the empty tree is unity) and right-hand sides of an
ordinary differential equation and perturbations thereof (where the coefficient
of the empty tree is zero) up to a prescribed [`order`](@ref).

Generally, this kind of `struct` should be constructed via [`bseries`](@ref)
or one of the other functions returning a B-series, e.g.,
[`modified_equation`](@ref) or [`modifying_integrator`](@ref).
"""
struct TruncatedBSeries{T<:AbstractRootedTree, V} <: AbstractDict{T, V}
  coef::OrderedDict{T, V}
end

TruncatedBSeries{T, V}() where {T, V} = TruncatedBSeries{T, V}(OrderedDict{T, V}())

# general interface methods of `AbstractDict` for `TruncatedBSeries`
@inline Base.iterate(series::TruncatedBSeries) = iterate(series.coef)
@inline Base.iterate(series::TruncatedBSeries, state) = iterate(series.coef, state)

@inline Base.empty(series::TruncatedBSeries) = TruncatedBSeries(empty(series.coef))
@inline Base.empty(series::TruncatedBSeries, ::Type{T}, ::Type{V}) where {T, V} = TruncatedBSeries(empty(series.coef, T, V))
@inline Base.empty!(series::TruncatedBSeries) = (empty!(series.coef); series)

@inline Base.length(series::TruncatedBSeries) = length(series.coef)

@inline Base.getindex(series::TruncatedBSeries, t::AbstractRootedTree) = getindex(series.coef, t)
@inline Base.setindex!(series::TruncatedBSeries, val, t::AbstractRootedTree) = setindex!(series.coef, val, t)

@inline Base.get(series::TruncatedBSeries, t::AbstractRootedTree, default) = get(series.coef, t, default)
@inline Base.get(f::F, series::TruncatedBSeries, t::AbstractRootedTree) where {F<:Function} = get(f, series.coef, t)

@inline Base.getkey(series::TruncatedBSeries, t::AbstractRootedTree, default) = getkey(series.coef, t, default)

@inline Base.delete!(series::TruncatedBSeries, t::AbstractRootedTree) = (delete!(series.coef, t); series)

@inline Base.pop!(series::TruncatedBSeries, t::AbstractRootedTree) = pop!(series.coef, t)
@inline Base.pop!(series::TruncatedBSeries, t::AbstractRootedTree, default) = pop!(series.coef, t, default)

@inline Base.sizehint!(series::TruncatedBSeries, n) = sizehint!(series.coef, n)

# internal interface of B-series
# evaluation_type(::TruncatedBSeries) = EagerEvaluation() # is already the default for dicts

# TODO: TruncatedBSeries
# This assumes that the `coef` are always stored with increasing `order` of the
# rooted trees and that the underlying data format in OrderedCollections.jl does
# not change. A general fallback would be
#   RootedTrees.order(series::TruncatedBSeries) = maximum(order, keys(series.coef))
# but that is O(n) instead of O(1). Since we do not consider the constructor as
# public API but as internal implementation detail, users violating this assumption
# are outside of the public API and may run into self-made problems.
"""
    order(series::TruncatedBSeries)

The maximal `order` of a rooted tree with non-vanishing coefficient in the
truncated B-series `series`.
"""
RootedTrees.order(series::TruncatedBSeries) = order(series.coef.keys[end])

"""
    ExactSolution{V}()

Lazy representation of the B-series of the exact solution of an ordinary
differential equation using coefficients of type at least as representative as
`V`.
"""
struct ExactSolution{V} end

Base.getindex(::ExactSolution{V}, t::AbstractRootedTree) where {V} = convert(V, 1//1) / γ(t)

# general interface methods of iterators for `ExactSolution`
Base.IteratorSize(::Type{<:ExactSolution}) = Base.SizeUnknown()
Base.eltype(::Type{ExactSolution{V}}) where {V} = V
Base.valtype(exact::ExactSolution) = valtype(typeof(exact))
Base.valtype(::Type{ExactSolution{V}}) where {V} = V

function Base.iterate(exact::ExactSolution)
  iterator = RootedTreeIterator(1)
  t, state = iterate(iterator)
  (exact[t], (state, iterator))
end

function Base.iterate(exact::ExactSolution, state_iterator)
  state, iterator = state_iterator
  t_state = iterate(iterator, state)
  if t_state === nothing
    iterator = RootedTreeIterator(iterator.order + 1)
    t_state = iterate(iterator)
  end
  t, state = t_state
  (exact[t], (state, iterator))
end

# internal interface of B-series
# @inline evaluation_type(::ExactSolution) = LazyEvaluation() # this is the default assumption

"""
    ExactSolution(series_integrator)

A representation of the B-series of the exact solution of an ODE using the same
type of coefficients as the B-series `series_integrator`.
"""
function ExactSolution(series_integrator)
  ExactSolution(series_integrator, evaluation_type(series_integrator))
end

function ExactSolution(series_integrator, ::LazyEvaluation)
  ExactSolution{valtype(series_integrator)}()
end

function ExactSolution(series_integrator, ::EagerEvaluation)
  exact = ExactSolution{valtype(series_integrator)}()
  T = keytype(series_integrator)
  V = valtype(series_integrator)
  series = TruncatedBSeries{T, V}()
  series_keys = keys(series_integrator)

  iter = iterate(series_keys)
  if iter !== nothing
    t, t_state = iter
    if isempty(t)
      series[t] = one(V)
      iter = iterate(series_keys, t_state)
    end
  end

  while iter !== nothing
    t, t_state = iter
    series[t] = exact[t]
    iter = iterate(series_keys, t_state)
  end
  series
end

"""
    bseries(rk::RungeKuttaMethod, order)
    bseries(A::AbstractMatrix, b::AbstractVector, c::AbstractVector, order)

Compute the B-series of the Runge-Kutta method `rk` with Butcher coefficients
`A, b, c` up to a prescribed integer `order`.

!!! note "Normalization by elementary differentials"
    The coefficients of the B-series returned by this method need to be
    multiplied by a power of the time step divided by the `symmetry` of the
    rooted tree and multiplied by the corresponding elementary differential
    of the input vector field ``f``.
    See also [`evaluate`](@ref).
"""
function bseries(rk::RungeKuttaMethod, order)
  V_tmp = eltype(rk)
  if V_tmp <: Integer
    # If people use integer coefficients, they will likely want to have results
    # as exact as possible. However, general terms are not integers. Thus, we
    # use rationals instead.
    V = Rational{V_tmp}
  else
    V = V_tmp
  end
  series = TruncatedBSeries{RootedTree{Int, Vector{Int}}, V}()

  series[rootedtree(Int[])] = one(V)
  for o in 1:order
    for t in RootedTreeIterator(o)
      series[copy(t)] = elementary_weight(t, rk)
    end
  end

  return series
end

function bseries(A::AbstractMatrix, b::AbstractVector, c::AbstractVector,
                 order)
  rk = RungeKuttaMethod(A, b, c)
  bseries(rk, order)
end

# TODO: bseries(rk::RungeKuttaMethod)
# should create a lazy version, optionally a memoized one

"""
    bseries(ark::AdditiveRungeKuttaMethod, order)

Compute the B-series of the additive Runge-Kutta method `ark` up to a prescribed
integer `order`.

!!! note "Normalization by elementary differentials"
    The coefficients of the B-series returned by this method need to be
    multiplied by a power of the time step divided by the `symmetry` of the
    colored rooted tree and multiplied by the corresponding elementary
    differential of the input vector fields ``f^\\nu``.
    See also [`evaluate`](@ref).
"""
function bseries(ark::AdditiveRungeKuttaMethod, order)
  if length(ark.rks) != 2
    throw(ArgumentError("Only AdditiveRungeKuttaMethod with a dual splitting are supported. Got an ARK with $(length(ark.rks)) methods."))
  end

  V_tmp = eltype(ark)
  if V_tmp <: Integer
    # If people use integer coefficients, they will likely want to have results
    # as exact as possible. However, general terms are not integers. Thus, we
    # use rationals instead.
    V = Rational{V_tmp}
  else
    V = V_tmp
  end
  series = TruncatedBSeries{BicoloredRootedTree{Int, Vector{Int}, Vector{Bool}}, V}()

  series[rootedtree(Int[], Bool[])] = one(V)
  for o in 1:order
    for t in BicoloredRootedTreeIterator(o)
      series[copy(t)] = elementary_weight(t, ark)
    end
  end

  return series
end

# TODO: bseries(ark::AdditiveRungeKuttaMethod)
# should create a lazy version, optionally a memoized one

# TODO: Documentation, Base.show, export etc.
"""
    MultirateInfinitesimalSplitMethod(A, D, G, c)

# References

- Knoth, Oswald, and Joerg Wensch.
  "Generalized split-explicit Runge-Kutta methods for the
  compressible Euler equations".
  Monthly Weather Review 142, no. 5 (2014): 2067-2081.
  [DOI: 10.1175/MWR-D-13-00068.1](https://doi.org/10.1175/MWR-D-13-00068.1)

!!! warning "Experimental code"
    This code is considered to be experimental at the moment
    and can change any time.
"""
struct MultirateInfinitesimalSplitMethod{T, PolyMatT<:AbstractMatrix{<:Polynomial{T}}, MatT<:AbstractMatrix{T}, VecT<:AbstractVector{T}} <: RootedTrees.AbstractTimeIntegrationMethod
  A::PolyMatT
  D::MatT
  G::MatT
  c::VecT
end

# TODO: Deduce `c` from other parameters?
function MultirateInfinitesimalSplitMethod(A::AbstractMatrix{<:Polynomial},
                                           D::AbstractMatrix,
                                           G::AbstractMatrix,
                                           c::AbstractVector)
  T = promote_type(eltype(eltype(A)), eltype(D), eltype(G), eltype(c))
  PolyT = typeof(zero(first(A)) + zero(T))
  _A = PolyT.(A)
  _D = T.(D)
  _G = T.(G)
  _c = T.(c)
  return MultirateInfinitesimalSplitMethod(_A, _D, _G, _c)
end

Base.eltype(mis::MultirateInfinitesimalSplitMethod{T}) where {T} = T

"""
    bseries(mis::MultirateInfinitesimalSplitMethod, order)

Compute the B-series of the multirate infinitesimal split method `mis`
up to a prescribed integer `order`.

!!! note "Normalization by elementary differentials"
    The coefficients of the B-series returned by this method need to be
    multiplied by a power of the time step divided by the `symmetry` of the
    colored rooted tree and multiplied by the corresponding elementary
    differential of the input vector fields ``f^\\nu``.
    See also [`evaluate`](@ref).
"""
function bseries(mis::MultirateInfinitesimalSplitMethod, order)
  V_tmp = eltype(mis)
  if V_tmp <: Integer
    # If people use integer coefficients, they will likely want to have results
    # as exact as possible. However, general terms are not integers. Thus, we
    # use rationals instead.
    V = Rational{V_tmp}
  else
    V = V_tmp
  end

  prototype_scalar = TruncatedBSeries{BicoloredRootedTree{Int, Vector{Int}, Vector{Bool}}, V}()
  prototype_scalar[rootedtree(Int[], Bool[])] = one(V)

  poly_one = Polynomial{V, :x}([one(V)])
  poly_zero = zero(poly_one)
  prototype_polynomial = TruncatedBSeries{BicoloredRootedTree{Int, Vector{Int}, Vector{Bool}}, typeof(poly_one)}()
  prototype_polynomial[rootedtree(Int[], Bool[])] = poly_one

  A = mis.A
  D = mis.D
  G = mis.G
  ns = size(A, 2)
  Z = Vector{typeof(prototype_polynomial)}(undef, ns + 1)
  η = Vector{typeof(prototype_scalar)}(undef, ns + 1)
  for i in 1:ns+1
    Z[i] = copy(prototype_polynomial)
    η[i] = copy(prototype_scalar)
  end

  for o in 1:order
    for t_ in BicoloredRootedTreeIterator(o)
      t = copy(t_)
      for i in 1:ns+1
        phi = poly_zero
        for j in 1:i-1
          r = poly_one
          for subtree in RootedTrees.SubtreeIterator(t)
            if subtree.color_sequence[1]
              r = r * Z[i][subtree]
            else
              r = r * η[j][subtree]
            end
          end

          v = [one(V)]
          for k in 0:length(A[i, j])-1
             phi = phi + Polynomial{V, :x}(v) * A[i, j][k] * r
             v = [0; v]
          end
        end

        for j in 1:i-1
          phi = phi + G[i,j] * (η[j][t] - η[1][t])
        end

        phi = Polynomials.integrate(phi)
        phi[0] = 0
        for j in 1:i-1
          phi[0] = phi[0] + D[i,j] * (η[j][t] - η[1][t])
        end

        Z[i][t] = phi
        η[i][t] = phi(1)
      end
    end
  end

  series = η[end]
  return series
end

# TODO: bseries(mis::MultirateInfinitesimalSplitMethod)
# should create a lazy version, optionally a memoized one

"""
    substitute(b, a, t::AbstractRootedTree)

Compute the coefficient correspoding to the tree `t` of the B-series that is
formed by substituting the B-series `b` into the B-series `a`. It is assumed
that the B-series `b` has the coefficient zero of the empty tree.

# References

Section 3.2 of
- Philippe Chartier, Ernst Hairer, Gilles Vilmart (2010)
  Algebraic Structures of B-series.
  Foundations of Computational Mathematics
  [DOI: 10.1007/s10208-010-9065-1](https://doi.org/10.1007/s10208-010-9065-1)
"""
function substitute(b, a, t::AbstractRootedTree)
  result = zero(first(values(a)) * first(values(b)))

  for (forest, skeleton) in PartitionIterator(t)
    update = a[skeleton]
    update isa Rational && iszero(update) && continue
    for tree in forest
      update *= b[tree]
    end
    result += update
  end

  return result
end

"""
    substitute(b, a)

Substitute the B-series `b` into the B-series `a`. It is assumed that the
B-series `b` has the coefficient zero of the empty tree.

In the notation of Chartier, Hairer and Vilmart (2010), we have
`substitute(b, a) = b ★ a`.

# References

Section 3.2 of
- Philippe Chartier, Ernst Hairer, Gilles Vilmart (2010)
  Algebraic Structures of B-series.
  Foundations of Computational Mathematics
  [DOI: 10.1007/s10208-010-9065-1](https://doi.org/10.1007/s10208-010-9065-1)
"""
function substitute(b, a)
  series_keys = keys(b)
  series = empty(b)

  t = first(series_keys)
  @assert isempty(t)
  series[t] = a[t]

  for t in Iterators.drop(series_keys, 1)
    coefficient = substitute(b, a, t)
    series[t] = coefficient
  end

  return series
end

"""
    compose(b, a, t::RootedTree)

Compute the coefficient correspoding to the tree `t` of the B-series that is
formed by composing the B-series `a` with the B-series `b`. It is assumed that
the B-series `b` has the coefficient unity of the empty tree.

# References

Section 3.1 of
- Philippe Chartier, Ernst Hairer, Gilles Vilmart (2010)
  Algebraic Structures of B-series.
  Foundations of Computational Mathematics
  [DOI: 10.1007/s10208-010-9065-1](https://doi.org/10.1007/s10208-010-9065-1)
"""
function compose(b, a, t::RootedTree)
  result = zero(first(values(a)) * first(values(b)))

  for (forest, subtree) in SplittingIterator(t)
    update = a[subtree]
    update isa Rational && iszero(update) && continue
    for tree in forest
      update *= b[tree]
    end
    result += update
  end

  return result
end

"""
    compose(b, a; normalize_stepsize=false)

Compose the B-series `a` with the B-series `b`. It is assumed that the B-series
`b` has the coefficient unity of the empty tree.

In the notation of Chartier, Hairer and Vilmart (2010), we have
`compose(b, a) = b ⋅ a`.

If `normalize_stepsize = true`, the coefficients of the returned B-series will
are divied by `2^order(t)` for each rooted tree `t`. This normalizes the step
size so that the resulting numerical integrator B-series uses the same step size
as the input series (instead of a doubled step size).

# References

Section 3.1 of
- Philippe Chartier, Ernst Hairer, Gilles Vilmart (2010)
  Algebraic Structures of B-series.
  Foundations of Computational Mathematics
  [DOI: 10.1007/s10208-010-9065-1](https://doi.org/10.1007/s10208-010-9065-1)
"""
function compose(b, a; normalize_stepsize=false)
  series_keys = keys(b)
  series = empty(b)

  for t in series_keys
    coefficient = compose(b, a, t)
    if normalize_stepsize
      coefficient /= 2^order(t)
    end
    series[t] = coefficient
  end

  return series
end

"""
    evaluate(f, u, dt, series, reduce_order_by=0)

Evaluate the B-series `series` specialized to the ordinary differential equation
``u'(t) = f(u(t))`` with vector field `f` and dependent variables `u` for a
time step size `dt`.

Here, `u` is assumed to be a vector of symbolic variables and `f` is assumed
to be a vector of expressions in these variables for plain B-series. For
B-series with colored trees, `f` must be a tuple of vectors of expressions in
the variables `u`. Currently, symbolic variables from

- [SymEngine.jl](https://github.com/symengine/SymEngine.jl),
- [SymPy.jl](https://github.com/JuliaPy/SymPy.jl), and
- [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl)

are supported.

The powers of `dt` can be controlled by `reduce_order_by` to make them different
from the usual `order(t)` for a rooted tree `t`. This can be useful in the
context of [`modified_equation`](@ref)s or [`modifying_integrator`](@ref)s,
where the B-series coeffcients are those of ``h fₕ``, i.e., they contain an
additional power of `dt`. In this case, the B-series of the vector field can
be obtained using `reduce_order_by = 1`.

# References

Section 3.2 of
- Philippe Chartier, Ernst Hairer, Gilles Vilmart (2010)
  Algebraic Structures of B-series.
  Foundations of Computational Mathematics
  [DOI: 10.1007/s10208-010-9065-1](https://doi.org/10.1007/s10208-010-9065-1)
"""
function evaluate(f, u, dt, series, reduce_order_by=0)
  evaluate(f, u, dt, series, evaluation_type(series), reduce_order_by)
end

function evaluate(f, u, dt, series, ::EagerEvaluation, reduce_order_by)
  differentials = elementary_differentials(f, u, order(series))

  # An additive decomposition is indicated by a tuple of vectors. A single
  # vector field is assumed to be a vector, as everywhere else.
  if f isa NTuple{N, AbstractVector} where {N}
    result = zero(first(f))
  else
    result = zero(f)
  end

  for t in keys(series)
    # Otherwise, SymPy.jl might result in
    #   DomainError with -1:
    #   Cannot raise an integer x to a negative power -1.
    #   Convert input to float.
    power = order(t) - reduce_order_by
    if power > 0
      dt_term = dt^power
    elseif power == 0
      dt_term = one(dt)
    else
      dt_term = inv(dt^(-power))
    end

    result += dt_term / symmetry(t) * series[t] * differentials[t]
  end
  result
end

"""
    modified_equation(series_integrator)

Compute the B-series of the modified equation of the time integration method
with B-series `series_integrator`.

Given an ordinary differential equation (ODE) ``u'(t) = f(u(t))`` and a
Runge-Kutta method, the idea is to interpret the numerical solution with
given time step size as exact solution of a modified ODE ``u'(t) = fₕ(u(t))``.
This method returns the B-series of ``h fₕ``.

!!! note "Normalization by elementary differentials"
    The coefficients of the B-series returned by this method need to be
    multiplied by a power of the time step divided by the `symmetry` of the
    rooted tree and multiplied by the corresponding elementary differential
    of the input vector field ``f``.
    See also [`evaluate`](@ref).

# References

Section 3.2 of
- Philippe Chartier, Ernst Hairer, Gilles Vilmart (2010)
  Algebraic Structures of B-series.
  Foundations of Computational Mathematics
  [DOI: 10.1007/s10208-010-9065-1](https://doi.org/10.1007/s10208-010-9065-1)
"""
function modified_equation(series_integrator)
  modified_equation(series_integrator, evaluation_type(series_integrator))
end

function modified_equation(series_integrator, ::EagerEvaluation)
  V = valtype(series_integrator)

  # B-series of the exact solution
  # We could just use the lazy version
  #   series_ex = ExactSolution{V}()
  # However, we need to access elements of `series_ex` more than once in the
  # subsitution below. Thus, it's cheaper to compute every entry only once and
  # re-use it later.
  series_ex = ExactSolution(series_integrator)

  # Prepare B-series of the modified equation
  series_keys = keys(series_integrator)
  series = empty(series_integrator)
  for t in series_keys
    series[t] = zero(V)
  end

  iter = iterate(series_keys)
  if iter !== nothing
    t, t_state = iter
    if isempty(t)
      iter = iterate(series_keys, t_state)
      if iter !== nothing
        t, t_state = iter
      end
    end

    series[t] = series_integrator[t]
    iter = iterate(series_keys, t_state)
  end

  # Recursively solve
  #   substitute(series, series_ex, t) == series_integrator[t]
  # This works because
  #   substitute(series, series_ex, t) = series[t] + lower order terms
  # Since the `keys` are ordered, we don't need to use nested loops of the form
  #   for o in 2:order
  #     for _t in RootedTreeIterator(o)
  #       t = copy(_t)
  # which are slightly less efficient due to additional computations and
  # allocations.
  while iter !== nothing
    t, t_state = iter
    series[t] += series_integrator[t] - substitute(series, series_ex, t)
    iter = iterate(series_keys, t_state)
  end

  return series
end

"""
    modified_equation(rk::RungeKuttaMethod, order)
    modified_equation(A::AbstractMatrix, b::AbstractVector, c::AbstractVector,
                      order)

Compute the B-series of the [`modified_equation`](@ref) of the Runge-Kutta
method `rk` with Butcher coefficients `A, b, c` up to the prescribed `order`.

!!! note "Normalization by elementary differentials"
    The coefficients of the B-series returned by this method need to be
    multiplied by a power of the time step divided by the `symmetry` of the
    rooted tree and multiplied by the corresponding elementary differential
    of the input vector field ``f``.
    See also [`evaluate`](@ref).
"""
function modified_equation(rk::RungeKuttaMethod, order)
  # B-series of the Runge-Kutta method
  series = bseries(rk, order)
  modified_equation(series)
end

function modified_equation(A::AbstractMatrix, b::AbstractVector,
                           c::AbstractVector, order)
  rk = RungeKuttaMethod(A, b, c)
  modified_equation(rk, order)
end

"""
    modified_equation(f, u, dt, series_integrator)

Compute the B-series of the [`modified_equation`](@ref) of the time integration
method with B-series `series_integrator` with respect to the ordinary
differential equation ``u'(t) = f(u(t))`` with vector field `f` and dependent
variables `u` for a time step size `dt`.

Here, `u` is assumed to be a vector of symbolic variables and `f` is assumed
to be a vector of expressions in these variables for plain B-series. For
B-series with colored trees, `f` must be a tuple of vectors of expressions in
the variables `u`. Currently, symbolic variables from

- [SymEngine.jl](https://github.com/symengine/SymEngine.jl),
- [SymPy.jl](https://github.com/JuliaPy/SymPy.jl), and
- [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl)

are supported.
"""
function modified_equation(f, u, dt, series_integrator)
  series = modified_equation(series_integrator)
  reduce_order_by = 1 # reduce the powers of `dt` by one since dt*f is given
                      # by recursively solving `substitute`
  evaluate(f, u, dt, series, reduce_order_by)
end

"""
    modified_equation(f, u, dt, rk::RungeKuttaMethod, order)
    modified_equation(f, u, dt,
                      A::AbstractMatrix, b::AbstractVector, c::AbstractVector,
                      order)

Compute the B-series of the [`modified_equation`](@ref) of the Runge-Kutta
method `rk` with Butcher coefficients `A, b, c` up to the prescribed `order` with
respect to the ordinary differential equation ``u'(t) = f(u(t))`` with vector
field `f` and dependent variables `u` for a time step size `dt`.

Here, `u` is assumed to be a vector of symbolic variables and `f` is assumed
to be a vector of expressions in these variables for plain B-series. For
B-series with colored trees, `f` must be a tuple of vectors of expressions in
the variables `u`. Currently, symbolic variables from

- [SymEngine.jl](https://github.com/symengine/SymEngine.jl),
- [SymPy.jl](https://github.com/JuliaPy/SymPy.jl), and
- [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl)

are supported.
"""
function modified_equation(f, u, dt, rk::RungeKuttaMethod, order)
  series_integrator = bseries(rk, order)
  modified_equation(f, u, dt, series_integrator)
end

function modified_equation(f, u, dt,
                           A::AbstractMatrix, b::AbstractVector,
                           c::AbstractVector, order)
  rk = RungeKuttaMethod(A, b, c)
  modified_equation(f, u, dt, rk, order)
end

"""
    modifying_integrator(series_integrator)

Compute the B-series of a "modifying integrator" equation of the time
integration method with B-series `series_integrator`.

Given an ordinary differential equation (ODE) ``u'(t) = f(u(t))`` and a
Runge-Kutta method, the idea is to find a modified ODE ``u'(t) = fₕ(u(t))``
such that the numerical solution with given time step size is the exact solution
of the original ODE. This method returns the B-series of ``h fₕ``.

!!! note "Normalization by elementary differentials"
    The coefficients of the B-series returned by this method need to be
    multiplied by a power of the time step divided by the `symmetry` of the
    rooted tree and multiplied by the corresponding elementary differential
    of the input vector field ``f``.
    See also [`evaluate`](@ref).

# References

Section 3.2 of
- Philippe Chartier, Ernst Hairer, Gilles Vilmart (2010)
  Algebraic Structures of B-series.
  Foundations of Computational Mathematics
  [DOI: 10.1007/s10208-010-9065-1](https://doi.org/10.1007/s10208-010-9065-1)
"""
function modifying_integrator(series_integrator)
  modifying_integrator(series_integrator, evaluation_type(series_integrator))
end

function modifying_integrator(series_integrator, ::EagerEvaluation)
  V = valtype(series_integrator)

  # B-series of the exact solution
  # Since we access each coefficient of this B-series only once below, it's
  # okay to use the lazy version of the exact solution.
  series_ex = ExactSolution{V}()

  # Prepare B-series of the modifying integrator equation
  series_keys = keys(series_integrator)
  series = empty(series_integrator)
  for t in series_keys
    series[t] = zero(V)
  end

  iter = iterate(series_keys)
  if iter !== nothing
    t, t_state = iter
    if isempty(t)
      iter = iterate(series_keys, t_state)
      if iter !== nothing
        t, t_state = iter
      end
    end

    series[t] = series_integrator[t]
    iter = iterate(series_keys, t_state)
  end

  # Recursively solve
  #   substitute(series, series_integrator, t) == series_ex[t]
  # This works because
  #   substitute(series, series_integrator, t) = series[t] + lower order terms
  # Since the `keys` are ordered, we don't need to use nested loops of the form
  #   for o in 2:order
  #     for _t in RootedTreeIterator(o)
  #       t = copy(_t)
  # which are slightly less efficient due to additional computations and
  # allocations.
  while iter !== nothing
    t, t_state = iter
    series[t] += series_ex[t] - substitute(series, series_integrator, t)
    iter = iterate(series_keys, t_state)
  end

  return series
end

"""
    modifying_integrator(rk::RungeKuttaMethod, order)
    modifying_integrator(A::AbstractMatrix, b::AbstractVector, c::AbstractVector,
                         order)

Compute the B-series of the [`modifying_integrator`](@ref) equation of the
Runge-Kutta method with Butcher coefficients `A, b, c` up to the prescribed
`order`.

!!! note "Normalization by elementary differentials"
    The coefficients of the B-series returned by this method need to be
    multiplied by a power of the time step divided by the `symmetry` of the
    rooted tree and multiplied by the corresponding elementary differential
    of the input vector field ``f``.
    See also [`evaluate`](@ref).
"""
function modifying_integrator(rk::RungeKuttaMethod, order)
  # B-series of the Runge-Kutta method
  series = bseries(rk, order)
  modifying_integrator(series)
end

function modifying_integrator(A::AbstractMatrix, b::AbstractVector,
                              c::AbstractVector, order)
  rk = RungeKuttaMethod(A, b, c)
  modifying_integrator(rk, order)
end

"""
    modifying_integrator(f, u, dt, series_integrator)

Compute the B-series of the [`modifying_integrator`](@ref) equation of the
time integration method with B-series `series_integrator` with respect to the
ordinary differential equation ``u'(t) = f(u(t))`` with vector field `f` and
dependent variables `u` for a time step size `dt`.

Here, `u` is assumed to be a vector of symbolic variables and `f` is assumed
to be a vector of expressions in these variables for plain B-series. For
B-series with colored trees, `f` must be a tuple of vectors of expressions in
the variables `u`. Currently, symbolic variables from

- [SymEngine.jl](https://github.com/symengine/SymEngine.jl),
- [SymPy.jl](https://github.com/JuliaPy/SymPy.jl), and
- [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl)

are supported.
"""
function modifying_integrator(f, u, dt, series_integrator)
  series = modifying_integrator(series_integrator)
  reduce_order_by = 1 # reduce the powers of `dt` by one since dt*f is given
                      # by recursively solving `substitute`
  evaluate(f, u, dt, series, reduce_order_by)
end

"""
    modifying_integrator(f, u, dt, rk::RungeKuttaMethod, order)
    modifying_integrator(f, u, dt,
                         A::AbstractMatrix, b::AbstractVector, c::AbstractVector,
                         order)

Compute the B-series of the [`modifying_integrator`](@ref) equation of the
Runge-Kutta method with Butcher coefficients `A, b, c` up to the prescribed
`order` with respect to the ordinary differential equation ``u'(t) = f(u(t))``
with vector field `f` and dependent variables `u` for a time step size `dt`.

Here, `u` is assumed to be a vector of symbolic variables and `f` is assumed
to be a vector of expressions in these variables for plain B-series. For
B-series with colored trees, `f` must be a tuple of vectors of expressions in
the variables `u`. Currently, symbolic variables from

- [SymEngine.jl](https://github.com/symengine/SymEngine.jl),
- [SymPy.jl](https://github.com/JuliaPy/SymPy.jl), and
- [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl)

are supported.
"""
function modifying_integrator(f, u, dt, rk::RungeKuttaMethod, order)
  series_integrator = bseries(rk, order)
  modifying_integrator(f, u, dt, series_integrator)
end

function modifying_integrator(f, u, dt,
                              A::AbstractMatrix, b::AbstractVector,
                              c::AbstractVector, order)
  rk = RungeKuttaMethod(A, b, c)
  modifying_integrator(f, u, dt, rk, order)
end

"""
    elementary_differentials(f::AbstractVector, u, order)

Compute all elementary differentials of the vector field `f` with independent
variables `u` up to the given `order`. The return value can be indexed by
rooted trees to obtain the corresponding elementary differential.
"""
function elementary_differentials(f::AbstractVector, u, order)
  order >= 1 || throw(ArgumentError("The `order` must be at least one (got $order)"))
  differentials = OrderedDict{RootedTree{Int, Vector{Int}}, typeof(f)}()

  t = rootedtree(Int[])
  differentials[t] = one.(f)

  t = rootedtree([1])
  differentials[t] = f

  # Compute all necessary partial derivatives at first
  derivatives = Array{eltype(f)}[]
  push!(derivatives, f)
  for o in 1:(order-1)
    d = similar(f, eltype(f), (length(f), ntuple(_ -> length(u), o)...))
    _compute_partial_derivatives!(d, f, u, derivatives[o])
    push!(derivatives, d)
  end

  # Compute all elementary differentials
  for o in 2:order
    for _t in RootedTreeIterator(o)
      t = copy(_t)
      differentials[t] = elementary_differential(f, t, differentials, derivatives)
    end
  end

  return differentials
end

function _compute_partial_derivatives!(d, f, u, lower_derivatives)
  for idx in CartesianIndices(d)
    idx_tuple = Tuple(idx)
    u_idx = Base.tail(idx_tuple)

    # All this B-series analysis only really makes sense for smooth functions.
    # Hence, we can use the symmetry of the partial derivatives to speed-up
    # the computations - Hermann Amandus Schwarz helps us again!
    issorted(u_idx) || continue

    # Next, we re-use already computed partial derivatives. Thus, we do not
    # compute the full set of derivatives as in
    #   f_idx = first(idx_tuple)
    #   partial_derivative = f[f_idx]
    #   for i in u_idx
    #     partial_derivative = compute_derivative(partial_derivative, u[i])
    #   end
    # but we re-use the already computed `lower_derivatives`.
    idx_known = Base.front(idx_tuple)
    idx_new = last(idx_tuple)
    partial_derivative = compute_derivative(lower_derivatives[idx_known...], u[idx_new])

    d[idx] = partial_derivative
  end
end

"""
    compute_derivative(expression, variable)

Internal function specialized on symbolic variables and expressions from

- [SymEngine.jl](https://github.com/symengine/SymEngine.jl),
- [SymPy.jl](https://github.com/JuliaPy/SymPy.jl), and
- [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl)

if these packages are loaded (via Requires.jl).
"""
function compute_derivative end

function __init__()
  @require SymEngine="123dc426-2d89-5057-bbad-38513e3affd8" begin
    using .SymEngine: SymEngine

    function compute_derivative(expression::SymEngine.Basic, variable::SymEngine.Basic)
      SymEngine.diff(expression, variable)
    end

    latexify_default_dt(::Type{SymEngine.Basic}) = SymEngine.symbols("h")
  end

  @require SymPy="24249f21-da20-56a4-8eb1-6a02cf4ae2e6" begin
    using .SymPy: SymPy

    function compute_derivative(expression::SymPy.Sym, variable::SymPy.Sym)
      SymPy.diff(expression, variable)
    end

    latexify_default_dt(::Type{SymPy.Sym}) = SymPy.symbols("h", real=true)
  end

  @require Symbolics="0c5d862f-8b57-4792-8d23-62f2024744c7" begin
    using .Symbolics: Symbolics

    function compute_derivative(expression::Symbolics.Num, variable::Symbolics.Num)
      Symbolics.expand_derivatives(Symbolics.Differential(variable)(expression))
    end

    # nested macro calls do not work so we need to define this part in another file
    include("latexify_symbolics.jl")
  end
end

function elementary_differential(f, t, differentials, derivatives)
  subtr = RootedTrees.subtrees(t)
  result = similar(f)

  input_differentials = ntuple(n -> differentials[subtr[n]], length(subtr))
  input_derivative = derivatives[length(subtr) + 1]

  _compute_elementary_differential!(result, input_differentials, input_derivative)

  return result
end

@generated function _compute_elementary_differential!(result, input_differentials, input_derivative::AbstractArray{T,N}) where {T,N}
  quote
    for i in eachindex(result)
      val = zero(eltype(result))

      Base.Cartesian.@nloops $(N-1) j input_derivative begin
        # A naive version uses the full input of all possible combinations of
        # partial derivatives. If these were all computed in
        # `_compute_partial_derivatives!` above, we could just use
        #   term = Base.Cartesian.@ncall $(N-1) getindex input_derivative i j
        # and get rid of `sorted_j` here. However, all this B-series analysis
        # only really makes sense for smooth functions. Hence, we can use the
        # symmetry of the partial derivatives to speed-up the computations
        # - Hermann Amandus Schwarz helps us again!
        # TODO: Shall we use the non-allocating version `TupleTools.sort∘tuple`
        # instead of `sort!∘Base.vect`?
        sorted_j = Base.Cartesian.@ncall $(N-1) sort!∘Base.vect j
        term = Base.Cartesian.@ncall $(N-1) getindex input_derivative i d -> sorted_j[d]
        Base.Cartesian.@nexprs $(N-1) d -> term *= input_differentials[d][j_d]
        val += term
      end

      result[i] = val
    end
  end
end

"""
    elementary_differentials(fs::NTuple{2, AbstractVector}, u, order)

Compute all elementary differentials of the sum of the two vector fields `f`
with independent variables `u` up to the given `order`. The return value
can be indexed by (bi-) colored rooted trees to obtain the corresponding
elementary differential.
"""
function elementary_differentials(f::NTuple{2, AbstractVector}, u, order)
  order >= 1 || throw(ArgumentError("The `order` must be at least one (got $order)"))
  N = 2 # length(fs)

  # Empty bicolored tree
  t = rootedtree(Int[], Bool[])
  differentials = OrderedDict{typeof(t), promote_type(map(typeof, f)...)}()
  differentials[t] = one.(f[1])

  # Bicolored trees with a single node
  t = rootedtree([1], Bool[0])
  differentials[t] = f[1]

  t = rootedtree([1], Bool[1])
  differentials[t] = f[2]

  # Compute all necessary partial derivatives at first
  derivatives = ntuple(n -> Array{eltype(f[n])}[], N)
  for n in 1:N
    push!(derivatives[n], f[n])
  end
  for o in 1:(order-1)
    for n in 1:N
      d = similar(f[n], eltype(f[n]), (length(f[n]), ntuple(_ -> length(u), o)...))
      _compute_partial_derivatives!(d, f[n], u, derivatives[n][o])
      push!(derivatives[n], d)
    end
  end

  # Compute all elementary differentials
  for o in 2:order
    for _t in BicoloredRootedTreeIterator(o)
      t = copy(_t)
      differentials[t] = elementary_differential(f, t, differentials, derivatives)
    end
  end

  return differentials
end

function elementary_differential(f::NTuple{2, AbstractVector},
                                 t::ColoredRootedTree,
                                 differentials, derivatives)
  subtr = RootedTrees.subtrees(t)
  n = RootedTrees.color_to_index(root_color(t))
  result = similar(f[n])

  input_differentials = ntuple(i -> differentials[subtr[i]], length(subtr))
  input_derivative = derivatives[n][length(subtr) + 1]

  _compute_elementary_differential!(result, input_differentials, input_derivative)

  return result
end

include("latexify.jl")

end # module


text:
using ArgParse
s = ArgParseSettings()
@add_arg_table s begin
    "data_bundles"
        help = "Location of the ChromNet data bundles (i.e. ENCODE_buildX.ChromNet.jld)"
        arg_type = ASCIIString
        nargs = '+'
        required = true
    "--output", "-o"
        help = "Where to save the network JSON file (default is STDOUT)"
        default = "-"
    "--save-cor"
        help = "Save a copy of the raw data correlation matrix"
        metavar = "COR_FILE"
    "--threshold"
        help = "At what absolute value threshold do we keep links"
        arg_type = Float64
        default = 0.3
    "--group-score-threshold"
        help = "Below what group score do we discard connecting links [0.0-1.0]"
        arg_type = Float64
        default = 0.8
        metavar = "THRESHOLD"
    "--quiet", "-q"
        help = "Don't print anything"
        action = :store_true
end
args = parse_args(s)

using ChromNet
using JSON
using JLD

# open the data bundle and load the metadata
dataBundles = [jldopen(x) for x in args["data_bundles"]]
metadata = Dict()
lastAssembly = nothing
for bundle in dataBundles

    # ensure all the assemblies match
    a = JLD.exists(bundle, "assembly") ? read(bundle, "assembly") : "GRCh38_alt"
    if lastAssembly != nothing
        @assert lastAssembly == a "Assemby mismatch! ("*lastAssembly*", $a)"
    end
    lastAssembly = a

    metadata = merge(metadata, read(bundle, "metadata"))
end
ids = vcat([read(bundle, "ids") for bundle in dataBundles]...)

# compute the joint correlation matrix
C = streaming_cor(dataBundles, quiet=args["quiet"])
if typeof(args["save-cor"]) != Void
    f = open(args["save-cor"], "w")
    println(f, join(ids, '\t'))
    writedlm(f, C)
    close(f)
end
map(close, dataBundles)
cov2cor!(C)

# compute groups using hclust
groups = build_groups(C, ids)

# compute groupgm matrix and create network
# we add a small amount of regularization to avoid potential singularity issues in user-provided data
G,idsG = build_groupgm(inv(C + 1e-8I), ids, groups)
network = build_network(G, groups, metadata,
    threshold=args["threshold"],
    groupScoreThreshold=args["group-score-threshold"],
    quiet=args["quiet"]
)

# send the network in JSON format to the output stream
outStream = STDOUT
if args["output"] != "-"
    outStream = open(args["output"], "w")
end
println(outStream, json(network))

ext: jl


text:
# creates the new directory
isdir("Picks") ? rm("Picks", recursive=true) : nothing ; mkdir("Picks")

# grabs all images(jpg format) in the current dir.
imgs =[i for i = readdir(".") if split(i,".")[end] == "jpg" ]

# random picks 30 images
picks = Set{String}()
while length(picks) < 30
	push!(picks, rand(imgs))
end

#
foreach(x -> cp(x, "Picks/$x"), picks)


text:
## loss
abstract type LossFunction <: Function end

struct Loss{F <: Array{T} where {T <: LossFunction}}
    functions::F
end

## Diff
abstract type SemDiff end

abstract type DiffFunction <: Function end

struct SemAnalyticDiff{F <: Tuple, F2 <: Union{Tuple, Nothing}} <: SemDiff
    algorithm
    options
    functions::F
    hessian_functions::F2
end

## Obs
abstract type SemObs end

## Imply
abstract type Imply end
abstract type ImplyDefinition <: Imply end

## Interaction
abstract type SemInteraction end

## SEModel
abstract type AbstractSem end

struct Sem{O <: SemObs, I <: Imply, L <: Loss, D <: SemDiff} <: AbstractSem
    observed::O
    imply::I
    loss::L
    diff::D
end

struct CollectionSem{
    O <: Vector{O} where O <: SemObs,
    I <: Vector{I} where I <: Imply,
    L <: Vector{L} where L <: Loss,
    D <: Vector{D} where D <: SemDiff} <: AbstractSem
    observed_vec::O
    imply_vec::I
    loss_vec::L
    diff_vec::D
end

struct MGSem{V <: Vector{AS} where {
    AS <: AbstractSem}, D <: Vector{Vector{Int64}}} <: AbstractSem
    sem_vec::V
    par_subsets::D
end

struct FIMLSem{} <: AbstractSem
end


text:
export segyhdr
# ============================================================================
# Utility functions not for export
function ibmfloat(x::UInt32)
  local fra = ntoh(x)
  local sgn = UInt8((fra & 0x80000000)>>31)
  fra <<= 1
  local exp = Int16(fra >> 25) - Int16(64)
  fra <<= 7
  y = (sgn == 0x00 ? 1.0 : -1.0) * 16.0^exp * signed(fra >> 8)/16777216
  return y
end
#=
  This is actual IBM hexadecimal float, and correctly parses the examples from
  https://en.wikipedia.org/wiki/IBM_hexadecimal_floating_point ; see tests.

  The version in JuliaSeis is wrong because IBM float has a range too wide for
  IEEE single float in Julia.

  The description in the SEG Y manual does something strange with the last bit.

  My last line is computationally expensive; working out the radix shift and
  using >> would be much better.
=#

function trid(i::Int16, fs::Float64, fc::Float64)
  S = ["DH", "HZ", "H1", "H2", "HZ", "HT", "HR"]
  return string(getbandcode(fs, fc=fc), S[i])
end

function mk_lc(lv::Int32)
  L1 = UInt8(div(lv, Int32(36)))
  L2 = UInt8(rem(lv, Int32(36)))
  c1 = 0x30 + L1 + (L1 > 0x09 ? 0x07 : 0x00)
  c2 = 0x30 + L2 + (L2 > 0x09 ? 0x07 : 0x00)
  return String([c1, c2])
end

function auto_coords(xy::Array{Int32, 1}, sc::Array{Int16, 1})
  xy == Int32[0,0] && return (0.0, 0.0)
  lon = Float64(xy[1])
  lat = Float64(xy[2])
  coord_scale = sc[1]
  coord_units = sc[2]
  (coord_scale < 0) && (coord_scale = -1 / coord_scale)
  lat *= coord_scale
  lon *= coord_scale

  #=
  Conversion formula for coord_units == 1 "borrowed" from PASSSOFT segy2sac.
  Assumes (possibly wrong) x-y coordinate origin & spherical Earth.
  coord_units == 3 and coord_units == 4 aren't supported; never seen them.
  =#
  if coord_units == 1
      iflg = lon < 0 ? -1 : 1
      c = 111194.6976
      D = sqrt(lat^2 + lon^2) / c
      lat /= c
      d = cosd(D)
      lon = iflg*acosd(d / cosd(lat))
  else
    lat /= 3600.0
    lon /= 3600.0
  end
  return lat, lon
end

function do_trace(f::IO,
                  passcal::Bool,
                  full::Bool,
                  ll::UInt8,
                  swap::Bool,
                  fh::Array{Int16,1}
                  )

  # First part of trace header is quite standard
  buf = BUF.buf
  ints = BUF.int32_buf
  shorts = BUF.int16_buf
  lat = 0.0
  lon = 0.0
  checkbuf_8!(buf, max(180, length(buf)))
  fast_readbytes!(f, buf, 180)

  intbuf = reinterpret(Int32, buf)
  copyto!(ints, 1, intbuf, 1, 7)
  copyto!(ints, 8, intbuf, 10, 8)
  copyto!(ints, 16, intbuf, 19, 4)

  shortbuf = reinterpret(Int16, buf)
  copyto!(shorts, 1, shortbuf, passcal ? 1 : 15, 4) # shorts[1:4]
  copyto!(shorts, 5, shortbuf, 35, 2)               # shorts[5:6]
  copyto!(shorts, 7, shortbuf, 45, 46)              # shorts[7:52]

  if passcal
    fast_readbytes!(f, buf, 40)

    scale_fac     = fastread(f, Float32)
    inst_no       = fastread(f, UInt16)
    shorts[62]    = fastread(f, Int16)
    ints[21]      = fastread(f, Int32)
    ints[22]      = fastread(f, Int32)
    ints[23]      = fastread(f, Int32)

    setindex!(shorts, shortbuf[10], 53)
    copyto!(shorts, 54, shortbuf, 13, 8)
    setindex!(ints, intbuf[6], 20)

    if swap
      shorts     .= ntoh.(shorts)
      ints       .= ntoh.(ints)
      scale_fac   = bswap(scale_fac)
      inst_no     = bswap(inst_no)
    end

    chars = buf[1:18]
    dt    = getindex(ints,20)
    n     = getindex(shorts, 20)
    fmt   = getindex(shorts, 54)
    nx    = (n == typemax(Int16) ? getindex(ints,21) : Int32(n))
    T     = (fmt == one(Int16) ? Int32 : Int16)
    nb    = checkbuf!(buf, nx, T)

    fast_readbytes!(f, buf, nb)

    # trace processing
    y = reinterpret(T, buf)
    x = Array{Float32,1}(undef, nx)
    if swap
      for i = 1:nx
        x[i] = bswap(y[i])
      end
    else
      copyto!(x, 1, y, 1, nx)
      # faster than reprocessing with fillx_ for long files
    end

    z = getindex(shorts, 5)
    δ = getindex(shorts, 21)
    fs        = sμ / Float64(δ == one(Int16) ? dt : δ)
    gain      = 1.0 / (Float64(scale_fac)  * 10.0^(-1.0*Float64(shorts[24]) / 10.0) / Float64(shorts[23]))
    lat, lon  = auto_coords(ints[18:19], shorts[6:7])
    if abs(lat) < 0.1 && abs(lon) < 0.1
      lat *= 3600.0
      lon *= 3600.0
    end
    el        = Float64(ints[9]) * Float64(abs(z))^(z < zero(Int16) ? -1.0 : 1.0)

    # Create ID
    id_arr = zeros(UInt8,8)
    i = one(Int8)
    o = one(Int8)
    while i < Int8(18)
      if chars[i] == 0x00
        chars[i] = 0x20
      end
      i = i+o
    end
    fill_id!(id_arr, chars, one(Int16), Int16(6), Int16(2), Int16(6))
    id_arr[1] = 0x2e
    id_arr[8] = 0x2e
    deleteat!(id_arr, id_arr.==0x00)

    # Channel string is tedious; use of one-char channel names like "Z"
    inds = Int64[]
    for i in 15:18
      if chars[i] > 0x20
        push!(inds, i)
      end
    end
    ch_arr = isempty(inds) ? UInt8[] : chars[inds]
    nc = length(ch_arr)
    if nc == 1
      c = uppercase(Char(ch_arr[1]))
      if c in ('Z','N','E')
        cha = string(getbandcode(fs), 'H', c)
      else
        cha = "YYY"
      end
    elseif nc == 0
      cha = "YYY"
    else
      cha = String(ch_arr[1:min(3, nc)])
    end
    id = String(id_arr) * cha

  else
    (dt, nx, fmt) = fh
    fast_readbytes!(f, buf, 38)
    ints[26]      = fastread(f, Int32)
    shorts[60]    = fastread(f, Int16)
    ints[27]      = fastread(f, Int32)
    shorts[61]    = fastread(f, Int16)
    shorts[62]    = fastread(f, Int16)
    ints[28]      = fastread(f, Int32)
    ints[29]      = fastread(f, Int32)

    copyto!(ints, 20, intbuf, 1, 5)
    copyto!(shorts, 53, shortbuf, 11, 2)          # shorts[53:54]
    ints[25] = getindex(intbuf, 7)
    copyto!(shorts, 55, shortbuf, 15, 5)          # shorts[55:59]

    T = getindex(segy_ftypes, fmt)
    if (T == Any)
      close(f)
      error(string("Trace data code = ", fmt, "  unsupported!"))
    end
    nb = checkbuf!(buf, nx, T)
    fast_readbytes!(f, buf, nb)

    # trace processing
    x = Array{Float32,1}(undef, nx)
    if T == Int16
      fillx_i16_be!(x, buf, nx, 0)
    elseif T == Int32
      fillx_i32_be!(x, buf, nx, 0)
    elseif T == Int8
      fillx_i8!(x, buf, nx, 0)
    elseif T == Float32
      x .= bswap.(reinterpret(Float32, buf))[1:nx]
    elseif T == UInt32
      y = Array{UInt32,1}(undef, nx)
      fillx_u32_le!(y, buf, nx, 0) # _le because ibmfloat bswaps
      x = ibmfloat.(y)
    end

    if swap
    shorts  .= ntoh.(shorts)
    ints    .= ntoh.(ints)
    end
    z       = getindex(shorts, 5)
    δ       = getindex(shorts, 21)
    fs      = sμ / Float64(δ)

    # not sure about meaning of "dB" in gain constants
    gain    = Float64(ints[25]) * 10.0^(shorts[55] + (shorts[23] + shorts[24])/10.0)
    el      = Float64(ints[9]) * Float64(abs(z))^(z<Int16(0) ? -1.0 : 1.0)

    # Create ID
    sta = string(reinterpret(UInt16, shorts[57]))
    lc = ll > 0x00 ? mk_lc(ints[ll]) : ""
    cha = Int16(10) < shorts[1] < Int16(18) ? trid(shorts[1]-Int16(10), fs, 1.0) : "YYY"
    id = string(".", sta[1:min(lastindex(sta),5)], ".", lc, ".", cha)

  end
  if full == true
    shorts_k = String["trace_id_code", "n_summed_z", "n_summed_h", "data_use",
      "z_sc", "h_sc", "h_units_code", "v_weather",
      "v_subweather", "t_src_uphole", "t_rec_uphole", "src_static_cor",
      "rec_static_cor", "total_static", "t_lag_a", "t_lag_b",
      "t_delay", "t_mute_st", "t_mute_en", "nx",
      "delta", "gain_type", "gain_const", "init_gain",
      "correlated", "sweep_st", "sweep_en", "sweep_len",
      "sweep_type", "sweep_tap_st", "sweep_tap_en", "tap_type",
      "f_alias", "slope_alias", "f_notch", "slope_notch",
      "f_low_cut", "f_high_cut", "slope_low_cut", "slope_high_cut",
      "year", "day", "hour", "minute",
      "second", "time_code", "trace_wt_fac", "geophone_roll_p1",
      "geophone_first_tr", "geophone_last_tr", "gap_size", "overtravel"]
    append!(shorts_k, passcal ?
        String["total_static_hi", "data_form", "ms", "trigyear",
        "trigday", "trighour", "trigminute", "trigsecond",
        "trigms", "not_to_be_used"] :
        String["shot_scalar", "trace_units_code", "trans_exp", "trans_units_code",
        "device_id", "trace_time_sc", "src_type_code", "src_energy_dir",
        "src_exp", "src_units_code"])

    ints_k = String["trace_seq_line", "trace_seq_file", "rec_no", "channel_no",
      "energy_src_pt", "ensemble_no", "trace_in_ensemble", "src-rec_dist",
      "rec_ele", "src_ele", "src_dep", "rec_datum_ele",
      "src_datum_ele", "src_water_dep", "rec_water_dep", "src_x",
      "src_y", "rec_x", "rec_y"]
    append!(ints_k, passcal ? String["samp_rate", "num_samps", "max", "min"] :
      String["cdp_x", "cdp_y", "inline_3d", "crossline_3d", "shot_point",
      "trans_mant", "unassigned_1", "unassigned_2"] )
    misc = Dict{String,Any}()
    [misc[shorts_k[i]] = shorts[i] for i in 1:length(shorts_k)]
    [misc[ints_k[i]] = ints[i] for i in 1:length(ints_k)]
    misc["scale_fac"] = gain
    if passcal
      sta = String(chars[1:6])
      misc["inst_no"] = inst_no
      misc["sensor_serial"] = String(chars[7:14])
    end
    misc["station_name"] = sta
    misc["channel_name"] = cha
  end

  # Trace info
  ts = mktime(shorts[41], shorts[42], shorts[43], shorts[45], shorts[45], zero(Int16)) +
        shorts[53] + 1000*sum(shorts[15:17])
  loc = GeoLoc()
  loc.lat = lat
  loc.lon = lon
  loc.el = el

  C = SeisChannel()
  setfield!(C, :name, id)
  setfield!(C, :id, id)
  setfield!(C, :loc, loc)
  setfield!(C, :gain, gain)
  setfield!(C, :fs, fs)
  mk_t!(C, length(x), ts)
  setfield!(C, :x, x)
  if passcal == false
    setfield!(C, :units, get(segy_units, Int16(shorts[56]), ""))
  end
  if full == true
    setfield!(C, :misc, misc)
  end
  return C
end

function read_segy_file!( S::GphysData,
                          fname::String,
                          ll::UInt8,
                          passcal::Bool,
                          memmap::Bool,
                          full::Bool,
                          swap::Bool,
                          strict::Bool)

  f = memmap ? IOBuffer(Mmap.mmap(fname)) : open(fname, "r")
  trace_fh = Array{Int16, 1}(undef, 3)
  if passcal == true
    C = do_trace(f, true, full, ll, swap, trace_fh)
    add_chan!(S, C, strict)
    close(f)
  else
    shorts  = getfield(BUF, :int16_buf)

    # File headers
    filehdr       = fastread(f, 3200)
    jobid         = bswap(fastread(f, Int32))
    lineid        = bswap(fastread(f, Int32))
    reelid        = bswap(fastread(f, Int32))
    fast_readbytes!(f, BUF.buf, 48)
    fillx_i16_be!(shorts, BUF.buf, 24, 0)
    fastskip(f, 240)

    # Some early sample files had these in little-endian byte order
    for i = 25:27
      shorts[i] = bswap(fastread(f, Int16))
    end
    fastskip(f, 94)

    # Process file header
    nh = max(zero(Int16), getindex(shorts, 27))
    if full == false
      fastskip(f, 3200*nh)
    else
      exthdr = Array{String,1}(undef, nh)
      [exthdr[i] = fastread(f, 3200) for i in 1:nh]
      fhd = Dict{String,Any}(
              zip(String["ntr", "naux", "filedt", "origdt", "filenx",
                         "orignx", "fmt", "cdpfold", "trasort", "vsum",
                         "swst", "swen0", "swlen", "swtyp", "tapnum",
                         "swtapst", "swtapen", "taptyp", "corrtra", "bgainrec",
                         "amprec", "msys", "zupdn", "vibpol", "segyver",
                         "isfixed", "n_exthdr"], shorts[1:27])
                            )
      fhd["jobid"]  = jobid
      fhd["lineid"] = lineid
      fhd["reelid"] = reelid
      fhd["filehdr"] = filehdr
      fhd["exthdr"] = exthdr
    end

    trace_fh[1] = getindex(shorts,3)
    trace_fh[2] = getindex(shorts,5)
    trace_fh[3] = getindex(shorts,7)

    # Channel headers
    nt = shorts[1]
    for i = 1:nt
      # "swap" is always true for valid SEG Y data
      C = do_trace(f, false, full, ll, true, trace_fh)
      j = add_chan!(S, C, strict)
      if full == true
        merge!(S.misc[j], fhd)
      end
    end
    close(f)
  end
  resize!(BUF.buf, 65535)
  return S
end

# ============================================================================

"""
    segyhdr(f[; passcal=false, ll=LL, swap=false])

Print formatted, sorted SEG-Y headers of file `f` to stdout. Use keyword
`passcal=true` for PASSCAL/NMT modified SEG Y; use `swap=true` for big-endian
PASSCAL. See SeisIO `read_data` documentation for `ll` codes.
"""
function segyhdr(fname::String; ll::UInt8=0x00, passcal::Bool=false, swap::Bool=false)
  if passcal
    seis = read_data("passcal", fname::String, full=true, ll=ll, swap=swap)
  else
    seis = read_data("segy", fname::String, ll=ll, full=true)
  end
  if passcal
    printstyled(stdout, @sprintf("%20s: %s\n", "PASSCAL SEG-Y FILE", realpath(fname)), color=:green, bold=true)
    D = getindex(getfield(seis, :misc),1)
    for k in sort(collect(keys(D)))
      @printf(stdout, "%20s: %s\n", k, string(get(D, k, "")))
    end
  else
    p = 1; w = 32; W = displaysize(stdout)[2]-2
    S = fill("", length(seis.misc[1])+1)
    printstyled(stdout, @sprintf("%20s: %s\n", "SEG-Y FILE", realpath(fname)), color=:green, bold=true)
    for i = 1:seis.n
      if p > 1
        s = @sprintf("%20i/%i", i, seis.n)
      else
        s = @sprintf(" Trace # %11i/%i", i, seis.n)
      end
      S[1] *= s
      S[1] *= " "^(w-length(s))
      D = getindex(getfield(seis, :misc),1)
      for (j,k) in enumerate(sort(collect(keys(D))))
        if k == "exthdr" || k == "filehdr"
          val = get(seis.misc[i], k, "")
          if isempty(val)
            s = "(empty)"
          else
            s = length(val) > 8 ? String(val[1:8])*"…" : String(val)
          end
        else
          s = string(get(seis.misc[i], k, ""))
        end
        filler = " "^max(0, 10-length(s))
        S[j+1] *= @sprintf("%20s: %s%s", k, s, filler)
      end
      if (p+2*w > W) || (i == seis.n)
        printstyled(stdout, S[1]*"\n", color=:yellow, bold=true)
        [println(S[j]) for j=2:length(S)]
        println("")
        S = fill("", length(seis.misc[i])+1)
        p = 1
      else
        p += w
      end
    end
    printstyled(stdout, @sprintf("%20s\n","END OF RECORD"), color=208, bold=true)
  end
  return nothing
end


text:
using Base.Test
using StochasticProcesses

let p = PoissonProcess(5), t=linspace(0, 1, 10)
  sim(p, t, 3)
  cumsim(p, t, 3)
  distribution(p, 10)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "libntl"
version = v"10.5.0"

# Collection of sources required to build libntl
sources = [
    "https://www.shoup.net/ntl/ntl-10.5.0.tar.gz" =>
    "b90b36c9dd8954c9bc54410b1d57c00be956ae1db5a062945822bbd7a86ab4d2",
    "./bundled"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd ntl-10.5.0/
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/configure.patch
cd $WORKSPACE/srcdir
CC_OLD=$CC
CC=$CC_FOR_BUILD
CXX_OLD=$CXX
CXX=$CXX_FOR_BUILD
LIBTOOL_OLD=$LIBTOOL
LIBTOOL=/usr/bin/libtool
LD_OLD=$LD
LD=$LD_FOR_BUILD
LDFLAGS_OLD=$LDFLAGS
LDFLAGS=""
wget https://gmplib.org/download/gmp/gmp-6.1.2.tar.bz2
tar -xvf gmp-6.1.2.tar.bz2
cd gmp-6.1.2
CC=$CC_FOR_BUILD CXX=$CXX_FOR_BUILD ./configure --enable-cxx --enable-shared
LD=$LD_FOR_BUILD LDFLAGS="" CC=$CC_FOR_BUILD CXX=$CXX_FOR_BUILD make -j${nproc}
make install
CC=$CC_OLD
CXX=$CXX_OLD
LIBTOOL=$LIBTOOL_OLD
LD=$LD_OLD
LDFLAGS=$LDFLAGS_OLD
cd $WORKSPACE/srcdir
export CPP_FLAGS=-I$prefix/include
export LD_LIBRARY_PATH=$target/lib:$LD_LIBRARY_PATH
export LDFLAGS=-Wl,-rpath,$prefix/lib
cd $WORKSPACE/srcdir
cd ntl-10.5.0/
cd src
./configure PREFIX=$prefix DEF_PREFIX=$prefix SHARED=on NTL_THREADS=off NTL_EXCEPTIONS=off NTL_GMP_LIP=on NATIVE=off CXXFLAGS="-I$prefix/include -march=core2" TUNE=x86 CXX_FOR_BUILD="$CXX_FOR_BUILD" LD_FOR_BUILD="$LD_FOR_BUILD" CXX="$CXX" HOST=$target
make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    MacOS(:x86_64),
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libntl", Symbol(""))
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaMath/GMPBuilder/releases/download/v6.1.2-2/build_GMP.v6.1.2.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
module UHDBindings

using Printf
using Preferences
using Libdl

# ----------------------------------------------------
# --- Loading UHD library
# ----------------------------------------------------
# We use Yggdrasil to load the artifact responsible for UHD.
# It is possible to use a local install instead of the proposed UHD version
"""
Change UHD driver provider. Support "yggdrasil" to use shipped jll file or "local" to use custom installed library
set_provider("yggdrasil")
or
set_provider("local")
"""
function set_provider(new_provider::String)
    if !(new_provider in ("yggdrasil","default", "local"))
        throw(ArgumentError("Invalid provider: \"$(new_provider)\""))
    end
    # Set it in our runtime values, as well as saving it to disk
    @set_preferences!("provider" => new_provider)
    @info("New provider set; restart your Julia session for this change to take effect!")
end
function get_provider()
    return @load_preference("provider","yggdrasil")
end
const uhd_provider = get_provider()
@static  if uhd_provider == "yggdrasil" || uhd_provider =="local"
    # --- Using Yggdrasil jll file
    using USRPHardwareDriver_jll
    const libUHD = USRPHardwareDriver_jll.libuhd
end
@static if uhd_provider == "local"
    # --- Using local install, assuming it works
    libUHD_system_h = dlopen("libuhd", false);
    const libUHD = dlpath(libUHD_system_h)
end
# ----------------------------------------------------
# --- Bindings, structure and low level functions
# ----------------------------------------------------
# --- Including the file
include("Bindings.jl");
# Export necessary structures for high API level managment
export Timestamp
export UHDBinding
# Exporting printing macros
export @infotx, @warntx;
export @inforx, @warnrx;

# ----------------------------------------------------
# --- Finding devices
# ----------------------------------------------------
include("Find.jl")
export uhd_find_devices

# ----------------------------------------------------
# --- Receiver Configuration
# ----------------------------------------------------
# All structures and functions for the Rx side
include("Rx.jl");
# Structures
export UHDRx
# Export functions
export recv,recv!;
export populateBuffer!
export getError, getTimestamp
export restartStreamer

# ----------------------------------------------------
# --- Transmitter Configuration
# ----------------------------------------------------
# All structures and functions for the Tx side
include("Tx.jl");
# Structures
export UHDTx
# Export functions
export send;
"""
Init the core parameter of the radio in Tx or in Rx mode and initiate RF parameters
# --- Syntax
openUHD(mode,sysImage,carrierFreq,samplingRate,txGain,antenna="RX2")
# --- Input parameters
- mode 			: String to open radio in "Tx" (transmitter) or in "Rx" (receive) mode
- carrierFreq	: Desired Carrier frequency [Union{Int,Float64}]
- samplingRate	: Desired bandwidth [Union{Int,Float64}]
- gain		: Desired Gain [Union{Int,Float64}]
- antenna		: Desired Antenna alias Dict{Symbol,Vector[String]} (e.g. Dict(:Rx=>["RX2"],:Tx=>["TRX"]))
Keywords=
- args	  : String with the additionnal load parameters (for instance, path to the FPHGA image) [String]
# --- Output parameters
- uhd		  	: UHD object [UHDBinding]
"""
function openUHD(carrierFreq, samplingRate, gain;args="",channels=[0],antennas=Dict(:Rx=>["RX2"],:Tx=>["TRX"]),cpu_format="fc32",otw_format="sc16",subdev="",nbAntennaRx=1,nbAntennaTx=1,bypassStreamer=false)
	# ----------------------------------------------------
	# --- Handler
	# ----------------------------------------------------
	addressUSRP = Ref{uhd_usrp_handle}();
	# --- Cal the init
    @assert_uhd uhd_usrp_make(addressUSRP,args)
	# --- Get the usable object
	pointerUSRP = addressUSRP[];
    # ----------------------------------------------------
    # --- Manage dictionnaries for custom option passing
    # ----------------------------------------------------
    aTx = parseDict(antennas,:Tx,["TRX"])
    aRx = parseDict(antennas,:Rx,["RX2"])
	# ----------------------------------------------------
	# --- Set Rx stage
	# ----------------------------------------------------
    rx = openUHDRx(pointerUSRP,carrierFreq, samplingRate, gain;channels,antennas=aRx,args,cpu_format,otw_format,subdev,nbAntennaRx,bypassStreamer)
	# ----------------------------------------------------
	# --- Set Tx stage
	# ----------------------------------------------------
    tx = openUHDTx(pointerUSRP,carrierFreq, samplingRate, gain;channels,antennas=aTx,args,cpu_format,otw_format,subdev,nbAntennaTx,bypassStreamer)
	# ----------------------------------------------------
	# --- Create radio
	# ----------------------------------------------------
	uhdBinding = UHDBinding(addressUSRP,rx,tx);
	return uhdBinding;
end
export openUHD;

# Be able to get a default value
parseDict(d,key,val=0) = (haskey(d,key) ? d[key] : val)

"""
Close the USRP device (Rx or Tx mode) and release all associated objects

# --- Syntax

close(uhd)
# --- Input parameters
- uhd	: UHD object [UHDBinding]
# --- Output parameters
- []
"""
function Base.close(uhdBinding::UHDBinding)
	# --- Close Rx and Tx streams
	close(uhdBinding.rx);
	close(uhdBinding.tx);
	# --- Close the USRP main object
    @assert_uhd uhd_usrp_free(uhdBinding.addressUSRP);
	# --- Print a flag
	print("\n");
	@info "USRP device is now closed.";
end

"""
Print the radio configuration

# --- Syntax

print(radio,chan=0)
# --- Input parameters
- radio		: UHD object [Union{UHDBinding,UHDTx,UHDRx}]
- chan : Channel index to display (default 0)
# --- Output parameters
- []
"""
function Base.print(radio::UHDBinding,chan=0)
	# --- Print the configuration of Tx and Rx
	print(radio.rx,chan);
	print(radio.tx,chan);
end

# All functions are defined @UHDTx or @UHDRx level, we should define configuration functions @UHDBinding level
# In such case, we will both apply configuration at Rx and Tx sides.

"""
Update sampling rate of current radio device, and update radio object with the new obtained sampling frequency. If the input parameter is the UHDBinding object, the desired sampling frequency will be applied on both Rx and Tx sides.
If the input is a [UHDRx] or a [UHDTx] object, it updates only the Rx or Tx sampling frequency

# --- Syntax

updateSamplingRate!(radio,samplingRate,chan=0)
# --- Input parameters
- radio	  : UHD device [Union{UHDBinding,UHDRx,UHDTx}]
- samplingRate	: New desired sampling rate
- chan : Channel index to use (default 0)
# --- Output parameters
-
"""
function updateSamplingRate!(radio::UHDBinding,samplingRate,chan=0)
	@sync updateSamplingRate!(radio.rx,samplingRate,chan);
	@sync updateSamplingRate!(radio.tx,samplingRate,chan);
end

"""
Update carrier frequency of current radio device, and update radio object with the new obtained carrier frequency. If the input parameter is the UHDBinding object, the desired carrier frequency will be applied on both Rx and Tx sides.
If the input is a [UHDRx] or a [UHDTx] object, it updates only the Rx or Tx carrier frequency

# --- Syntax

updateCarrierFreq!(radio,carrierFreq)
# --- Input parameters
- radio	  : UHD device [Union{UHDBinding,UHDRx,UHDTx}]
- carrierFreq	: New desired carrier frequency
- chan : Channel index to use (default 0)
# --- Output parameters
-
"""
function updateCarrierFreq!(radio::UHDBinding,carrierFreq,chan=0)
	@sync updateCarrierFreq!(radio.rx,carrierFreq,chan);
	@sync updateCarrierFreq!(radio.tx,carrierFreq,chan);
end

"""
Update gain of current radio device, and update radio object with the new obtained gain. If the input parameter is the UHDBinding object, the desired gain will be applied on both Rx and Tx sides.
If the input is a [UHDRx] or a [UHDTx] object, it updates only the Rx or Tx gain

# --- Syntax

updateGain!(radio,gain)
# --- Input parameters
- radio	  : UHD device [Union{UHDBinding,UHDRx,UHDTx}]
- gain	: New desired gain
- chan : Channel index to use (default 0)
# --- Output parameters
-
"""
function updateGain!(radio::UHDBinding,gain,chan=0)
	@sync updateGain!(radio.rx,gain,chan);
	@sync updateGain!(radio.tx,gain,chan);
end

# When given to UHDBinding, recv and send will dispatch to the appropriate substructure
# Recv
recv(radio::UHDBinding,nbSamples)  = recv(radio.rx,nbSamples);
recv!(sig::Vector{Vector{T}},radio::UHDBinding;kwargs...) where T = recv!(sig,radio.rx;kwargs...);
recv!(sig::Vector{T},radio::UHDBinding;kwargs...) where T = recv!([sig],radio.rx;kwargs...);
# Send
send(radio::UHDBinding,params...) = send(radio.tx,params...);

# Config can laos be dispatch
uhd_usrp_create_stream(radio::UHDBinding;kwargs...) = uhd_usrp_create_stream(radio.rx;kwargs...)

""""
Returns the  internal radio buffer size
"""
function getBufferSize(radio)
    return radio.rx.packetSize
end

# ----------------------------------------------------
# --- Common functions and structures
# ----------------------------------------------------
export updateSamplingRate!
export updateGain!
export updateCarrierFreq!
export print;
export close;
export getBufferSize

end # module


text:
module DailyMed

export applicationnumbers, drugclasses, drugnames, ndcs, rxcuis, spls, spls_setid, history, media, ncds, packaging, uniis

using HTTP
using EzXML
using RxNav

const BASEURL = "https://dailymed.nlm.nih.gov/dailymed/services/v2/"

const METATAGS = ["total_elements", "elements_per_page", "total_pages", "current_page",
    "current_url", "previous_page", "previous_page_url", "next_page", "next_page_url", "db_published_date"]

"""
    `allsetids(resource)`

Given a resource which is either an RxCUI id or the name of a drug, return a vector
of all of the DailyMed Set ID identifiers for the resouce. These can then be used for
the functions `history`, `media`, `ndcs`, or `packaging`.
"""
function allsetids(resource::String)
    rxcu = RxNav.is_in_rxcui_format(resource) ? resource : RxNav.rcui(resource)
    setids = [x.setid for x in first(spls(extra = ["rxcui" => rxcu]))]
    return setids
end
allsetids(resource) = allsetids(string(resource))

"""
    `dailymed(restfunc, extra)`

Get and partially parse data from the url formed by `BASEURL * restfunc *`` expanded `extra` args
"""
function dailymed(restfunc, extra)
    url = BASEURL * restfunc
    if any(x -> x[2] != "", extra)
        url *= "?"
        for p in extra
            if p[2] != ""
                url *= (url[end] == '?' ? "&" : "") * "$(p[1])=$(p[2])"
            end
        end
    end
    try
        req = HTTP.request("GET", url)
        doc = parsexml(String(req.body)).root
        meta = findfirst("//metadata", doc)
        metatags = Dict{String, String}()
        for tag in METATAGS
            if (x = findfirst(tag, meta)) != nothing
                metatags[tag] = nodecontent(x)
            end
        end
        return doc, metatags
    catch y
        # warn if it is likely to have been our error, just return if it was internal to REST code on server
        !contains(string(y), "Internal Server Error") && @warn y
        return parsexml("<root></root>"), Dict()
    end
end

"""
    `applicationnumbers(; extra = [])`

Returns a list of all NDA numbers.

extra is optional. If provided it should be a `Dict` or list of string `Pair`s,
and can be "application_number", "marketing_category_code", "setid", "pagesize", "page"
"""
function applicationnumbers(; extra = [])
    appnums, d, metadict = String[], Dict(extra), Dict()
    while true
        doc, metadict = dailymed("applicationnumbers.xml", d)
        append!(appnums, nodecontent.(findall("//application_number", doc)))
        nextpage = tryparse(Int, get(metadict, "next_page", ""))
        nextpage == nothing && break
        d["page"] = string(nextpage)
    end
    return appnums, metadict
end

"""
    `drugclasses(; extra = [])`

Returns a list of all drug classes associated with at least one SPL in the
Pharmacologic Class Indexing Files.

extra is optional. If provided it should be a `Dict` or list of string `Pair`s,
and can be "drug_class_code", "drugclass_coding_system", "code_class_type",
"class_name", "unii_code", "pagesize", "page"
"""
function drugclasses(; extra = [])
    dctups, d, metadict = NamedTuple[], Dict(extra), Dict()
    while true
        doc, metadict = dailymed("drugclasses.xml", d)
        for dcl in findall("//drugclass", doc)
            push!(dctups, (name = nodecontent(findfirst("name", dcl)),
                           code = nodecontent(findfirst("code", dcl))))
        end
        nextpage = tryparse(Int, get(metadict, "next_page", ""))
        nextpage == nothing && break
        d["page"] = string(nextpage)
    end
    return dctups, metadict
end

"""
    `drugnames(; extra = [])`

Returns a list of all drug names. A very large string vector is returned, and the metadata

extra is optional. If provided it should be a `Dict` or list of string `Pair`s,
and can be "name_type", "manufacturer", "pagesize", "page"
"""
function drugnames(; extra = [])
    dnames, d, metadict = String[], Dict(extra), Dict()
    while true
        doc, metadict = dailymed("drugnames.xml", d)
        append!(dnames, nodecontent.(findall("//drug_name", doc)))
        nextpage = tryparse(Int, get(metadict, "next_page", ""))
        nextpage == nothing && break
        d["page"] = string(nextpage)
    end
    return dnames, metadict
end

"""
    `function ndcs(; extra = [])`

Returns a list of all NDC codes.

`extra` is optional. If provided it should be a `Dict` or list of string `Pair`s,
and can be "pagesize" or "page"
"""
function ndcs(; extra = [])
    codes, d, metadict = String[], Dict(extra), Dict()
    while true
        doc, metadict = dailymed("ndcs.xml", d)
        append!(codes, nodecontent.(findall("//ndc", doc)))
        nextpage = tryparse(Int, get(metadict, "next_page", ""))
        nextpage == nothing && break
        d["page"] = string(nextpage)
    end
    return codes, metadict
end

"""
    function rxcuis(; extra = [])

Returns a list of all product-level RxCUIs.

extra is optional. If provided it should be a `Dict` or list of string `Pair`s,
and can be "rxtty", "rxstring", "rxcui", "pagesize", "page"
"""
function rxcuis(; extra = [])
    dctups, d, metadict = NamedTuple[], Dict(extra), Dict()
    while true
        doc, metadict = dailymed("rxcuis.xml", d)
        for dcl in findall("//rxconcept", doc)
            push!(dctups, (rxcui = nodecontent(findfirst("rxcui", dcl)),
                rxstring = nodecontent(findfirst("rxstring", dcl)),
                rxtty = nodecontent(findfirst("rxtty", dcl))))
        end
        nextpage = tryparse(Int, get(metadict, "next_page", ""))
        nextpage == nothing && break
        d["page"] = string(nextpage)
    end
    return dctups, metadict
end

"""
    function spls(; extra = [])

Returns a list of all SPLs.

extra is optional. If provided it should be a `Dict` or list of string `Pair`s,
and can be "application_number", "boxed_warning", "dea_schedule_code", "doctype",
"drug_class_code", "drugclass_coding_system", "drug_name", "name_type", "labeler",
"manufacturer", "marketing_category_code", "ndc", "published_date",
"published_date_comparison", "rxcui", "setid", "unii_code", "pagesize", "page"
"""
function spls(; extra = [])
    dctups, d, metadict = NamedTuple[], Dict(extra), Dict()
    while true
        doc, metadict = dailymed("spls.xml", d)
        for dcl in findall("//spl", doc)
            push!(dctups, (setid = nodecontent(findfirst("setid", dcl)),
                spl_version = nodecontent(findfirst("spl_version", dcl)),
                title = nodecontent(findfirst("title", dcl)),
                published_date = nodecontent(findfirst("published_date", dcl))))
        end
        nextpage = tryparse(Int, get(metadict, "next_page", ""))
        nextpage == nothing && break
        d["page"] = string(nextpage)
    end
    return dctups, metadict
end

"""
    function function spls_setid(setid)

Returns an SPL document for specific SET ID.
"""
function spls_setid(setid)
    url, metadict = BASEURL * "spls/$setid.xml", Dict()
    try
        req = HTTP.request("GET", url)
        return String(req.body), metadict
    catch y
        return "<xml>$y</xml>", metadict
    end
end

"""
    function history(setid; extra)

Returns version history for specific SET ID.

extra is optional. If provided it should be a `Dict` or list of string `Pair`s,
and can be "pagesize", "page"
"""
function history(setid; extra = [])
    dctups, d, metadict = NamedTuple[], Dict(extra), Dict()
    while true
        doc, metadict = dailymed("spls/$(setid)/history.xml", d)
        for dcl in findall("//history/history_entry", doc)
            push!(dctups, (spl_version = nodecontent(findfirst("spl_version", dcl)),
                           published_date = nodecontent(findfirst("published_date", dcl))))
        end
        nextpage = tryparse(Int, get(metadict, "next_page", ""))
        nextpage == nothing && break
        d["page"] = string(nextpage)
    end
    return dctups, metadict
end

"""
    function media(setid; extra = [])

Returns links to all media for specific SET ID.

`extra` is optional. If provided it should be a `Dict` or list of string `Pair`s,
and can be "pagesize", "page"
"""
function media(setid; extra = [])
    dctups, d, metadict = NamedTuple[], Dict(extra), Dict()
    while true
        doc, metadict = dailymed("spls/$(setid)/media.xml", d)
        for dcl in findall("//file", doc)
            push!(dctups, (name = nodecontent(findfirst("name", dcl)),
                mime_type = nodecontent(findfirst("mime_type", dcl)),
                url = nodecontent(findfirst("url", dcl))))
        end
        nextpage = tryparse(Int, get(metadict, "next_page", ""))
        nextpage == nothing && break
        d["page"] = string(nextpage)
    end
    return dctups, metadict
end

"""
    function ndcs(setid; extra = [])

Returns all ndcs for specific SET ID.

`extra` is optional. If provided it should be a `Dict` or list of string `Pair`s,
and can be "pagesize", "page"
"""
function ndcs(setid; extra = [])
    nds, d, metadict = String[], Dict(extra), Dict()
    while true
        doc, metadict = dailymed("spls/$setid/ndcs.xml", d)
        append!(nds, nodecontent.(findall("//ndcs/ndc", doc)))
        nextpage = tryparse(Int, get(metadict, "next_page", ""))
        nextpage == nothing && break
        d["page"] = string(nextpage)
    end
    return nds, metadict
end

"""
    function packaging(setid; extra = [])

Return the XML string for the packaging of the item with the given setid.
The packaging XML is highly variable in labeling and may be deeply nested, so an array
or tuple is not computed, but instead the XML itself is returned.

`extra` is optional. If provided it should be a `Dict` or list of string `Pair`s,
and can be "pagesize", "page"
"""
function packaging(setid; extra = [])
    allpages, d, metadict = "", Dict(extra), Dict()
    try
        while true
            doc, metadict = dailymed("spls/$setid/packaging.xml", d)
            allpages *= string(findfirst("//products", doc))
            nextpage = tryparse(Int, get(metadict, "next_page", ""))
            nextpage == nothing && break
            d["page"] = string(nextpage)
        end
    catch y
        @warn y
    end
    return allpages, metadict
end

"""
    function uniis(; extra = [])

Returns a list of all UNIIs.

`extra` is optional. If provided it should be a `Dict` or list of string `Pair`s,
and can be "active_moiety", "drug_class_code", "drug_class_coding_system",
"rxcui", "unii_code", "pagesize", "page"
"""
function uniis(; extra = [])
    dctups, d, metadict = NamedTuple[], Dict(extra), Dict()
    while true
        doc, metadict = dailymed("uniis.xml", d)
        for dcl in findall("//unii", doc)
            push!(dctups, (unii_code = nodecontent(findfirst("unii_code", dcl)),
                active_moiety = nodecontent(findfirst("active_moiety", dcl))))
        end
        nextpage = tryparse(Int, get(metadict, "next_page", ""))
        nextpage == nothing && break
        d["page"] = string(nextpage)
    end
    return dctups, metadict
end

end # module


text:
nfactors(d::ProductDomain) = length(components(d))
factors(d::ProductDomain) = components(d)
factor(d::ProductDomain,k::Integer) = component(d,k)

canonicaldomain(d::ProductDomain) = ProductDomain(map(canonicaldomain,factors(d))...)

# product domains are their own canonical domain
for OP in (:fromcanonical,:tocanonical)
    @eval begin
        $OP(d::ProductDomain, x::Vec) = Vec(map($OP,factors(d),x)...)
        $OP(d::ProductDomain, x::Vec{2}) = Vec($OP(first(factors(d)), first(x)), $OP(last(factors(d)), last(x)))
    end
end

function pushappendpts!(ret, xx, pts)
    if isempty(pts)
        push!(ret,Vec(xx...))
    else
        for x in pts[1]
            pushappendpts!(ret,(xx...,x...),pts[2:end])
        end
    end
    ret
end

function checkpoints(d::ProductDomain)
    pts = checkpoints.(factors(d))
    ret=Vector{Vec{sum(dimension.(factors(d))),float(promote_type(eltype.(eltype.(factors(d)))...))}}(undef, 0)

    pushappendpts!(ret,(),pts)
    ret
end

function points(d::ProductDomain,n::Tuple)
    @assert length(factors(d)) == length(n)
    pts=map(points,factors(d),n)
    ret=Vector{Vec{length(factors(d)),mapreduce(eltype,promote_type,factors(d))}}(undef, 0)
    pushappendpts!(ret,Vec(x),pts)
    ret
end

reverseorientation(d::ProductDomain) = ProductDomain(map(reverseorientation, factors(d)))

domainscompatible(a::ProductDomain,b::ProductDomain) =
                        length(factors(a))==length(factors(b)) &&
                        all(map(domainscompatible,factors(a),factors(b)))


text:
module ActuaryUtilities

using Dates
using ForwardDiff
using QuadGK
using Yields
import Roots
import StatsBase

include("financial_math.jl")
include("risk_measures.jl")
include("derivatives.jl")

# https://timholy.github.io/SnoopCompile.jl/stable/snoopi_deep_parcel/
if Base.VERSION >= v"1.4.2"
    include("precompile.jl")
    _precompile_()
end

"""
    Years_Between(d1::Date, d2::Date)

Compute the number of integer years between two dates, with the
first date typically before the second. Will return negative number if
first date is after the second. Use third argument to indicate if calendar
anniversary should count as a full year.

# Examples
```jldoctest
julia> d1 = Date(2018,09,30);

julia> d2 = Date(2019,09,30);

julia> d3 = Date(2019,10,01);

julia> years_between(d1,d3)
1
julia> years_between(d1,d2,false) # same month/day but `false` overlap
0
julia> years_between(d1,d2) # same month/day but `true` overlap
1
julia> years_between(d1,d2) # using default `true` overlap
1
```
"""
function years_between(d1::Date, d2::Date, overlap = true)
    iy, im, id = Dates.year(d1), Dates.month(d1), Dates.day(d1)
    vy, vm, vd = Dates.year(d2), Dates.month(d2), Dates.day(d2)
    dur = vy - iy
    if vm == im
        if overlap
            if vd >= id
                dur += 1
            end
        else
            if vd > id
                dur += 1
            end
        end
    elseif vm > im
        dur += 1
    end

    return dur - 1
end

"""
    duration(d1::Date, d2::Date)

Compute the duration given two dates, which is the number of years
since the first date. The interval `[0,1)` is defined as having
duration `1`. Can return negative durations if second argument is before the first.

```jldoctest
julia> issue_date  = Date(2018,9,30);

julia> duration(issue_date , Date(2019,9,30) )
2
julia> duration(issue_date , issue_date)
1
julia> duration(issue_date , Date(2018,10,1) )
1
julia> duration(issue_date , Date(2019,10,1) )
2
julia> duration(issue_date , Date(2018,6,30) )
0
julia> duration(Date(2018,9,30),Date(2017,6,30))
-1
```

"""
function duration(issue_date::Date, proj_date::Date)
    return years_between(issue_date, proj_date, true) + 1
end

"""
    accum_offset(x; op=*, init=1.0)

A shortcut for the common operation wherein a vector is scanned with an operation, but has an initial value and the resulting array is offset from the traditional accumulate.

This is a common pattern when calculating things like survivorship given a mortality vector and you want the first value of the resulting vector to be `1.0`, and the second value to be `1.0 * x[1]`, etc.

Two keyword arguments:

- `op` is the binary (two argument) operator you want to use, such as `*` or `+`
- `init` is the initial value in the returned array

# Examples

```julia=repl
julia> accum_offset([0.9, 0.8, 0.7])
3-element Array{Float64,1}:
 1.0
 0.9
 0.7200000000000001

julia> accum_offset(1:5) # the product of elements 1:n, with the default `1` as the first value
5-element Array{Int64,1}:
  1
  1
  2
  6
 24

julia> accum_offset(1:5,op=+)
5-element Array{Int64,1}:
  1
  2
  4
  7
 11

```

"""
function accum_offset(x; op=*,init=1.0)
    xnew = similar(x)
    xnew[1] = init
    for i in 2:length(x)
      xnew[i] = op(xnew[i-1],x[i-1])
    end
    return xnew
end

export years_between, duration,
    irr, internal_rate_of_return,
    pv, present_value, price, present_values,
    breakeven,moic,
    accum_offset,
    Macaulay,Modified,DV01,KeyRatePar,KeyRateZero,KeyRate,duration, convexity,
    VaR,ValueAtRisk,CTE,ConditionalTailExpectation,ExpectedShortfall,
    eurocall, europut

end # module

max_stars_repo_name: alecloudenback/ActuaryUtilities

max_issues_repo_name: alecloudenback/ActuaryUtilities

max_forks_repo_name: alecloudenback/ActuaryUtilities


text:
# Autogenerated wrapper script for P4est_jll for x86_64-w64-mingw32
export libp4est, libsc

using Zlib_jll
JLLWrappers.@generate_wrapper_header("P4est")
JLLWrappers.@declare_library_product(libp4est, "libp4est-0.dll")
JLLWrappers.@declare_library_product(libsc, "libsc-0.dll")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_library_product(
        libp4est,
        "bin\\libp4est-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsc,
        "bin\\libsc-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# array.jl - arrays of images

export jim

using UnitfulRecipes
using Plots: plot!
import Plots # gui
using MosaicViews: mosaic
#using MIRTjim: prompt

"""
    jim(z::AbstractArray{<:AbstractArray{<:Number}} ; kwargs...)

Display an array of images.
Same arguments and options as display of a 3D stack of images.
"""
function jim(z::AbstractArray{<:AbstractArray{<:Number}} ;
    gui::Bool = jim_def[:gui],
    prompt::Bool = jim_def[:prompt],
    line3plot = jim_def[:line3plot],
    line3type = jim_def[:line3type],
    ncol::Int = jim_def[:ncol],
    nrow::Int = jim_def[:nrow],
    padval = nothing_else(jim_def[:padval], _mingood(z)),
    mosaic_npad::Int = jim_def[:mosaic_npad],
    fft0::Bool = jim_def[:fft0],
    x::AbstractVector{<:Number} = axes(z[1],1),
    y::AbstractVector{<:Number} = axes(z[1],2),
    xy::Tuple = (),
    xticks = _ticks(x),
    yticks = _ticks(y),
    yflip::Bool = nothing_else(jim_def[:yflip], minimum(y) >= zero(y[1])),
    kwargs...
)

    all(x -> size(x) == size(z[1]), z) || throw("size mismatch")
    all(x -> axes(x) == axes(z[1]), z) || throw("axes mismatch")

    # determine mosaic layout
    n3 = length(z) * prod(size(z[1])[3:end])
    n1,n2 = size(z[1])[1:2]
    if ncol == 0 && nrow == 0
        nrow = n3 * n1 / n2 # wider images means fewer columns
        nrow *= 600/400 # typical plot aspect ratio
        nrow = max(floor(Int, sqrt(nrow)), 1)
        ncol = ceil(Int, n3 / nrow)
    end
    if ncol == 0
        ncol = ceil(Int, n3 / nrow)
    end
    if nrow == 0
        nrow = ceil(Int, n3 / ncol)
    end

    if !yflip # handling yflip with array mosaic is tricky
        z = [reverse(x, dims=2) for x in z]
        yflip = !yflip
    end

    zz = mosaic(vec(z) ; fillvalue = padval, ncol, nrow, npad = mosaic_npad)
    fft0 && @warn("fft0 option ignored for 3D")

    # adjusted x,y for mosaic
    x = range(x[1]; length=size(zz,1), step=x[2]-x[1])
    zyflip = false
    if y[[1,end]] == axes(z[1],2)[[1,end]] # typical 1:ny case
        y = range(y[1]; length=size(zz,2), step=y[2]-y[1])
    else # ig.y setting
        y = range(y[end]; length=size(zz,2), step=y[1]-y[2])
        yflip = false
        zz = reverse(zz, dims=2)
        zyflip = true
        y = reverse(y)
    end
    xy = (x,y)

    p = jim(zz ; transpose = false, xy, xticks, yticks, yflip,
        gui=false, prompt=false, kwargs...)

    if n3 > 1 && line3plot # lines around each subimage
        n1 += mosaic_npad
        n2 += mosaic_npad
        m1 = (1+size(zz,1)) ÷ n1 # add one because of mosaicview non-edge
        m2 = (1+size(zz,2)) ÷ n2
        fx = o -> x[1] + o * (x[2] - x[1])
        fy = o -> y[1] + o * (y[2] - y[1])
        if zyflip
            fy = o -> y[end] - o * (y[2] - y[1])
        end
        function plot_box!(ox, oy)
            plot!(p,
                fx.(ox .+ [0,1,1,0,0] * n1 .- 1),
                fy.(oy .+ [0,0,1,1,0] * n2 .- 1),
                line = jim_def[:line3type], label="",
            )
        end
        for ii=0:n3-1 # boxes around data only
            i1 = mod(ii, m1)
            i2 = ii ÷ m1
            plot_box!(i1*n1, i2*n2)
        end
    end

    plot!()
    gui && Plots.gui()
    prompt && MIRTjim.prompt()
    return p
end


text:
#=

Author: Victoria Gregory

=#
using LaTeXStrings

# Set up the problem and initial guess, solve by VFI
sp = SearchProblem(;w_grid_size=100, pi_grid_size=100)
v_init = zeros(sp.n_w, sp.n_pi) + sp.c / (1 - sp.bet)
f(x) = bellman_operator(sp, x)
v = compute_fixed_point(f, v_init)
policy = get_greedy(sp, v)

# Make functions for the linear interpolants of these
vf = extrapolate(interpolate((sp.w_grid, sp.pi_grid), v, Gridded(Linear())),
                 Flat())
pf = extrapolate(interpolate((sp.w_grid, sp.pi_grid), policy,
                 Gridded(Linear())), Flat())

function plot_value_function(;w_plot_grid_size::Integer=100,
                             pi_plot_grid_size::Integer=100)
    pi_plot_grid = linspace(0.001, 0.99, pi_plot_grid_size)
    w_plot_grid = linspace(0, sp.w_max, w_plot_grid_size)
    Z = [vf[w_plot_grid[j], pi_plot_grid[i]]
            for j in 1:w_plot_grid_size, i in 1:pi_plot_grid_size]
    p = contour(pi_plot_grid, w_plot_grid, Z, levels=12, alpha=0.6, fill=true)
    plot!(xlabel=L"$\pi$", ylabel="wage", xguidefont=font(12))
    return p
end

function plot_policy_function(;w_plot_grid_size::Integer=100,
                             pi_plot_grid_size::Integer=100)
    pi_plot_grid = linspace(0.001, 0.99, pi_plot_grid_size)
    w_plot_grid = linspace(0, sp.w_max, w_plot_grid_size)
    Z = [pf[w_plot_grid[j], pi_plot_grid[i]]
            for j in 1:w_plot_grid_size, i in 1:pi_plot_grid_size]
    p = contour(pi_plot_grid, w_plot_grid, Z, levels=1, alpha=0.6, fill=true)
    plot!(xlabel=L"$\pi$", ylabel="wage", xguidefont=font(12), cbar=false)
    annotate!(0.4, 1.0, "reject")
    annotate!(0.7, 1.8, "accept")
    return p
end