text:
module Stat

import Gadfly
import StatsBase
import Contour
using Color
using Compose
using DataArrays
using DataStructures
using Hexagons
using Loess

import Gadfly: Scale, Coord, element_aesthetics, default_scales, isconcrete,
               nonzero_length, setfield!
import KernelDensity
import Distributions: Uniform, Distribution, qqbuild
import Iterators: chain, cycle, product, partition, distinct

include("bincount.jl")

# Apply a series of statistics.
#
# Args:
#   stats: Statistics to apply in order.
#   scales: Scales used by the plot.
#   aes: A Aesthetics instance.
#
# Returns:
#   Nothing, modifies aes.
#
function apply_statistics(stats::Vector{Gadfly.StatisticElement},
                          scales::Dict{Symbol, Gadfly.ScaleElement},
                          coord::Gadfly.CoordinateElement,
                          aes::Gadfly.Aesthetics)
    for stat in stats
        apply_statistic(stat, scales, coord, aes)
    end
    nothing
end

immutable Nil <: Gadfly.StatisticElement
end

const nil = Nil

immutable Identity <: Gadfly.StatisticElement
end

function apply_statistic(stat::Identity,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    nothing
end

const identity = Identity

immutable HistogramStatistic <: Gadfly.StatisticElement
    minbincount::Int
    maxbincount::Int
    orientation::Symbol

    function HistogramStatistic(; bincount=nothing,
                                  minbincount=3,
                                  maxbincount=150,
                                  orientation::Symbol=:vertical)
        if bincount != nothing
            new(bincount, bincount, orientation)
        else
            new(minbincount, maxbincount, orientation)
        end
    end
end

element_aesthetics(::HistogramStatistic) = [:x]

function default_scales(stat::HistogramStatistic)
    if stat.orientation == :vertical
        return [Gadfly.Scale.y_continuous()]
    else
        return [Gadfly.Scale.x_continuous()]
    end
end

const histogram = HistogramStatistic

function apply_statistic(stat::HistogramStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    if stat.orientation == :horizontal
        var = :y
        othervar = :x
        minvar = :ymin
        maxvar = :ymax
        drawmaxvar = :xdrawmax
        labelvar = :x_label
    else
        var = :x
        othervar = :y
        minvar = :xmin
        maxvar = :xmax
        drawmaxvar = :ydrawmax
        labelvar = :y_label
    end

    Gadfly.assert_aesthetics_defined("HistogramStatistic", aes, var)

    values = getfield(aes, var)
    if stat.minbincount > stat.maxbincount
        error("Histogram minbincount > maxbincount")
    end

    if isempty(getfield(aes, var))
        setfield!(aes, minvar, Float64[1.0])
        setfield!(aes, maxvar, Float64[1.0])
        setfield!(aes, var, Float64[1.0])
        setfield!(aes, othervar, Float64[0.0])
        return
    end

    if haskey(scales, var) && isa(scales[var], Scale.DiscreteScale)
        isdiscrete = true
        x_min = minimum(values)
        x_max = maximum(values)
        d = x_max - x_min + 1
        bincounts = zeros(Int, d)
        for x in values
            bincounts[x - x_min + 1] += 1
        end
    else
        isdiscrete = false
        value_set = collect(Set(values[Bool[Gadfly.isconcrete(v) for v in values]]))
        sort!(value_set)

        if  length(value_set) / length(values) < 0.9
            d, bincounts, x_max = choose_bin_count_1d_discrete(
                        values, value_set, stat.minbincount, stat.maxbincount)
        else
            d, bincounts, x_max = choose_bin_count_1d(
                        values, stat.minbincount, stat.maxbincount)
        end
    end

    x_min = Gadfly.concrete_minimum(values)
    binwidth = isdiscrete ? 1 : (x_max - x_min) / d

    if aes.color === nothing
        setfield!(aes, othervar, Array(Float64, d))
        if isdiscrete
            setfield!(aes, var, collect(Int, 1:d))
            setfield!(aes, othervar, bincounts)
        else
            setfield!(aes, minvar, Array(isdiscrete ? Int : Float64, d))
            setfield!(aes, maxvar, Array(isdiscrete ? Int : Float64, d))
            setfield!(aes, var, Array(isdiscrete ? Int : Float64, d))
            for j in 1:d
                getfield(aes, minvar)[j] = x_min + (j - 1) * binwidth
                getfield(aes, maxvar)[j] = x_min + j * binwidth
                getfield(aes, var)[j] = x_min + (j - 0.5) * binwidth
                getfield(aes, othervar)[j] = bincounts[j]
            end
        end
    else
        groups = Dict()
        for (x, c) in zip(values, cycle(aes.color))
            if !Gadfly.isconcrete(x)
                continue
            end

            if !haskey(groups, c)
                groups[c] = Float64[x]
            else
                push!(groups[c], x)
            end
        end

        if isdiscrete
            setfield!(aes, var, Array(Int, d * length(groups)))
        else
            setfield!(aes, minvar, Array(isdiscrete ? Int : Float64, d * length(groups)))
            setfield!(aes, maxvar, Array(isdiscrete ? Int : Float64, d * length(groups)))
            setfield!(aes, var, Array(isdiscrete ? Int : Float64, d * length(groups)))
        end

        setfield!(aes, othervar, Array(Float64, d * length(groups)))
        colors = Array(RGB{Float32}, d * length(groups))

        x_min = Gadfly.concrete_minimum(values)
        stack_height = zeros(Int, d)
        for (i, (c, xs)) in enumerate(groups)
            fill!(bincounts, 0)
            for x in xs
                if !Gadfly.isconcrete(x)
                    continue
                end
                if isdiscrete
                    bincounts[int(x)] += 1
                else
                    bin = max(1, min(d, int(ceil((x - x_min) / binwidth))))
                    bincounts[bin] += 1
                end
            end
            stack_height += bincounts[1:d]

            if isdiscrete
                for j in 1:d
                    idx = (i-1)*d + j
                    getfield(aes, var)[idx] = j
                    getfield(aes, othervar)[idx] = bincounts[j]
                    colors[idx] = c
                end
            else
                for j in 1:d
                    idx = (i-1)*d + j
                    getfield(aes, minvar)[idx] = x_min + (j - 1) * binwidth
                    getfield(aes, maxvar)[idx] = x_min + j * binwidth
                    getfield(aes, var)[idx] = x_min + (j - 0.5) * binwidth
                    getfield(aes, othervar)[idx] = bincounts[j]
                    colors[idx] = c
                end
            end

        end

        drawmax = float64(maximum(stack_height))
        aes_drawmax = getfield(aes, drawmaxvar)
        if aes_drawmax === nothing || aes_drawmax < drawmax
            setfield!(aes, drawmaxvar, drawmax)
        end

        aes.color = PooledDataArray(colors)
    end

    if haskey(scales, othervar)
        data = Gadfly.Data()
        setfield!(data, othervar, getfield(aes, othervar))
        Scale.apply_scale(scales[othervar], [aes], data)
    else
        setfield!(aes, labelvar, Scale.identity_formatter)
    end
end

immutable DensityStatistic <: Gadfly.StatisticElement
    # Number of points sampled
    n::Int

    function DensityStatistic(n=300)
        new(n)
    end
end

const density = DensityStatistic

element_aesthetics(::DensityStatistic) = [:x, :y]

default_scales(::DensityStatistic) = [Gadfly.Scale.y_continuous()]

function apply_statistic(stat::DensityStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("DensityStatistic", aes, :x)

    if aes.color === nothing
        if !isa(aes.x[1], Real)
            error("Kernel density estimation only works on Real types.")
        end

        x_f64 = collect(Float64, aes.x)
        # When will stat.n ever be <= 1? Seems pointless
        # certainly its length will always be 1
        window = stat.n > 1 ? KernelDensity.default_bandwidth(x_f64) : 0.1
        f = KernelDensity.kde(x_f64, bandwidth=window, npoints=stat.n)
        aes.x = collect(Float64, f.x)
        aes.y = f.density
    else
        groups = Dict()
        for (x, c) in zip(aes.x, cycle(aes.color))
            if !haskey(groups, c)
                groups[c] = Float64[x]
            else
                push!(groups[c], x)
            end
        end

        colors = Array(RGB{Float32}, 0)
        aes.x = Array(Float64, 0)
        aes.y = Array(Float64, 0)
        for (c, xs) in groups
            window = stat.n > 1 ? KernelDensity.default_bandwidth(xs) : 0.1
            f = KernelDensity.kde(xs, bandwidth=window, npoints=stat.n)
            append!(aes.x, f.x)
            append!(aes.y, f.density)
            for _ in 1:length(f.x)
                push!(colors, c)
            end
        end
        aes.color = PooledDataArray(colors)
    end
    aes.y_label = Gadfly.Scale.identity_formatter
end

immutable Histogram2DStatistic <: Gadfly.StatisticElement
    xminbincount::Int
    xmaxbincount::Int
    yminbincount::Int
    ymaxbincount::Int

    function Histogram2DStatistic(; xbincount=nothing,
                                    xminbincount=3,
                                    xmaxbincount=150,
                                    ybincount=nothing,
                                    yminbincount=3,
                                    ymaxbincount=150)
        if xbincount != nothing
            xminbincount = xbincount
            xmaxbincount = xbincount
        end

        if ybincount != nothing
            yminbincount = ybincount
            ymaxbincount = ybincount
        end

        new(xminbincount, xmaxbincount, yminbincount, ymaxbincount)
    end
end

element_aesthetics(::Histogram2DStatistic) = [:x, :y, :color]

default_scales(::Histogram2DStatistic) = [Gadfly.Scale.continuous_color()]

const histogram2d = Histogram2DStatistic

function apply_statistic(stat::Histogram2DStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)

    Gadfly.assert_aesthetics_defined("Histogram2DStatistic", aes, :x, :y)

    x_min, x_max = Gadfly.concrete_minimum(aes.x), Gadfly.concrete_maximum(aes.x)
    y_min, y_max = Gadfly.concrete_minimum(aes.y), Gadfly.concrete_maximum(aes.y)

    if haskey(scales, :x) && isa(scales[:x], Scale.DiscreteScale)
        x_categorial = true
        xminbincount = x_max - x_min + 1
        xmaxbincount = xminbincount
    else
        x_categorial = false
        xminbincount = stat.xminbincount
        xmaxbincount = stat.xmaxbincount
    end

    if haskey(scales, :y) && isa(scales[:y], Scale.DiscreteScale)
        y_categorial = true
        yminbincount = y_max - y_min + 1
        ymaxbincount = yminbincount
    else
        y_categorial = false
        yminbincount = stat.yminbincount
        ymaxbincount = stat.ymaxbincount
    end

    dy, dx, bincounts = choose_bin_count_2d(aes.x, aes.y,
                                            xminbincount, xmaxbincount,
                                            yminbincount, ymaxbincount)

    wx = x_categorial ? 1 : (x_max - x_min) / dx
    wy = y_categorial ? 1 : (y_max - y_min) / dx

    n = 0
    for cnt in bincounts
        if cnt > 0
            n += 1
        end
    end

    if x_categorial
        aes.x = Array(Int64, n)
    else
        aes.xmin = Array(Float64, n)
        aes.xmax = Array(Float64, n)
    end

    if y_categorial
        aes.y = Array(Int64, n)
    else
        aes.ymin = Array(Float64, n)
        aes.ymax = Array(Float64, n)
    end

    k = 1
    for i in 1:dy, j in 1:dx
        cnt = bincounts[i, j]
        if cnt > 0
            if x_categorial
                aes.x[k] = x_min + (j - 1)
            else
                aes.xmin[k] = x_min + (j - 1) * wx
                aes.xmax[k] = x_min + j * wx
            end

            if y_categorial
                aes.y[k] = y_min + (i - 1)
            else
                aes.ymin[k] = y_min + (i - 1) * wy
                aes.ymax[k] = y_min + i * wy
            end
            k += 1
        end
    end
    @assert k - 1 == n

    if !haskey(scales, :color)
        error("Histogram2DStatistic requires a color scale.")
    end
    color_scale = scales[:color]
    if !(typeof(color_scale) <: Scale.ContinuousColorScale)
        error("Histogram2DStatistic requires a continuous color scale.")
    end

    aes.color_key_title = "Count"

    data = Gadfly.Data()
    data.color = Array(Int, n)
    k = 1
    for cnt in transpose(bincounts)
        if cnt > 0
            data.color[k] = cnt
            k += 1
        end
    end

    if x_categorial
        aes.x = PooledDataArray(aes.x)
    end

    if y_categorial
        aes.y = PooledDataArray(aes.y)
    end

    Scale.apply_scale(color_scale, [aes], data)
    nothing
end

# Find reasonable places to put tick marks and grid lines.
immutable TickStatistic <: Gadfly.StatisticElement
    in_vars::Vector{Symbol}
    out_var::String

    # fixed ticks, or nothing
    ticks::Union(Nothing, AbstractArray)
end

function xticks(ticks::Union(Nothing, AbstractArray)=nothing)
    TickStatistic([:x, :xmin, :xmax, :xdrawmin, :xdrawmax], "x", ticks)
end

function yticks(ticks::Union(Nothing, AbstractArray)=nothing)
    TickStatistic(
        [:y, :ymin, :ymax, :middle, :lower_hinge, :upper_hinge,
         :lower_fence, :upper_fence, :ydrawmin, :ydrawmax], "y", ticks)
end

# Apply a tick statistic.
#
# Args:
#   stat: statistic.
#   aes: aesthetics.
#
# Returns:
#   nothing
#
# Modifies:
#   aes
#
function apply_statistic(stat::TickStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)

    if isa(coord, Coord.SubplotGrid)
        error("TickStatistic cannot be applied to subplot coordinates.")
    end

    in_group_var = symbol(string(stat.out_var, "group"))
    minval, maxval = nothing, nothing
    in_values = {}
    categorical = (:x in stat.in_vars && Scale.iscategorical(scales, :x)) ||
                  (:y in stat.in_vars && Scale.iscategorical(scales, :y))

    for var in stat.in_vars
        vals = getfield(aes, var)
        if vals != nothing && eltype(vals) != Function
            if minval == nothing
                minval = first(vals)
            end
            if maxval == nothing
                maxval = first(vals)
            end
            T = promote_type(typeof(minval), typeof(maxval))
            T = promote_type(T, eltype(vals))
            minval = convert(T, minval)
            maxval = convert(T, maxval)

            if stat.out_var == "x"
                dsize = aes.xsize === nothing ? [nothing] : aes.xsize
            elseif stat.out_var == "y"
                dsize = aes.ysize === nothing ? [nothing] : aes.ysize
            else
                dsize = [nothing]
            end

            size = aes.size === nothing ? [nothing] : aes.size

            minval, maxval = apply_statistic_typed(minval, maxval, vals, size, dsize)
            push!(in_values, vals)
        end
    end

    if isempty(in_values)
        return
    end

    in_values = chain(in_values...)

    # consider forced tick marks
    if stat.ticks != nothing
        minval = min(minval, minimum(stat.ticks))
        maxval = max(maxval, maximum(stat.ticks))
    end

    # TODO: handle the outliers aesthetic

    n = Gadfly.concrete_length(in_values)

    # take into account a forced viewport in cartesian coordinates.
    if typeof(coord) == Coord.Cartesian
        if stat.out_var == "x"
            if !is(coord.xmin, nothing)
                minval = min(minval, coord.xmin)
            end
            if !is(coord.xmax, nothing)
                maxval = max(maxval, coord.xmax)
            end
        elseif stat.out_var == "y"
            if !is(coord.ymin, nothing)
                minval = min(minval, coord.ymin)
            end
            if !is(coord.ymax, nothing)
                maxval = max(maxval, coord.ymax)
            end
        end
    end

    # check the x/yviewmin/max pesudo-aesthetics
    if stat.out_var == "x"
        if aes.xviewmin != nothing
            minval = aes.xviewmin
        end
        if aes.xviewmax != nothing
            maxval = aes.xviewmax
        end
    elseif stat.out_var == "y"
        if aes.yviewmin != nothing
            minval = aes.yviewmin
        end
        if aes.yviewmax != nothing
            maxval = aes.yviewmax
        end
    end

    # all the input values in order.
    if stat.ticks != nothing
        grids = ticks = stat.ticks
        viewmin = minval
        viewmax = maxval
        tickvisible = fill(true, length(ticks))
        tickscale = fill(1.0, length(ticks))
    elseif categorical
        ticks = Set()
        for val in in_values
            push!(ticks, val)
        end
        ticks = Float64[t for t in ticks]
        sort!(ticks)
        grids = (ticks .- 0.5)[2:end]
        viewmin = minimum(ticks)
        viewmax = maximum(ticks)
        tickvisible = fill(true, length(ticks))
        tickscale = fill(1.0, length(ticks))
    else
        minval, maxval = promote(minval, maxval)

        ticks, viewmin, viewmax =
            Gadfly.optimize_ticks(minval, maxval, extend_ticks=true)
        grids = ticks
        multiticks = Gadfly.multilevel_ticks(viewmin - (viewmax - viewmin),
                                             viewmax + (viewmax - viewmin))
        tickcount = length(ticks) + sum([length(ts) for ts in values(multiticks)])
        tickvisible = Array(Bool, tickcount)
        tickscale = Array(Float64, tickcount)
        i = 1
        for t in ticks
            tickscale[i] = 1.0
            tickvisible[i] = viewmin <= t <= viewmax
            i += 1
        end

        for (scale, ts) in multiticks
            for t in ts
                push!(ticks, t)
                tickvisible[i] = false
                tickscale[i] = scale
                i += 1
            end
        end
    end

    # We use the first label function we find for any of the aesthetics. I'm not
    # positive this is the right thing to do, or would would be.
    labeler = getfield(aes, symbol(string(stat.out_var, "_label")))

    setfield!(aes, symbol(string(stat.out_var, "tick")), ticks)
    setfield!(aes, symbol(string(stat.out_var, "grid")), grids)
    setfield!(aes, symbol(string(stat.out_var, "tick_label")), labeler)
    setfield!(aes, symbol(string(stat.out_var, "tickvisible")), tickvisible)
    setfield!(aes, symbol(string(stat.out_var, "tickscale")), tickscale)

    viewmin_var = symbol(string(stat.out_var, "viewmin"))
    if getfield(aes, viewmin_var) === nothing ||
       getfield(aes, viewmin_var) > viewmin
        setfield!(aes, viewmin_var, viewmin)
    end

    viewmax_var = symbol(string(stat.out_var, "viewmax"))
    if getfield(aes, viewmax_var) === nothing ||
       getfield(aes, viewmax_var) < viewmax
        setfield!(aes, viewmax_var, viewmax)
    end

    nothing
end

function apply_statistic_typed{T}(minval::T, maxval::T, vals, size, dsize)
#     for (val, s, ds) in zip(vals, cycle(size), cycle(dsize))
    lensize  = length(size)
    lendsize = length(dsize)
    for (i, val) in enumerate(vals)
        if !Gadfly.isconcrete(val) || !isfinite(val)
            continue
        end

        s = size[mod1(i, lensize)]
        ds = dsize[mod1(i, lendsize)]

        minval, maxval = minvalmaxval(minval, maxval, convert(T, val), s, ds)
    end
    minval, maxval
end

function apply_statistic_typed{T}(minval, maxval, vals::DataArray{T}, size, dsize)
    lensize  = length(size)
    lendsize = length(dsize)
    for i = 1:length(vals)
        if vals.na[i]
            continue
        end

        val::T = vals.data[i]
        s = size[mod1(i, lensize)]
        ds = dsize[mod1(i, lendsize)]

        minval, maxval = minvalmaxval(minval, maxval, val, s, ds)
    end
    minval, maxval
end

function minvalmaxval{T}(minval::T, maxval::T, val, s, ds)
    if val < minval || !isfinite(minval)
        minval = val
    end

    if val > maxval || !isfinite(maxval)
        maxval = val
    end

    if s != nothing
        minval = min(minval, val - s)::T
        maxval = max(maxval, val + s)::T
    end

    if ds != nothing
        minval = min(minval, val - ds)::T
        maxval = max(maxval, val + ds)::T
    end

    minval, maxval
end

immutable BoxplotStatistic <: Gadfly.StatisticElement
end

element_aesthetics(::BoxplotStatistic) = [:x, :y]

const boxplot = BoxplotStatistic

function apply_statistic(stat::BoxplotStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    if aes.y === nothing
        Gadfly.assert_aesthetics_defined("BoxplotStatistic", aes,
            :x, :lower_hinge, :upper_hinge, :lower_fence, :upper_fence)

        aes_color = aes.color === nothing ? [nothing] : aes.color
        groups = {}
        for (x, c) in zip(aes.x, cycle(aes_color))
            push!(groups, (x, c))
        end

        if !is(aes.color, nothing)
            aes.color = PooledDataArray(ColorValue[c for (x, c) in groups],
                                        levels(aes.color))
        end

        return
    end

    if aes.x === nothing
        aes_x = [1]
        aes.x_label = x -> fill("", length(x))
    else
        aes_x = aes.x
    end
    aes_color = aes.color === nothing ? [nothing] : aes.color

    T = isempty(aes.y) ? eltype(aes.y) : typeof(aes.y[1] / 1)
    groups = DefaultOrderedDict(() -> T[])

    for (x, y, c) in zip(cycle(aes_x), aes.y, cycle(aes_color))
        push!(groups[(x, c)], y)
    end

    if aes.y != nothing
        m = length(groups)
        aes.x = Array(eltype(aes.x), m)
        aes.middle = Array(T, m)
        aes.lower_hinge = Array(T, m)
        aes.upper_hinge = Array(T, m)
        aes.lower_fence = Array(T, m)
        aes.upper_fence = Array(T, m)
        aes.outliers = Vector{T}[]

        for (i, ((x, c), ys)) in enumerate(groups)
            sort!(ys)

            aes.x[i] = x
            aes.lower_hinge[i], aes.middle[i], aes.upper_hinge[i] =
                    quantile!(ys, [0.25, 0.5, 0.75])
            iqr = aes.upper_hinge[i] - aes.lower_hinge[i]

            idx = searchsortedfirst(ys, aes.lower_hinge[i] - 1.5iqr)
            aes.lower_fence[i] = ys[idx]

            idx = searchsortedlast(ys, aes.upper_hinge[i] + 1.5iqr)
            aes.upper_fence[i] = ys[idx]

            push!(aes.outliers,
                 filter(y -> y < aes.lower_fence[i] || y > aes.upper_fence[i], ys))
        end
    end

    if isa(aes_x, PooledDataArray)
        aes.x = PooledDataArray(aes.x, aes_x.pool)
    end

    if !is(aes.color, nothing)
        aes.color = PooledDataArray(RGB{Float32}[c for (x, c) in keys(groups)],
                                    levels(aes.color))
    end

    nothing
end

immutable SmoothStatistic <: Gadfly.StatisticElement
    method::Symbol
    smoothing::Float64

    function SmoothStatistic(; method::Symbol=:loess, smoothing::Float64=0.75)
        new(method, smoothing)
    end
end

const smooth = SmoothStatistic

element_aesthetics(::SmoothStatistic) = [:x, :y]

function apply_statistic(stat::SmoothStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)

    Gadfly.assert_aesthetics_defined("Stat.smooth", aes, :x, :y)
    Gadfly.assert_aesthetics_equal_length("Stat.smooth", aes, :x, :y)

    if !(stat.method in [:loess,:lm])
        error("The only Stat.smooth methods currently supported are loess and lm.")
    end

    num_steps = 750

    if aes.color === nothing
        x_min, x_max = minimum(aes.x), maximum(aes.x)

        if x_min == x_max
            error("Stat.smooth requires more than one distinct x value")
        end

        local xs, ys

        try
            xs = convert(Vector{Float64}, aes.x)
            ys = convert(Vector{Float64}, aes.y)
        catch e
            error("Stat.loess and Stat.lm require that x and y be bound to arrays of plain numbers.")
        end

        # loess can't predict points <x_min or >x_max. Make sure that doesn't
        # happen through a floating point fluke
        nudge = 1e-5 * (x_max - x_min)
        aes.x = collect((x_min + nudge):((x_max - x_min) / num_steps):(x_max - nudge))

        if stat.method == :loess
            aes.y = predict(loess(xs, ys, span=stat.smoothing), aes.x)
        elseif stat.method == :lm
            lmcoeff = linreg(xs,ys)
            aes.y = lmcoeff[2].*aes.x .+ lmcoeff[1]
        end
    else
        groups = Dict()
        aes_color = aes.color === nothing ? [nothing] : aes.color
        for (x, y, c) in zip(aes.x, aes.y, cycle(aes_color))
            if !haskey(groups, c)
                groups[c] = (Float64[], Float64[])
            end

            try
                push!(groups[c][1], x)
                push!(groups[c][2], y)
            catch
                error("Stat.loess and Stat.lm require that x and y be bound to arrays of plain numbers.")
            end
        end

        aes.x = Array(Float64, length(groups) * num_steps)
        aes.y = Array(Float64, length(groups) * num_steps)
        colors = Array(RGB{Float32}, length(groups) * num_steps)

        for (i, (c, (xs, ys))) in enumerate(groups)
            x_min, x_max = minimum(xs), maximum(xs)
            if x_min == x_max
                error("Stat.smooth requires more than one distinct x value")
            end
            nudge = 1e-5 * (x_max - x_min)
            steps = collect((x_min + nudge):((x_max - x_min) / num_steps):(x_max - nudge))

            if stat.method == :loess
                smoothys = predict(loess(xs, ys, span=stat.smoothing), steps)
            elseif stat.method == :lm
                lmcoeff = linreg(xs,ys)
                smoothys = lmcoeff[2].*steps .+ lmcoeff[1]
            end

            for (j, (x, y)) in enumerate(zip(steps, smoothys))
                aes.x[(i - 1) * num_steps + j] = x
                aes.y[(i - 1) * num_steps + j] = y
                colors[(i - 1) * num_steps + j] = c
            end
        end
        aes.color = PooledDataArray(colors)
    end
end

immutable HexBinStatistic <: Gadfly.StatisticElement
    xbincount::Int
    ybincount::Int

    function HexBinStatistic(; xbincount=50, ybincount=50)
        new(xbincount, ybincount)
    end
end

const hexbin = HexBinStatistic

function apply_statistic(stat::HexBinStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    xmin, xmax = minimum(aes.x), maximum(aes.x)
    ymin, ymax = minimum(aes.y), maximum(aes.y)
    xspan, yspan = xmax - xmin, ymax - ymin

    xsize = xspan / stat.xbincount
    ysize = yspan / stat.ybincount

    counts = Dict{(Any, Any), Int}()
    for (x, y) in zip(aes.x, aes.y)
        h = convert(HexagonOffsetOddR, pointhex(x - xmin + xspan/2,
                                                y - ymin + yspan/2,
                                                xsize, ysize))
        idx = (h.q, h.r)
        if !haskey(counts, idx)
            counts[idx] = 1
        else
            counts[idx] += 1
        end
    end

    N = length(counts)
    aes.x = Array(Float64, N)
    aes.y = Array(Float64, N)
    data = Gadfly.Data()
    data.color = Array(Int, N)
    k = 1
    for (idx, cnt) in counts
        x, y = center(HexagonOffsetOddR(idx[1], idx[2]), xsize, ysize,
                      xmin - xspan/2, ymin - yspan/2)
        aes.x[k] = x
        aes.y[k] = y
        data.color[k] = cnt
        k += 1
    end
    aes.xsize = [xsize]
    aes.ysize = [ysize]

    color_scale = scales[:color]
    if !(typeof(color_scale) <: Scale.ContinuousColorScale)
        error("HexBinGeometry requires a continuous color scale.")
    end

    Scale.apply_scale(color_scale, [aes], data)
end

function default_scales(::HexBinStatistic)
    return [Gadfly.Scale.continuous_color()]
end

immutable StepStatistic <: Gadfly.StatisticElement
    direction::Symbol

    function StepStatistic(; direction::Symbol=:hv)
        return new(direction)
    end
end

const step = StepStatistic

function element_aesthetics(::StepStatistic)
    return [:x, :y]
end

function apply_statistic(stat::StepStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("StepStatistic", aes, :x)
    Gadfly.assert_aesthetics_defined("StepStatistic", aes, :y)
    Gadfly.assert_aesthetics_equal_length("StepStatistic", aes, :x, :y)

    points = collect(zip(aes.x, aes.y))
    sort!(points, by=first)
    n = length(points)
    x_step = Array(eltype(aes.x), 2n - 1)
    y_step = Array(eltype(aes.y), 2n - 1)

    for i in 1:(2n-1)
        if isodd(i)
            x_step[i] = points[div(i-1,2)+1][1]
            y_step[i] = points[div(i-1,2)+1][2]
        elseif stat.direction == :hv
            x_step[i] = points[div(i-1,2)+2][1]
            y_step[i] = points[div(i-1,2)+1][2]
        else
            x_step[i] = points[div(i-1,2)+1][1]
            y_step[i] = points[div(i-1,2)+2][2]
        end
    end

    aes.x = x_step
    aes.y = y_step
end

immutable FunctionStatistic <: Gadfly.StatisticElement
    # Number of points to evaluate the function at
    num_samples::Int

    function FunctionStatistic(; num_samples=250)
        return new(num_samples)
    end
end

const func = FunctionStatistic

function default_scales(::FunctionStatistic)
    return [Gadfly.Scale.x_continuous(), Gadfly.Scale.y_continuous()]
end

function element_aesthetics(::FunctionStatistic)
    return [:y, :xmin, :xmax, :ymin, :ymax]
end

function apply_statistic(stat::FunctionStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("FunctionStatistic", aes, :y)
    Gadfly.assert_aesthetics_defined("FunctionStatistic", aes, :xmin)
    Gadfly.assert_aesthetics_defined("FunctionStatistic", aes, :xmax)
    Gadfly.assert_aesthetics_equal_length("FunctionStatistic", aes, :xmin, :xmax)

    aes.x = Array(Float64, length(aes.y) * stat.num_samples)
    ys = Array(Float64, length(aes.y) * stat.num_samples)

    i = 1
    for (f, xmin, xmax) in zip(aes.y, cycle(aes.xmin), cycle(aes.xmax))
        for x in linspace(xmin, xmax, stat.num_samples)
            aes.x[i] = x
            ys[i] = f(x)
            i += 1
        end
    end

    # color was bound explicitly
    if aes.color != nothing
        func_color = aes.color
        aes.color = DataArray(eltype(aes.color), length(aes.y) * stat.num_samples)
        groups = DataArray(Int, length(aes.y) * stat.num_samples)
        for i in 1:length(aes.y)
            aes.color[1+(i-1)*stat.num_samples:i*stat.num_samples] = func_color[i]
            groups[1+(i-1)*stat.num_samples:i*stat.num_samples] = i
        end
        aes.group = PooledDataArray(groups)
    elseif length(aes.y) > 1 && haskey(scales, :color)
        data = Gadfly.Data()
        data.color = Array(String, length(aes.y) * stat.num_samples)
        for i in 1:length(aes.y)
            fname = "f<sub>$(i)</sub>"
            data.color[1+(i-1)*stat.num_samples:i*stat.num_samples] = fname
        end
        Scale.apply_scale(scales[:color], [aes], data)
    end

    data = Gadfly.Data()
    data.y = ys
    Scale.apply_scale(scales[:y], [aes], data)
end

immutable ContourStatistic <: Gadfly.StatisticElement
    levels
    samples::Int

    function ContourStatistic(; n=15, samples=150)
        new(n, samples)
    end

    function ContourStatistic(; levels=15, samples=150)
        new(levels, samples)
    end
end

element_aesthetics(::ContourStatistic) = [:z, :xmin, :xmax, :ymin, :ymax]

const contour = ContourStatistic

function default_scales(::ContourStatistic)
    return [Gadfly.Scale.z_func(), Gadfly.Scale.x_continuous(),
            Gadfly.Scale.y_continuous(),
            Gadfly.Scale.continuous_color_gradient()]
end

function apply_statistic(stat::ContourStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    xs = aes.x === nothing ? nothing : convert(Vector{Float64}, aes.x)
    ys = aes.y === nothing ? nothing : convert(Vector{Float64}, aes.y)

    if typeof(aes.z) <: Function
        if xs == nothing && aes.xmin != nothing && aes.xmax != nothing
            xs = linspace(aes.xmin[1], aes.xmax[1], stat.samples)
        end

        if ys == nothing && aes.ymin != nothing && aes.ymax != nothing
            ys = linspace(aes.ymin[1], aes.ymax[1], stat.samples)
        end

        zs = Float64[aes.z(x, y) for x in xs, y in ys]

    elseif typeof(aes.z) <: Matrix
        zs = convert(Matrix{Float64}, aes.z)

        if xs == nothing
            xs = float([1:size(zs)[1]])
        end
        if ys == nothing
            ys = float([1:size(zs)[2]])
        end
        if size(zs) != (length(xs), length(ys))
            error("Stat.contour requires dimension of z to be length(x) by length(y)")
        end
    else
        error("Stat.contour requires either a matrix or a function")
    end

    levels = Float64[]
    contour_xs = eltype(xs)[]
    contour_ys = eltype(ys)[]

    groups = PooledDataArray(Int[])
    group = 0
    for contour in Contour.contours(xs, ys, zs, stat.levels)
        for curve in contour.lines
            for v in curve.vertices
                push!(contour_xs, v[1])
                push!(contour_ys, v[2])
                push!(levels, contour.level)
                push!(groups, group)
            end
            group += 1
        end
    end

    aes.group = groups
    color_scale = get(scales, :color, Gadfly.Scale.continuous_color_gradient())
    Scale.apply_scale(color_scale, [aes], Gadfly.Data(color=levels))
    Scale.apply_scale(scales[:x], [aes],  Gadfly.Data(x=contour_xs))
    Scale.apply_scale(scales[:y], [aes], Gadfly.Data(y=contour_ys))
end

immutable QQStatistic <: Gadfly.StatisticElement
end

element_aesthetics(::QQStatistic) = [:x, :y]

const qq = QQStatistic

function default_scales(::QQStatistic)
    return [Gadfly.Scale.x_continuous(), Gadfly.Scale.y_continuous]
end

function apply_statistic(stat::QQStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)

    Gadfly.assert_aesthetics_defined("Stat.qq", aes, :x, :y)

    # NOTES:
    #
    # apply_scales happens before apply_statistics, so we need to handle in
    # apply_scales the Distributions that might be bound to x and y... By
    # analogy with Stat.func, we can add a check in apply_statistic which defers
    # application.  Stat.func though requires an ARRAY of Functions, and doesn't
    # work on naked functions bound to aes.y.  If we want to bind Distributions,
    # we'd need to extend the types that are allowed for aes.y/.x (e.g. change
    # type of Aesthetics fields x and y).  Right now these are of type
    # NumericalOrCategoricalAesthetic.  The .x and .y fields are the _only_
    # place where this type is used, but I'm not sure if there's a reason that
    # changing this typealias would be a bad idea...for now I've just used a
    # direct `Union(NumericalOrCategoricalAesthetic, Distribution)`.
    #
    # TODO:
    #
    # Grouping by color etc.?

    # a little helper function to convert either numeric or distribution
    # variables to a format suitable to input to qqbuild.
    toVecOrDist = v -> typeof(v) <: Distribution ? v : convert(Vector{Float64}, v)

    # check and convert :x and :y to proper types for input to qqbuild
    local xs, ys
    try
        (xs, ys) = map(toVecOrDist, (aes.x, aes.y))
    catch e
        error("Stat.qq requires that x and y be bound to either a Distribution or to arrays of plain numbers.")
    end

    qqq = qqbuild(xs, ys)

    aes.x = qqq.qx
    aes.y = qqq.qy

    # apply_scale to Distribution-bound aesthetics is deferred, so re-apply here
    # (but only for Distribution, numeric data is already scaled).  Only one of
    # :x or :y can be a Distribution since qqbuild will throw an error for two
    # Distributions.
    data = Gadfly.Data()
    if typeof(xs) <: Distribution
        data.x = aes.x
        Scale.apply_scale(scales[:x], [aes], data)
    elseif typeof(ys) <: Distribution
        data.y = aes.y
        Scale.apply_scale(scales[:y], [aes], data)
    end

end

end # module Stat


text:
# Features are peaks the may occur in a spectrum. They include characteristic X-rays,
# Compton shifted characteristic X-rays, and escape peaks.

struct EscapeArtifact
    xray::CharXRay
    escape::CharXRay
    function EscapeArtifact(xray::CharXRay, escape::CharXRay = n"Si K-L3")
        @assert energy(xray) > energy(escape) "The energy($xray) must be larger than energy($escape)"
        return new(xray, escape)
    end
end

NeXLCore.exists(::Type{EscapeArtifact}, cxr::CharXRay, esc::CharXRay = n"Si K-L3") =
    energy(cxr) - energy(esc) > 100.0

NeXLCore.element(esc::EscapeArtifact) = element(esc.xray)
NeXLCore.brightest(escs::AbstractVector{EscapeArtifact}) =
    brightest(map(ea -> ea.xray, escs))

Base.show(io::IO, esc::EscapeArtifact) = print(io, name(esc))
Base.show(io::IO, escs::AbstractVector{EscapeArtifact}) = print(io, name(escs))
NeXLCore.name(esc::EscapeArtifact) = "Esc[$(esc.xray)]"
NeXLCore.name(escs::AbstractVector{EscapeArtifact}) =
    "Esc[$(join(map(esc->"$(esc.xray)", escs),","))]"

struct ComptonArtifact
    xray::CharXRay
    angle::Float64
    function ComptonArtifact(prim::CharXRay, θ::Float64)
        @assert (θ >= 0.0) && (θ <= π) "The Compton scatter angle must be between 0.0 and π"
        return new(prim, θ)
    end
end

Base.show(io::IO, ca::ComptonArtifact) = print(io, name(ca))
NeXLCore.name(ca::ComptonArtifact) = "Compton[$(ca.xray), $(rad2deg(ca.angle))°]"
NeXLCore.name(cas::AbstractVector{ComptonArtifact}) =
    "Compton[$(name(map(ca->ca.xray, cas)))]"
NeXLCore.element(ca::ComptonArtifact) = element(xray)

NeXLCore.energy(ca::ComptonArtifact) =
    energy(ca.xray) * NeXLCore.comptonShift(ca.angle, energy(ca.xray))
NeXLCore.weight(esc::ComptonArtifact) = weight(NormalizeToUnity, esc.xray)

NeXLCore.energy(esc::EscapeArtifact) = energy(esc.xray) - energy(esc.escape)

"""
    weight(esc::EscapeArtifact, factor=0.01)

The weight of an EscapeArtifact which is factor * weight(esc.xray).
"""
NeXLCore.weight(esc::EscapeArtifact, factor = 0.01) = factor * weight(NormalizeToUnity, esc.xray)

"""
    SpectrumFeature

A union representing the different type of peak-like features (helpful and harmful) that can appear in a spectrum.
"""
SpectrumFeature = Union{CharXRay,EscapeArtifact,ComptonArtifact}

Base.show(io::IO, sfs::AbstractVector{<:SpectrumFeature}) =
    map(sf -> println(io, name(sf)), sfs)

Base.convert(
    ::AbstractVector{SpectrumFeature},
    x::AbstractVector{CharXRay},
)::AbstractVector{SpectrumFeature} = x

charFeature(
    elm::Element,
    tr::Transition...;
    minweight = 1.0e-3,
    maxE = 1.0e6,
)::AbstractVector{SpectrumFeature} = characteristic(elm, tr, minweight, maxE)

escapeFeature(
    elm::Element,
    trs::Transition...;
    minweight = 0.1,
    maxE = 1.0e6,
    escape = n"Si K-L3",
)::AbstractVector{SpectrumFeature} =
    map(tr -> EscapeArtifact(tr, escape), characteristic(elm, trs, minweight, maxE))

comptonFeature(
    elm::Element,
    trs::Transition...;
    minweight = 1.0e-3,
)::AbstractVector{SpectrumFeature} =
    map(tr -> ComptonArtifact(tr, escape), characteristic(elm, trs, minweight, maxE))

ext: jl

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
using DomainSets: element, elements, numelements, ×

"""
Recombine the outcome of several invocations of `I,E = quadrature(...)` into
a single `I,E` tuple, by summing the integral values and error estimates.
"""
recombine_outcome(IEs) = sum(IE[1] for IE in IEs), sum(IE[2] for IE in IEs)

## Singularities

"""
Split the domain into a list of domains, such that any singularity lies
on the boundary of one or more of the constituting parts.

Care is taken to avoid a singularity in the interior due to roundoff errors.
"""
splitdomain_sing(sing, domain) = (domain,)

splitdomain_sing(sing::PointSingularity, domain) = splitdomain_point(point(sing), domain)

function splitdomain_point(point, domain::MappedDomain)
    m = inverse_map(domain)
    map_domain.(Ref(m), splitdomain_point(m(point), superdomain(domain)))
end

function splitdomain_point(point, domain::Domain)
    if point ∈ domain
        @warn("Don't know how to split domain $(domain) at point $(point)")
        (domain,)
    else
        (domain,)
    end
end

"Split the domain according to a singularity at the point `x`."
function splitdomain_point(x, domain::AbstractInterval)
    T = promote_type(typeof(x),prectype(domain))
    tol = tolerance(T)
    a, b = extrema(domain)
    if (x > a+tol) && (x < b-tol)
        # x is in [a,b] and sufficiently far away from the boundary
        (a..x-tol, x+tol..b)
    elseif (x < a-tol) || (x > b+tol)
        # x is not in [a,b] and sufficiently far away from the boundary
        (domain,)
    else
        # x is close to a or b
        if abs(x-a) <= tol
            (max(a,x)+tol..b,)
        else
            (a..min(b,x)-tol,)
        end
    end
end

splitdomain_point(x, domain::ProductDomain) = splitdomain_point(x, domain, elements(domain)...)

function splitdomain_point(x, domain::ProductDomain, domain1, domain2)
    domains1 = splitdomain_point(x[1], domain1)
    domains2 = splitdomain_point(x[2], domain2)
    [d1 × d2 for d1 in domains1,d2 in domains2]
end

splitdomain_sing(sing::DiagonalSingularity, domain) =
    splitdomain_diagonal(domain)

splitdomain_diagonal(domain::ProductDomain) = splitdomain_diagonal(domain, elements(domain)...)

function splitdomain_diagonal(domain::ProductDomain, domain1::AbstractInterval, domain2::AbstractInterval)
    T = prectype(domain)
    diff1 = domain1 \ domain2
    diff2 = domain2 \ domain1
    overlap = domain1 ∩ domain2
    if width(overlap) < 100eps(T)
        # "if isempty(overlap)" is not sufficient, because the overlap may be numerically small
        [domain1 × domain2]
    else
        a, b = extrema(overlap)
        if isempty(diff1) && isempty(diff2)
            [LowerRightTriangle(a, b), UpperRightTriangle(a, b)]
        elseif isempty(diff1)
            [LowerRightTriangle(a, b), UpperRightTriangle(a, b), domain1 × diff2]
        elseif isempty(diff2)
            [diff1 × domain2, LowerRightTriangle(a, b), UpperRightTriangle(a, b)]
        else
            [diff1 × domain2, LowerRightTriangle(a, b), UpperRightTriangle(a, b), overlap × diff2]
        end
    end
end

const SplittingSingularity = Union{PointSingularity,DiagonalSingularity}

function quadrature_s(qs, integrand, domain, measure, sing::SplittingSingularity)
    domains = splitdomain_sing(sing, domain)
    if length(domains) > 1
        IEs = quadrature_m.(Ref(qs), Ref(integrand), domains, Ref(measure), Ref(sing))
        recombine_outcome(IEs)
    else
        quadrature_m(qs, integrand, domain, measure, sing)
    end
end

## Measures

function quadrature_m(qs, integrand, domain, δ::DiracMeasure, sing)
    x = point(δ)
    if x ∈ domain
        I = integrand(x)
        E = zero(I)
    else
        I, E = zero_result(integrand, eltype(domain))
    end
    I, E
end

function quadrature_m(qs, integrand, domain, μ::DiscreteMeasure, sing)
    x = points(μ)
    w = weights(μ)
    I,E = zero_result(integrand, promote_type(prectype(domain),eltype(w)))
    for i in 1:length(x)
        if x[i] ∈ domain
            I += w[i]*integrand(x[i])
        end
    end
    I, E
end

# Lebesgue measures can pass through unaltered
# Other types of measures we "process" first, and then we continue with quadrature_d
quadrature_m(qs, integrand, domain, measure::AbstractLebesgueMeasure, sing) =
    quadrature_d(qs, integrand, domain, measure, sing)

function quadrature_m(qs, integrand, domain, measure::Measure{T}, sing) where {T}
    prefactor, map, domain2, measure2, sing2 = process_measure(qs, domain, measure, sing)
    integrand2 = transform_integrand(integrand, prefactor, map)
    quadrature_d(qs, integrand2, domain2, measure2, sing2)
end

struct Identity end
(::Identity)(x) = x
const id = Identity()

transform_integrand(integrand, prefactor::Identity, map::Identity) = integrand
transform_integrand(integrand, prefactor::Identity, map) = t -> integrand(map(t))
transform_integrand(integrand, prefactor, map::Identity) = t -> prefactor(t) * integrand(t)
transform_integrand(integrand, prefactor, map) = t -> prefactor(t) * integrand(map(t))

process_measure(qs, domain, measure::Measure, sing) =
    process_measure_default(qs, domain, measure, sing)

# By default we replace all measures by the LebesgueMeasure on the space
function process_measure_default(qs, domain, measure::Measure{T}, sing) where {T}
    prefactor = t -> unsafe_weight(measure, t)
    prefactor, id, domain, LebesgueMeasure{T}(), sing
end

# Truncate an infinite domain to a finite one for numerical evaluation of Hermite integrals
function process_measure(qs::AdaptiveStrategy, domain::FullSpace{T}, measure::HermiteMeasure{T}, sing) where {T}
    U = sqrt(-log(eps(T)))
    hermite_weight, id, -U..U, LebesgueMeasure{T}(), sing
end

# apply the cosine map for integrals with the ChebyshevT weight, to avoid the singularities
function process_measure(qs::AdaptiveStrategy, domain::ChebyshevInterval, measure::ChebyshevTMeasure{T}, sing) where {T}
    # Transformation is: f(t) -> pi*f(cos(pi*t))
    Tpi = convert(T, pi)
    prefactor = t -> Tpi
    map = t -> cos(Tpi*t)
    prefactor, map, UnitInterval{T}(), UnitLebesgueMeasure{T}(), sing
end

# same as above, but on a subinterval
function process_measure(qs::AdaptiveStrategy, domain::AbstractInterval, measure::ChebyshevTMeasure{T}, sing) where {T}
    Tpi = convert(T, pi)
    prefactor = t -> Tpi
    map = t -> cos(Tpi*t)
    a, b = extrema(domain)
    a < -1.001 && throw(BoundsError(measure, a))
    b > 1.001 && throw(BoundsError(measure, b))
    # Set a and b to be within [-1,1] in order to avoid errors with acos below
    a = max(a, -1)
    b = min(b, 1)
    prefactor, map, acos(b)/pi..acos(a)/pi, LebesgueMeasure{T}(), sing
end

# apply the cosine map for integrals with the ChebyshevU weight as well
function process_measure(qs::AdaptiveStrategy, domain::ChebyshevInterval, measure::ChebyshevUMeasure{T}, sing) where {T}
    Tpi = convert(T, pi)
    prefactor = t -> Tpi * sin(Tpi*t)^2
    map = t -> cos(Tpi*t)
    prefactor, map, UnitInterval{T}(), UnitLebesgueMeasure{T}(), sing
end

function process_measure(qs::AdaptiveStrategy, domain::AbstractInterval, measure::ChebyshevUMeasure{T}, sing) where {T}
    Tpi = convert(T, pi)
    prefactor = t -> Tpi * sin(Tpi*t)^2
    map = t -> cos(Tpi*t)
    a, b = extrema(domain)
    a < -1.001 && throw(BoundsError(measure, a))
    b > 1.001 && throw(BoundsError(measure, b))
    # Set a and b to be within [-1,1] in order to avoid errors with acos below
    a = max(a, -1)
    b = min(b, 1)
    prefactor, map, acos(b)/pi..acos(a)/pi, LebesgueMeasure{T}(), sing
end

# The "best" rule for certain measures becomes a Gauss rule
quadrature_m(qs::BestRule, integrand, domain::ChebyshevInterval{T}, μ::LegendreMeasure, sing) where {T} =
    quadrature_m(Q_GaussLegendre(T, qs.n), integrand, domain, μ, sing)

quadrature_m(qs::BestRule, integrand, domain::ChebyshevInterval, μ::JacobiMeasure{T}, sing) where {T} =
    quadrature_m(Q_GaussJacobi(qs.n, μ.α, μ.β), integrand, domain, LegendreMeasure{T}(), sing)

quadrature_m(qs::BestRule, integrand, domain::HalfLine, μ::LaguerreMeasure{T}, sing) where {T} =
    quadrature_m(Q_GaussLaguerre(qs.n, μ.α), integrand, domain, LebesgueMeasure{T}(), sing)

quadrature_m(qs::BestRule, integrand, domain::FullSpace{T}, μ::HermiteMeasure{T}, sing) where {T} =
    quadrature_m(Q_GaussHermite(qs.n), integrand, domain, LebesgueMeasure{T}(), sing)

## Domains

quadrature_d(qs, integrand, domain::EmptySpace, measure, sing) =
    zero_result(integrand, prectype(domain))

"Convert a union of domains into a vector of domains without overlap"
function nonoverlapping_domains(domain::UnionDomain{T}) where {T}
    domains = Vector{Domain{T}}()
    push!(domains, element(domain, 1))
    for i in 2:numelements(domain)
        d = element(domain, i)
        for j=1:i-1
            d = d \ domains[j]
        end
        push!(domains, d)
    end
    domains
end

expand_domain(domain::Domain) = (domain,)

expand_domain(domain::UnionDomain) = nonoverlapping_domains(domain)

expand_domain(domain::ProductDomain) = expand_domain(domain, elements(domain)...)

function expand_domain(domain::ProductDomain, domain1::Domain, domain2::Domain)
    domains1 = expand_domain(domain1)
    domains2 = expand_domain(domain2)
    [d1 × d2 for d1 in domains1,d2 in domains2]
end

function expand_domain(domain::ProductDomain, domain1::Domain, domain2::Domain, domain3::Domain)
    domains1 = expand_domain(domain1)
    domains2 = expand_domain(domain2)
    domains3 = expand_domain(domain3)
    [×(d1, d2, d3) for d1 in domains1,d2 in domains2,d3 in domains3]
end

function expand_domain(domain::ProductDomain, domain1::Domain, domain2::Domain, domain3::Domain, domain4::Domain)
    domains1 = expand_domain(domain1)
    domains2 = expand_domain(domain2)
    domains3 = expand_domain(domain3)
    domains3 = expand_domain(domain3)
    [×(d1,d2,d3,d4) for d1 in domains1,d2 in domains2,d3 in domains3]
end

const ExpandableDomain = Union{UnionDomain,ProductDomain}

function quadrature_d(qs, integrand, domain::MappedDomain, measure, sing)
    m = inverse_map(domain)
    quadrature_d(qs, t->integrand(m(t))/jacdet(m,t), superdomain(domain), measure, sing)
end

function quadrature_d(qs, integrand, domain::ExpandableDomain, measure, sing)
    domains = expand_domain(domain)
    if length(domains) > 1
        # We deliberately invoke quadrature here, starting the chain again from the start.
        # This gives the earlier steps a change to simplify again, with the new domains
        IEs = quadrature.(Ref(qs), Ref(integrand), domains, Ref(measure), Ref(sing))
        recombine_outcome(IEs)
    else
        # Nothing changed, we move on in the chain with the original domain
        select_quad(qs, integrand, domain, measure, sing)
    end
end

function quadrature_d(qs, integrand, domain::LowerRightTriangle{T}, measure::LebesgueMeasure, sing) where {T}
    # Restriction to Lebesgue measure because we would have to map the measure too
    # TODO with a more general framework for change-of-variables in integrals
    # For now, we hardcode the change of variables and we use Duffy's trick.
    a = domain.a
    b = domain.b
    # For x from a to b, and y from a to x, we set y = a+(x-a)*u, where u goes from 0 to 1. We then have dy = (x-1)du.
    d = UnitInterval{T}()
    square = (a..b) × d
    # TODO: change the singularity to something sensible: the diagonal is mapped to the right side of the square
    quadrature(qs, x -> integrand(SVector(x[1],a+x[2]*(x[1]-a))) * (x[1]-a),
        square, measure, NoSingularity())
end

function quadrature_d(qs, integrand, domain::UpperRightTriangle{T}, measure::LebesgueMeasure, sing) where {T}
    a = domain.a
    b = domain.b
    # For x from a to b, and y from x to b, we set y = b-(b-x)*u, where u goes from 0 to 1. We then have dy = (b-x)du.
    d = UnitInterval{T}()
    square = (a..b) × d
    quadrature(qs, x -> integrand(SVector(x[1], b-(b-x[1])*(1-x[2])))*(b-x[1]),
        square, measure, NoSingularity())
end


text:
"""
optimize
"""

export optimize!

function optimize!(::JuMP.Model, ::AbstractOptimizeAttribute) <: Bool end
function optimize!(::JuMP.Model) <: Bool end

# Will use all available methods for recognizion, and solving if recognized.
# Else use the provided solver if any is attatched.
# NB, cannot change model while solving.
function optimize!(model::JuMP.Model)

    lpmodel = get_lpmodel(model)

    if recognize(lpmodel, DifferenceConstraints()) && recognize(lpmodel, AllIntegerConstraintBounds())
        if solve!(lpmodel, model, ShortestPath())
            print("\nProgram solved with Bellman-Ford shortest path.\n")
            return true     # Replace with status
        end
    else
        print("\nDid not reccognize any faster solutions.\n")
    end

    if !(isnothing(model.moi_backend.optimizer))
        print("\nUsing the attatched optimizer...\n")
        JuMP.optimize!(model)
        return false
    else
        print("\nThe model cannot be solved without an attatched solver.\n")
        return false
    end

end

max_stars_repo_name: henriean/solver-peeker

max_issues_repo_name: henriean/solver-peeker

max_forks_repo_name: henriean/solver-peeker


text:
# include("../../datastructs/src/datastructs_module.jl")
# include("../../numder/src/numder_module.jl")
# include("../../nonlineq/src/nonlineq_module.jl")

include("../../../std/datastructs/src/datastructs_module.jl")
include("../../common/numder/src/numder_module.jl")

# update the ode functions such that, they are of form: ode_fcn!(dqdt, t, q) where dqdt is a commonly used array -> avoid allocations

module ode
using LinearAlgebra
using ..numder
# using ..nonlineq
using ..datastructs
using PyPlot
PyPlot.pygui(true)

# currently commented out for later refactoring
# include("fixed_step_size_solvers.jl")

include("utils.jl")
include("time_interval.jl")
include("ode_problem.jl")
include("basic_options.jl")
include("ode_result.jl")
include("ode_stat.jl")
include("error_estimation.jl")
include("stepper_options.jl")
include("solvers.jl")

function create_stepper_temporaries(ode_solver::S, n_prob_dim::Int) where {S <: AbstractAdaptiveExplicitSolver}
    stepper_temp = StepperTemporaries(
        n_prob_dim,
        n_intermediates(ode_solver),
        n_temporaries(ode_solver))

    return stepper_temp
end

function create_stepper_temporaries(ode_solver::S, n_prob_dim::Int) where {S <: AbstractRosenbrockSolver}
    stepper_temp = RosenbrockStepperTemporaries(
        n_prob_dim,
        n_intermediates(ode_solver),
        n_temporaries(ode_solver))

    return stepper_temp
end

function solve(
    ode_problem::OdeProblem,
    ode_solver::S;
    basic_options::BasicOdeOptions,
    stepper_options::StepperOptions = StepperOptions(),
    error_estimation_options::ErrorEstimationOptions = ErrorEstimationOptions()
    ) where {S <: AbstractOdeSolver}

    error("Default implementation")
end

function solve(ode_problem::OdeProblem)
    ode_solver = Tsitouras45()
    return solve(ode_problem, ode_solver)
end

function solve(
    ode_problem::OdeProblem,
    ode_solver::S;
    basic_options::BasicOdeOptions = BasicOdeOptions(n_dim(ode_problem)),
    stepper_options::StepperOptions = StepperOptions(),
    error_estimation_options::ErrorEstimationOptions = ErrorEstimationOptions(),
    step_size_controller::SSC = PIStepSizeControl()
    ) where {S <: Union{AbstractAdaptiveExplicitSolver, AbstractRosenbrockSolver}, SSC <: AbstractStepSizeControl}

    (ode_res, step_stat) = integrate_adaptive(
                                    ode_solver,
                                    ode_problem,
                                    basic_options;
                                    stepper_options          = stepper_options,
                                    error_estimation_options = error_estimation_options,
                                    step_size_controller     = step_size_controller)
    shrink_to_fit!(ode_res)

    return (ode_res, step_stat)
end

# add step size selection based on rejection/acceptance
# add event handling
# etc...

function do_nothing(x...)
    # [value, isterminal, direction] = my_event_fcn(t, y)
    return [nothing, nothing, nothing]
end

function integrate_adaptive(
    ode_solver::S,
    ode_problem::OdeProblem,
    basic_options::BasicOdeOptions;
    stepper_options::StepperOptions = StepperOptions(),
    error_estimation_options::ErrorEstimationOptions = ErrorEstimationOptions(),
    step_size_controller::SSC = PIStepSizeControl(),
    event_functions = [do_nothing]
    ) where {S <: Union{AbstractAdaptiveExplicitSolver, AbstractRosenbrockSolver}, SSC <: AbstractStepSizeControl}

    n_prob_dim   = n_dim(ode_problem)
    t_old        = ode_problem.time_interval.t_start
    x_old        = ode_problem.init_cond

    ode_res      = OdeResult(x_old)

    solver_order = order(ode_solver) # solver dependent
    step_stat    = OdeStat()

    der_x_old = similar(x_old)
    ode_problem.ode_fcn(der_x_old, t_old, x_old)

    stepper_temp = create_stepper_temporaries(ode_solver, n_prob_dim)
    set_t_old!( stepper_temp, t_old)
    set_x_old!( stepper_temp, x_old)
    stepper_temp.der_x[:,1] = der_x_old

    stepper_temp.dt = init_step_size(
                        ode_problem.ode_fcn,
                        der_x_old,
                        t_old,
                        x_old,
                        solver_order,
                        stepper_options,
                        basic_options,
                        error_estimation_options) # @time

    error_est = ErrorEstimation()

    for iter = 1:basic_options.max_iter

		stepper!(ode_solver, ode_problem.ode_fcn, basic_options, stepper_temp) # solver dependent # @time

		error_est.current = error_estimation(
                                        stepper_temp.x_new,
                                        stepper_temp.x_old,
                                        stepper_temp.x_est,
                                        basic_options,
                                        error_estimation_options) # @time

		if is_step_acceptable(error_est)
			rejection_break_flag = handle_accepted_step!(ode_res, error_est, stepper_temp, step_stat)
            # handle_event_function!(event_functions, ode_res, stepper_temp)
            # here comes the event function
            # here comes the output function: could be useful for controller state estimation
        else
			rejection_break_flag = handle_rejected_step!(step_stat, basic_options) # @time
		end

		if rejection_break_flag || stepper_temp.t_new >= ode_problem.time_interval.t_end
			break
		end

        step_size_control(
            step_size_controller,
            ode_solver,
            stepper_temp,
            error_est,
            rejection_break_flag,
            stepper_options)

	end

    return (ode_res, step_stat)
end

function handle_accepted_step!(ode_res, error_est, stepper_temp, step_stat)

	reset_rejected!(    step_stat )
    increase_accepted!( step_stat )

    update!(error_est)

    push!(ode_res, stepper_temp.t_new, stepper_temp.x_new)

    set_t_old!(stepper_temp, stepper_temp.t_new)
    set_x_old!(stepper_temp, stepper_temp.x_new)
    swap_fsal!(stepper_temp)

    rejection_break_flag = false
    return rejection_break_flag
end

function handle_rejected_step!(step_stat, basic_options)

	increase_rejected!(step_stat)

    rejection_break_flag = false

	if step_stat.current_rejected_steps >= basic_options.max_rejection_streak
        rejection_break_flag = true
		println("Too many rejected steps")
	end

    return rejection_break_flag
end

function init_step_size(
    ode_fcn,
    der_x0,
    t0,
    x0,
    solver_order,
    stepper_options,
    basic_options,
    error_estimation_options)

    # allocations:
    # x1
    # der_x0
    # der_x1
    # der_diff

    der_x1 = similar(x0)

    # base case
    if stepper_options.init_step_size > 0.0
        h_init = stepper_options.init_step_size
        return h_init
    end

    x_est = datastructs.Zero{Float64}()

    d0 = error_estimation(x0,     x0,     x_est, basic_options, error_estimation_options)
    d1 = error_estimation(der_x0, der_x0, x_est, basic_options, error_estimation_options)

    # some black magic from the Hairer ode solving book
	if d0 <= 1e-5 || d1 <= 1e-5
		h0 = 1e-6;
	else
		h0 = 0.01 * d0 / d1;
	end

    # compute one step of explicit euler
	x1 = x0 .+ h0 .* der_x0;

	# approx the derivative norm
	t1     = t0 + h0;
	ode_fcn(der_x1, t1, x1);

	der_diff = der_x1 - der_x0;
	d2 = 1 / h0 * error_estimation(der_diff, der_diff, x_est, basic_options, error_estimation_options)

	if max(abs(d1), abs(d2)) <= 1e-15
		h1 = max(1e-6, h0 * 1e-3);
	else
		h1 = 1e-2 / max(abs(d1), abs(d2));
		h1 = h1^(1/(solver_order + 1));
	end

	h_init = min(100 * h0, h1);

    return h_init
end

end # ode

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
@doc raw"""
`AbstractNTDist` is an abstract type for a wrapper type around measure(s)
and distributions(s).
The main idea is that instead of `rand`, `mean` and other statistical tools
wrapped objects return `NamedTuple`s or `TupleVector` when having a collection
of them.

The following API has to be implemented:
Given `π::AbstractNTDist` and `Π::AbstractVector{<:AbstractNTDist}`

## Necessary
- `ntrand(rng, π)` -> NamedTuple
- `ntmean(π)` -> NamedTuple
- `MeasureBase.logdensity_def(π, x::NamedTuple)` -> Real
- `Distributions.kldivergence(π₀, π₁)` -> Real
## Optional
- `tvrand(rng, Π)` -> TupleVector
- `tvmean(Π)` -> TupleVector
"""
abstract type AbstractNTDist end

# We replace the `rand` approach by `ntrand`
Base.rand(rng::AbstractRNG, ::Type{T}, π::AbstractNTDist) where {T} = ntrand(rng, π)
Base.rand(rng::AbstractRNG, π::AbstractNTDist) = ntrand(rng, π)
function Base.rand(rng::AbstractRNG, π::AbstractNTDist, n::Int)
    return TupleVector([ntrand(rng, π) for i in 1:n])
end

# Simple wrapper around any measure or distribution object
@doc raw"""
    NTDist(d) -> NTDist{typeof(d),:ω}
    NTDist{S}(d) -> NTDist{typeof(d),S}

Wrapper around a single distribution to be compatible with the [`ntrand`](@ref), [`ntmean`](@ref) interface.
One can pass the wanted symbol via `S` while the default will be `:ω`.
"""
struct NTDist{Td,S} <: AbstractNTDist
    d::Td
end

NTDist(d) = NTDist{typeof(d),:ω}(d)
NTDist{S}(d) where {S} = NTDist{typeof(d),S}(d)

dist(π::NTDist) = π.d

Statistics.mean(π::NTDist) = ntmean(π)
Statistics.mean(Π::AbstractVector{<:NTDist}) = tvmean(Π)

MeasureBase.logdensity_def(π::NTDist, x::NamedTuple) = logpdf(dist(π), only(x))

ntrand(rng::AbstractRNG, d) = ntrand(rng, NTDist(d))
function ntrand(rng::AbstractRNG, π::NTDist{Td,S}) where {Td,S}
    return NamedTuple{(S,)}((rand(rng, dist(π)),))
end

function tvrand(rng::AbstractRNG, Π::AbstractVector{<:T}) where {S,Td,T<:NTDist{Td,S}}
    return TupleVector(NamedTuple{(S,)}((rand.(rng, dist.(Π)),)))
end

ntmean(d) = ntmean(NTDist(d))
ntmean(π::NTDist{Td,S}) where {Td,S} = NamedTuple{(S,)}((mean(dist(π)),))
function tvmean(Π::AbstractVector{<:T}) where {S,Td,T<:NTDist{Td,S}}
    return TupleVector(NamedTuple{(S,)}((mean.(dist.(Π)),)))
end

Distributions.kldivergence(π₀::NTDist, π₁::NTDist) = kldivergence(dist(π₀), dist(π₁))


text:
# an idea: create one method named pivots with kwargs
@doc raw"""
    floorpivots(ohlc::TimeArray)

Floor Trader Pivots

# Formula

```math
\begin{align*}

  R_3            & = \text{Pivot}_t + (R_2 - S_1) \\
  R_2            & = \text{Pivot}_t + (R_1 - S_1) \\
  R_1            & = 2 \text{Pivot}_t - P^{low}_{t-1} \\
  \text{Pivot}_t & = P^{typical}_{t-1} =
      \frac{P^{high}_{t-1} + P^{low}_{t-1} + P^{close}_{t-1}}{3} \\
  S_1            & = 2 \text{Pivot}_t - P^{high}_{t-1} \\
  S_2            & = \text{Pivot}_t - (R_1 - S_1) \\
  S_3            & = \text{Pivot}_t - (R_2 - S_1)

\end{align*}
```
"""
function floorpivots(ohlc::TimeArray{T,N}) where {T,N}
  p  = lag(typical(ohlc))
  s1 = 2 .*p .- lag(ohlc[:High])
  r1 = 2 .*p .- lag(ohlc[:Low])
  s2 = p .- (r1 .- s1)
  r2 = (p .- s1) .+ r1
  s3 = p .- (r2 .- s1)
  r3 = (p .- s1) .+ r2

  TimeArray(
    timestamp(s1),
    [values(s3) values(s2) values(s1) values(p) values(r1) values(r2) values(r3)],
    [:s3, :s2, :s1, :pivot, :r1, :r2, :r3], meta(ohlc))
end

@doc raw"""
    woodiespivots(ohlc::TimeArray)

Woodie's Pivot

# Formula

```math
\begin{align*}
  \text{Range}   & = P^{high}_{t-1} - P^{low}_{t-1} \\
  R_4            & = S_4 + \text{Range} \tag{not implemented} \\
  R_3            & = R_1 + \text{Range} \\
  R_2            & = \text{Pivot}_t + \text{Range} \\
  R_1            & = 2 \text{Pivot}_t - P^{low}_{t-1} \\
  \text{Pivot}_t & =
      frac{P^{high}_{t-1} + P^{low}_{t-1} + 2 P^{open}_t}{4} \\
  S_1            & = 2 \text{Pivot}_t - P^{high}_{t-1} \\
  S_2            & = \text{Pivot}_t - \text{Range} \\
  S_3            & = S1 - \text{Range} \\
  S_4            & = S3 - \text{Range} \tag{not implemented}
\end{align*}
```
"""
function woodiespivots(ohlc::TimeArray{T,N}) where {T,N}
  rng = lag(ohlc[:High]) .- lag(ohlc[:Low])

  p  = (lag(ohlc[:High]) .+ lag(ohlc[:Low]) .+ 2 .*ohlc[:Open]) ./ 4
  s1 = 2 .*p .- lag(ohlc[:High])
  r1 = 2 .*p .- lag(ohlc[:Low])
  s2 = p .- rng
  r2 = p .+ rng
  s3 = s1 .- rng
  r3 = r1 .+ rng
  # s4 = s3 - range # can't get an answer that matches online calculators
  # r4 = r3 + range

  TimeArray(
    timestamp(s1),
    [values(s3) values(s2) values(s1) values(p) values(r1) values(r2) values(r3)],
    [:s3, :s2, :s1, :pivot, :r1, :r2, :r3], meta(ohlc))
end
#
# function camarillapivots{T,V}(hi::Array{SeriesPair{T,V},1},
#                               lo::Array{SeriesPair{T,V},1},
#                               cl::Array{SeriesPair{T,V},1})
#   #code here
#   ##  R4 = (H - L) * 1.1/2 + C
#   ##  R3 = (H - L) * 1.1/4 + C
#   ##  R2 = (H - L) * 1.1/6 + C
#   ##  R1 = (H - L) * 1.1/12 + C
#   ##  S1 = C - (H - L) * 1.1/12
#   ##  S2 = C - (H - L) * 1.1/6
#   ##  S3 = C - (H - L) * 1.1/4
#   ##  S4 = C - (H - L) * 1.1/2
# end
#
# function demarkpivots{T,V}(op::Array{SeriesPair{T,V},1},
#                            hi::Array{SeriesPair{T,V},1},
#                            lo::Array{SeriesPair{T,V},1},
#                            cl::Array{SeriesPair{T,V},1})
#   #code here
#   ## If Close < Open then X = (H + (L * 2) + C)
#   ## If Close > Open then X = ((H * 2) + L + C)
#   ## If Close = Open then X = (H + L + (C * 2))
#   ##
#   ## R1 = X / 2 - L
#   ## PP = X / 4
#   ## S1 = X / 2 - H
# end

function market_profile(x)
  #code here
end


text:
abstract type TimeSeriesFileFormat end
abstract type TimeSeriesFormatPeriodAsColumn <: TimeSeriesFileFormat end
abstract type TimeSeriesFormatYMDPeriodAsColumn <: TimeSeriesFormatPeriodAsColumn end
abstract type TimeSeriesFormatDateTimePeriodAsColumn <: TimeSeriesFormatPeriodAsColumn end
abstract type TimeSeriesFormatPeriodAsHeader <: TimeSeriesFileFormat end
abstract type TimeSeriesFormatYMDPeriodAsHeader <: TimeSeriesFormatPeriodAsHeader end
abstract type TimeSeriesFormatComponentsAsColumnsNoTime <: TimeSeriesFileFormat end
abstract type TimeSeriesFormatDateTimeAsColumn <: TimeSeriesFileFormat end

"""
Return a TimeArray from a CSV file.

Pass component_name when the file does not have the component name in a column header.
"""
function read_time_series(
    ::Type{T},
    data_file::AbstractString,
    component_name = nothing;
    kwargs...,
) where {T <: TimeSeriesData}
    if !isfile(data_file)
        msg = "TimeSeries file doesn't exist : $(data_file)"
        throw(DataFormatError(msg))
    end

    file = CSV.File(data_file)
    @debug "Read CSV data from $(data_file)." _group = LOG_GROUP_TIME_SERIES

    format = get_time_series_format(file)
    @debug "$format detected for the time series" _group = LOG_GROUP_TIME_SERIES
    return read_time_series(format, T, file, component_name; kwargs...)
end

function read_time_series(metadata::TimeSeriesFileMetadata; kwargs...)
    return read_time_series(
        metadata.time_series_type,
        metadata.data_file,
        metadata.component_name;
        kwargs...,
    )
end

"""
Return the time series format used in the CSV file.
"""
function get_time_series_format(file::CSV.File)
    columns = propertynames(file)
    has_ymd = :Year in columns && :Month in columns && :Day in columns
    has_period = :Period in columns
    has_datetime = :DateTime in columns

    if has_period
        if has_datetime
            format = TimeSeriesFormatDateTimePeriodAsColumn
        else
            format = TimeSeriesFormatYMDPeriodAsColumn
        end
    elseif has_ymd
        format = TimeSeriesFormatYMDPeriodAsHeader
    elseif has_datetime
        format = TimeSeriesFormatDateTimeAsColumn
    else
        format = TimeSeriesFormatComponentsAsColumnsNoTime
    end

    if format in (TimeSeriesFormatYMDPeriodAsColumn, TimeSeriesFormatYMDPeriodAsHeader)
        if !has_ymd
            throw(DataFormatError("$(file.name) is missing required Year/Month/Day"))
        end
    end

    return format
end

"""
Return the column names with values (components).
"""
function get_value_columns(::Type{TimeSeriesFormatYMDPeriodAsColumn}, file::CSV.File)
    return [x for x in propertynames(file) if !in(x, (:Year, :Month, :Day, :Period))]
end

function get_value_columns(
    ::Type{T},
    file::CSV.File,
) where {
    T <: Union{TimeSeriesFormatDateTimePeriodAsColumn, TimeSeriesFormatDateTimeAsColumn},
}
    return [x for x in propertynames(file) if !in(x, (:DateTime, :Period))]
end

"""
Return the column names with values.
"""
function get_value_columns(
    ::Type{TimeSeriesFormatComponentsAsColumnsNoTime},
    file::CSV.File,
)
    return propertynames(file)
end

"""
Return the column names that specify the Period.
"""
function get_period_columns(::Type{TimeSeriesFormatPeriodAsColumn}, file::CSV.File)
    return [:Period]
end

function get_period_columns(::Type{TimeSeriesFormatYMDPeriodAsHeader}, file::CSV.File)
    return [x for x in propertynames(file) if !in(x, (:Year, :Month, :Day))]
end

"""
Return a vector of dicts of unique timestamps and their counts.
"""
function get_unique_timestamps(::Type{T}, file::CSV.File) where {T <: TimeSeriesFileFormat}
    timestamps = Vector{Dict{String, Any}}()
    new_timestamp = x -> Dict("timestamp" => x, "count" => 1)

    for i in 1:length(file)
        timestamp = get_timestamp(T, file, i)
        if i == 1
            push!(timestamps, new_timestamp(timestamp))
        else
            if timestamp == timestamps[end]["timestamp"]
                timestamps[end]["count"] += 1
            else
                push!(timestamps, new_timestamp(timestamp))
            end
        end
    end

    @assert_op length(timestamps) > 0
    for timestamp in timestamps[2:end]
        @assert_op timestamp["count"] == timestamps[1]["count"]
    end

    return timestamps
end

"""
Return a Dates.DateTime for the row in the CSV file.
"""
function get_timestamp(
    ::Type{TimeSeriesFormatYMDPeriodAsColumn},
    file::CSV.File,
    row_index::Int,
)
    return Dates.DateTime(file.Year[row_index], file.Month[row_index], file.Day[row_index])
end

function get_timestamp(
    ::Type{T},
    file::CSV.File,
    row_index::Int,
) where {
    T <: Union{TimeSeriesFormatDateTimePeriodAsColumn, TimeSeriesFormatDateTimeAsColumn},
}
    return Dates.DateTime(file.DateTime[row_index])
end

function get_timestamp(
    ::Type{TimeSeriesFormatYMDPeriodAsHeader},
    file::CSV.File,
    row_index::Int,
)
    return get_timestamp(TimeSeriesFormatYMDPeriodAsColumn, file, row_index)
end

"""
Return a RawTimeSeries from a CSV file.

Pass component_name when the file does not have the component name in a column header.
"""
function read_time_series(
    ::Type{T},
    ::Type{Deterministic},
    file::CSV.File,
    component_name = nothing;
    kwargs...,
) where {T <: TimeSeriesFormatDateTimeAsColumn}
    @debug "Read CSV data from $file." _group = LOG_GROUP_TIME_SERIES
    horizon = length(first(file)) - 1
    data = SortedDict{Dates.DateTime, Vector{Float64}}()
    # First element in the row is the time series. We use integer indexes not to rely on
    # column names
    for row in file
        vector = Vector{Float64}(undef, horizon)
        for i in 1:horizon
            vector[i] = row[i + 1]
        end
        data[Dates.DateTime(row.DateTime)] = vector
    end
    return RawTimeSeries(first(keys(data)), data, horizon)
end

"""
Return a TimeSeries.TimeArray representing the CSV file.

This version of the function only has component_name to match the interface. It is unused.
"""
function read_time_series(
    ::Type{T},
    ::Type{<:StaticTimeSeries},
    file::CSV.File,
    component_name = nothing;
    kwargs...,
) where {T <: Union{TimeSeriesFormatPeriodAsColumn, TimeSeriesFormatDateTimeAsColumn}}
    first_timestamp = get_timestamp(T, file, 1)
    value_columns = get_value_columns(T, file)
    vals = [(string(x) => getproperty(file, x)) for x in value_columns]
    series_length = length(vals[1][2])
    return RawTimeSeries(first_timestamp, Dict(vals...), series_length)
end

"""
This version of the function supports the format where there is no column header for
a component, so the component_name must be passed in.
"""
function read_time_series(
    ::Type{T},
    ::Type{<:StaticTimeSeries},
    file::CSV.File,
    component_name::AbstractString;
    kwargs...,
) where {T <: TimeSeriesFormatPeriodAsHeader}
    period_cols_as_symbols = get_period_columns(T, file)
    period = [parse(Int, string(x)) for x in period_cols_as_symbols]
    first_timestamp = get_timestamp(T, file, 1)
    vals = Vector{Float64}()
    for i in 1:length(file)
        for period in period_cols_as_symbols
            val = getproperty(file, period)[i]
            push!(vals, val)
        end
    end

    return RawTimeSeries(first_timestamp, Dict(component_name => vals), length(vals))
end

"""
This version of the function only has component_name to match the interface.
It is unused.

Set start_datetime as a keyword argument for the starting timestamp, otherwise the current
day is used.
"""
function read_time_series(
    ::Type{T},
    ::Type{<:StaticTimeSeries},
    file::CSV.File,
    component_name = nothing;
    kwargs...,
) where {T <: TimeSeriesFormatComponentsAsColumnsNoTime}
    first_timestamp = get(kwargs, :start_datetime, Dates.DateTime(Dates.today()))
    value_columns = get_value_columns(T, file)
    vals = [(string(x) => getproperty(file, x)) for x in value_columns]
    series_length = length(vals[1][2])
    return RawTimeSeries(first_timestamp, Dict(vals...), series_length)
end

"""
Return the number of steps specified by the period in the file.
"""
function get_num_steps(
    ::Type{T},
    file::CSV.File,
    period::AbstractArray,
) where {T <: TimeSeriesFileFormat}
    error("Unsupported time series file format")
end

"""
Return the number of steps specified by the period in the file.
"""
function get_num_steps(
    ::Type{T},
    file::CSV.File,
    period::AbstractArray,
) where {T <: Union{TimeSeriesFormatPeriodAsColumn, TimeSeriesFormatDateTimeAsColumn}}
    timestamps = get_unique_timestamps(T, file)
    return timestamps[1]["count"]
end

"""
Return the number of steps specified by the period in the file.
"""
function get_num_steps(
    ::Type{T},
    file::CSV.File,
    period::AbstractArray,
) where {T <: TimeSeriesFormatPeriodAsHeader}
    return num_steps = period[end]
end


text:
import Distributions.MvNormal

struct BinChoice <: DiscreteUnivariateDistribution end
struct NoNoise <: DiscreteUnivariateDistribution end

rand(::BinChoice) = sample([-1, 1])
rand(::NoNoise) = 0.0
rand(::NoNoise, i::Int) = zeros(i)

const BenchmarkXDists = Union{Normal, MvNormal, Uniform}
const BenchmarkwDists = Union{BinChoice, Normal, Uniform}
const BenchmarkNoiseDists = Union{NoNoise, Normal}

abstract type XCorrelation end
struct NoCorrelation <: XCorrelation end
struct MatrixCorrelation <: XCorrelation
    coeff::Float64
end

struct XData{T <: BenchmarkXDists}
    dist::T
    corr::XCorrelation
end
XData() = XData(MvNormal, NoCorrelation())

"""
    BenchmarkData(; Xdata::XData = XData(),
                         wdist::BenchmarkwDists = BinChoice(),
                         noisedist::BenchmarkNoiseDists = NoNoise(),
                         SNR::Float64 = 20.0,
                         n::Vector{Int} = 100,
                         nfeatures::Int = 100,
                         k::Int = 10
                    )

Constructs a data object that fully describes the regression task as far as
this library is concerned.

# Arguments
* Xdata: how synthetic independent variables will be constructed (an object
of type `XData` has a distribution and a correlation mechanism)
* wdist: the distribution that regression coefficients will be sampled from
* noisedist: the distribution that noise will be sampled from
* SNR: signal-to-noise ration ||Y|| / ||noise|| we will use to normalize Y
* n: range of "number of observations" we want to test
* nfeatures: number of features
* sparsity: underlying sparsity pattern
"""
mutable struct BenchmarkData
    Xdata::XData
    wdist::BenchmarkwDists
    noisedist::BenchmarkNoiseDists
    SNR::Float64
    n::Vector{Int}
    nfeatures::Int
    sparsity::Int
end
function BenchmarkData(; Xdata::XData = XData(),
                         wdist::BenchmarkwDists = BinChoice(),
                         noisedist::BenchmarkNoiseDists = NoNoise(),
                         SNR::Float64 = 20.0,
                         n::Vector{Int} = [100],
                         nfeatures::Int = 100,
                         sparsity::Int = 10
                    )
    BenchmarkData(Xdata, wdist, noisedist, SNR, n, nfeatures, sparsity)
end

"""
    XData(::Type{T}, c::XCorrelation, d::Int) where {T <: BenchmarkXDists}

Instantiates a distribution of type `T` and modifies it,so that we follow
 correlation type `c`. Returns type `XData`, caching the distribution and
correlation type we used.
"""
function XData(::Type{T}, c::XCorrelation, d::Int) where {T <: BenchmarkXDists}
    error("Need to implement method XData for $T and $c.")
end
function XData(::Type{MvNormal}, c::NoCorrelation, d::Int)
    XData(MvNormal(zeros(d), eye(d)), c)
end
function XData(::Type{MvNormal}, c::MatrixCorrelation, d::Int)
    Σ = c.coeff * ones(d, d)
    # TODO smarter way to make correlated random variables. can't cache huge covariance matrices.
    @inbounds for i = 1:d
        Σ[i, i] = 1.0
    end
    XData(MvNormal(zeros(d), Σ), c)
end
function XData(::Type{Uniform}, c::NoCorrelation, ::Int)
    XData(Uniform(), c)
end

function getw(d::Int, k::Int, wdist::BenchmarkwDists)
    w = zeros(d)
    support = sample(1:d, k, replace=false)
    w[support] = rand(wdist, k)
    w
end
function getw(bd::BenchmarkData)
    getw(bd.nfeatures, bd.sparsity, bd.wdist)
end

function getX(n::Int, d::Int, Xdata::XData)
    rand(Xdata.dist, n, d)
end
function getX(n::Int, ::Int, Xdata::XData{T}) where {T <: MvNormal}
    rand(Xdata.dist, n)'
end
function getX(bd::BenchmarkData, i::Int)
    getX(bd.n[i], bd.nfeatures, bd.Xdata)
end

function getnoise(bnd::NoNoise, ::Float64, Y::Vector{Float64})
    n = length(Y)
    zeros(n)
end
function getnoise(bnd::BenchmarkNoiseDists, SNR::Float64, Y::Vector{Float64})
    n = length(Y)
    noise = rand(bnd, n)
    norm(noise) < 1e-3 && return noise
    (SNR < 1e-3) && error("Your SNR is too small.")
    noise .*= (norm(Y) / ( SNR * norm(noise) ) )
end
function getnoise(bd::BenchmarkData, Y::Vector{Float64})
    getnoise(bd.noisedist, bd.SNR, Y)
end


text:
##
using Plots
# import VortexDistributions:unwrap
using VortexDistributions

##
# function unwrap(phase::Array{Float64,1})
#     uphase = copy(phase)
#     s1 = length(phase)
#         @inbounds (k=uphase[1] - uphase[s1];k > π) && (uphase[1] -= 2π*div(k,2pi,RoundNearest))
#         @inbounds (k=uphase[1] - uphase[s1];k < -π) && (uphase[1] += 2π*div(-k,2pi,RoundNearest))
#         for i in 2:s1
#         @inbounds (k=uphase[i] - uphase[i-1];k > π) && (uphase[i] -= 2π*div(k,2pi,RoundNearest))
#         @inbounds (k=uphase[i] - uphase[i-1];k < -π) && (uphase[i] += 2π*div(-k,2pi,RoundNearest))
#         end
#         return uphase
# end

##
N=1000
θ = LinRange(0,4*pi,N)
psi = exp.(im*θ) + 0.1*(randn(N)+im*randn(N))
ϕ = angle.(psi)
ϕu = unwrap(ϕ)
ϕui = similar(ϕu)
unwrap!(ϕui,ϕ)

#
plot(ϕ,label="phi",c=:red,alpha=0.4,lw=6)
plot!(ϕu,label="phiu",lw=2,c=:black)
plot!(ϕui,label="phiu",lw=2,c=:green)


text:
module Corpuscles

using DocStringExtensions
using DelimitedFiles
using Unitful
using UnitfulAtomic
using Printf

import Base

export Particle, PDGID, PythiaID, Geant3ID, particles

# Julia 1.0 compatibility
eachrow_(x) = (x[i, :] for i in 1:size(x)[1])

const _data_dir = abspath(joinpath(@__DIR__, "..", "data"))

function Base.parse(::Type{Rational{T}}, val::AbstractString) where {T <: Integer}
    !('/' in val) && return parse(T, val) // 1
    nums, denoms = split(val, '/', keepempty=false)
    num = parse(T, nums)
    denom = parse(T, denoms)
    return num//denom
end

abstract type ParticleID end

struct PDGID <: ParticleID
    value::Int32
end

struct Geant3ID <: ParticleID
    value
end

@deprecate GeantID Geant3ID true

struct PythiaID <: ParticleID
    value
end

@enum PDGStatus begin
    Common      = 0
    Rare        = 1
    Unsure      = 2
    Further     = 3
    Nonexistent = 4
end

@enum InvProperty begin
    Same = 0
    Barred = 1
    ChargeInv = 2
end

struct MeasuredValue{D}
    value::Quantity{T1,D,U1} where {T1 <: Real, U1 <: Unitful.Units}
    lower_limit::Quantity{T2,D,U2} where {T2 <: Real, U2 <: Unitful.Units}
    upper_limit::Quantity{T3,D,U3} where {T3 <: Real, U3 <: Unitful.Units}
end

function Base.isless(x::MeasuredValue, y::Quantity)
    x.value + x.upper_limit < y
end

function Base.isless(x::Quantity, y::MeasuredValue)
    x < y.value - y.lower_limit
end

function Base.isless(x::MeasuredValue, y::MeasuredValue)
    x.value + x.upper_limit < y.value - y.lower_limit
end

function Base.isapprox(x::Quantity, y::MeasuredValue)
    (x < y.value + y.upper_limit) & (x > y.value - y.lower_limit)
end

function Base.isapprox(y::MeasuredValue, x::Quantity)
    (x < y.value + y.upper_limit) & (x > y.value - y.lower_limit)
end

function Base.isapprox(x::MeasuredValue, y::MeasuredValue)
    # TODO: I think this can be done way more elegant, but functionality is now given
    retval = isapprox(x, y.value + y.upper_limit)
    retval |= isapprox(x, y.value - y.lower_limit)
    retval |= isapprox(x, y.value)
    retval |= isapprox(y, x.value + x.upper_limit)
    retval |= isapprox(y, x.value - x.upper_limit)
    retval |= isapprox(y, x.value)
    retval
end

const _energy_dim = Unitful.dimension(u"J")
const _charge_dim = Unitful.dimension(u"C")

struct Particle
    pdgid::PDGID
    mass::MeasuredValue{_energy_dim}
    width::Union{Missing, MeasuredValue{_energy_dim}}
    charge::Quantity{T,_charge_dim,U} where {T<:Real, U<: Unitful.Units}
    isospin::Union{Missing, Rational{Int8}}
    parity::Union{Missing, Int8}
    gparity::Union{Missing, Int8}
    cparity::Union{Missing, Int8}
    antiprop::InvProperty
    rank::Int8
    status::PDGStatus
    name::String
    quarks::String
    latex::String
end

function read_conversion_csv(filepath::AbstractString)
    file_content = readdlm(filepath, ',', AbstractString, skipstart=2, comments=true)
    conversions = parse.(Int, file_content[:,1:3])
end

const _id_conversion_cols = Dict(PDGID => 1, Geant3ID => 3, PythiaID => 2)
const _id_conversion_tbl = read_conversion_csv(joinpath(_data_dir, "conversions.csv"))

Particle(id::ParticleID) = catalog.particle_dict[Base.convert(PDGID, id)]
Particle(id::Integer) = Particle(PDGID(id))

struct IDException <: Exception
    var::AbstractString
end

Base.showerror(io::IO, e::IDException) = Printf.@printf(io, "ParticleID Error: %s", e.var)

function Base.convert(t::Type{X}, id::Y) where {X <: ParticleID, Y <: ParticleID}
    if isequal(X, Y)
        return id
    end
    val_col = _id_conversion_cols[t]
    key_col = _id_conversion_cols[Y]
    row = findfirst(x->isequal(x, id.value), _id_conversion_tbl[:,key_col])
    if iszero(id.value) || isequal(row, nothing)
        throw(IDException("No corresponding $t for $id found!"))
    end
    X(_id_conversion_tbl[row, val_col])
end

function read_parity(val::AbstractString)
    tmp = parse(Int8, val)
    if tmp == 5
        return missing
    else
        return tmp
    end
end

const ParticleDict = Dict{PDGID, Particle}

function read_particle_csv(filepath::AbstractString)
    file_content = readdlm(filepath, ',', AbstractString, comments=true)
    header = string.(file_content[1,:])
    dct_particles = ParticleDict()
    for row in eachrow_(file_content[2:end,:])
        pdgid       = PDGID(parse(Int64, row[1]))
        mass_value  = parse(Float64, row[2]) * u"MeV"
        mass_lower  = parse(Float64, row[3]) * u"MeV"
        mass_upper  = parse(Float64, row[4]) * u"MeV"
        mass = MeasuredValue{_energy_dim}(mass_value, mass_lower, mass_upper)
        width_value  = parse(Float64, row[5]) * u"MeV"
        width_lower  = parse(Float64, row[6]) * u"MeV"
        width_upper  = parse(Float64, row[7]) * u"MeV"
        width = MeasuredValue{}(width_value, width_lower, width_upper)
        isospin = if (row[8] in ["", "?"]) missing else parse(Rational{Int8}, row[8]) end
        gparity = read_parity(row[9])
        parity = read_parity(row[10])
        cparity = read_parity(row[11])
        antiprop = InvProperty(parse(Int8, row[12]))
        charge = parse(Int8, row[13]) // 3 * u"e_au"
        rank = parse(Int8, row[14])
        status = PDGStatus(parse(Int8, row[15]))
        name = row[16]
        quarks = row[17]
        latex = row[18]
        dct_particles[pdgid] = Particle(pdgid,
                                        mass,
                                        width,
                                        charge,
                                        isospin,
                                        parity,
                                        gparity,
                                        cparity,
                                        antiprop,
                                        rank,
                                        status,
                                        name,
                                        quarks,
                                        latex)
    end
    dct_particles
end

"""
$(SIGNATURES)

Function to get the available catalog files which are available within
the package and returns a list with the absolute filepaths.

# Examples
```julia-repl
julia> Corpuscles.available_catalog_files()
["/home/foobar/dev/Corpuscles.jl/data/particle2019.csv"]
```
"""
function available_catalog_files()
    dir_content = readdir(_data_dir)
    filter!(s->occursin(".csv",s), dir_content)
    filter!(s->!occursin("conversions", s), dir_content)
    joinpath.(_data_dir, dir_content)
end

const _catalogs = available_catalog_files()

const _default_year = "2020"
const _default_catalog = filter(s->occursin(_default_year,s), _catalogs)[end]

mutable struct Catalog
    particle_dict::ParticleDict
    particles::Vector{Particle}

    Catalog(d::ParticleDict) = new(d, collect(values(d)))
end

"""
$(SIGNATURES)

Returns the full list of particles from the currently selected catalog.
"""
particles() = catalog.particles

const catalog = Catalog(read_particle_csv(_default_catalog))

"""
$(SIGNATURES)

This function reads a given catalog file and sets it as reference

# Arguments
- `filepath::AbstractString`: filepath to the catalog file

# Examples
```julia-repl
julia> Corpuscles.use_catalog_file("/home/foobar/dev/Corpuscles.jl/data/particle2019.csv")
```
"""
function use_catalog_file(filepath::AbstractString)
    # TODO: there is surely a better way to manage this "closure"-like thing
    catalog.particle_dict = read_particle_csv(filepath)
    catalog.particles = collect(values(catalog.particle_dict))
    return
end

function Base.show(io::IO, m::MeasuredValue)
    if isapprox(m.upper_limit, m.lower_limit)
        print(io, "$(m.value) ± $(m.lower_limit)")
    else
        print(io, "$(m.value) + $(m.upper_limit) - $(m.lower_limit)")
    end
    return
end

function Base.show(io::IO, p::Particle)
    Printf.@printf(io, "Particle(%s) '%s'", p.pdgid.value, p.name)
end

function Base.print(io::IO, p::Particle)
    Printf.@printf(io, "%-8s %s\n", "Name:", p.name)
    Printf.@printf(io, "%-8s %s\n", "PDG ID:", p.pdgid.value)
    Printf.@printf(io, "%-8s %s\n", "LaTeX:", "\$$(p.latex)\$")
    Printf.@printf(io, "%-8s %s\n", "Status:", p.status)
    fields = Dict("Mass" => p.mass,
                  "Width" => p.width,
                  "Q (charge)" => p.charge,
                  "C (charge parity)" => p.cparity,
                  "P (space parity)" => p.parity,
                  "G (G-parity)" => p.gparity,
                  "Isospin" => p.isospin,
                  "Composition" => p.quarks)
    for (key, value) in fields
        if value isa MeasuredValue || !ismissing(value) && !isempty(value)
            Printf.@printf(io, "%s = %s\n",key, value)
        end
    end
end

end # module


text:
export Sequence

struct Sequence <: CompositeBlock{Any}
    blocks::Vector
end

Sequence(args...) = Sequence(collect(AbstractBlock, args))

Yao.subblocks(seq::Sequence) = filter(x->x isa AbstractBlock, seq.blocks)
Yao.chsubblocks(pb::Sequence, blocks::Vector) = Sequence(blocks)

function Yao.apply!(reg::ArrayReg, seq::Sequence)
    for x in seq.blocks
        reg |> x
    end
    reg
end

for PROP in [:lastindex, :firstindex, :getindex, :length, :eltype, :iterate, :eachindex, :popfirst!, :pop!]
    @eval Base.$PROP(c::Sequence, args...; kwargs...) = $PROP(c.blocks, args...; kwargs...)
end

function Base.:(==)(lhs::Sequence, rhs::Sequence)
    (length(lhs.blocks) == length(rhs.blocks)) && all(lhs.blocks .== rhs.blocks)
end

Base.copy(c::Sequence) = Sequence(copy(c.blocks))
Base.similar(c::Sequence) = Sequence(empty!(similar(c.blocks)))
Base.getindex(c::Sequence, index::Union{UnitRange, Vector}) = Sequence(getindex(c.blocks, index))

Base.setindex!(c::Sequence, val::AbstractBlock{N}, index::Integer) where N = (setindex!(c.blocks, val, index); c)
Base.insert!(c::Sequence, index::Integer, val::AbstractBlock{N}) where N = (insert!(c.blocks, index, val); c)
Base.push!(c::Sequence, m) where N = (push!(c.blocks, m); c)
Base.append!(c::Sequence, list::Vector) where N = (append!(c.blocks, list); c)
Base.append!(c1::Sequence, c2::Sequence) where N = (append!(c1.blocks, c2.blocks); c1)
Base.prepend!(c1::Sequence, list::Vector{<:AbstractBlock{N}}) where N = (prepend!(c1.blocks, list); c1)
Base.prepend!(c1::Sequence, c2::Sequence) where N = (prepend!(c1.blocks, c2.blocks); c1)

max_stars_repo_name: QuantumBFS/YaoExtensions

max_issues_repo_name: QuantumBFS/YaoExtensions

max_forks_repo_name: QuantumBFS/YaoExtensions


text:
export KeepActor, keep, getvalues

"""
    KeepActor{D}() where D

Keep actor provides a storage actor. It saves all incoming successful `next` events in a `values` array, throws an ErrorException on `error!` event
and does nothing on completion event.

# Examples
```jldoctest
using Rocket

source = from(1:5)
actor  = keep(Int)

subscribe!(source, actor)
show(getvalues(actor))

# output
[1, 2, 3, 4, 5]
```

See also: [`Actor`](@ref), [`keep`](@ref)
"""
struct KeepActor{T} <: Actor{T}
    values :: Vector{T}

    KeepActor{T}() where T = new(Vector{T}())
end

getvalues(actor::KeepActor) = actor.values

on_next!(actor::KeepActor{T}, data::T) where T = push!(actor.values, data)
on_error!(actor::KeepActor, err)               = error(err)
on_complete!(actor::KeepActor)                 = begin end

"""
    keep(::Type{T}) where T

# Arguments
- `::Type{T}`: Type of keep data

Creation operator for the `KeepActor` actor.

# Examples

```jldoctest
using Rocket

actor = keep(Int)
actor isa KeepActor{Int}

# output
true
```

See also: [`KeepActor`](@ref), [`AbstractActor`](@ref)
"""
keep(::Type{T}) where T = KeepActor{T}()

# Julia iterable interface

Base.IteratorSize(::Type{ <: KeepActor })   = Base.HasLength()
Base.IteratorEltype(::Type{ <: KeepActor }) = Base.HasEltype()

Base.IndexStyle(::Type{ <: KeepActor }) = Base.IndexLinear()

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

Base.iterate(actor::KeepActor)        = iterate(actor.values)
Base.iterate(actor::KeepActor, state) = iterate(actor.values, state)

Base.size(actor::KeepActor)        = (length(actor.values), )
Base.length(actor::KeepActor)      = length(actor.values)
Base.getindex(actor::KeepActor, I) = Base.getindex(actor.values, I)

Base.firstindex(actor::KeepActor) = firstindex(actor.values)
Base.lastindex(actor::KeepActor)  = lastindex(actor.values)


text:
# Non-deterministic Finite Automaton
# ==================================

struct NFANode
    edges::Vector{Tuple{Edge,NFANode}}
end

function NFANode()
    return NFANode(Tuple{Edge,NFANode}[])
end

function Base.show(io::IO, node::NFANode)
    print(io, summary(node), "(<", length(node.edges), " edges", '@', objectid(node), ">)")
end

struct NFA
    start::NFANode
    final::NFANode

    function NFA(start::NFANode, final::NFANode)
        @assert start !== final
        return new(start, final)
    end
end

# epsilon transition
const eps = ByteSet()

function iseps(e::Edge)
    return isempty(e.labels)
end

function re2nfa(re::RegExp.RE, predefined_actions::Dict{Symbol,Action}=Dict{Symbol,Action}())
    actions = Dict{Tuple{RegExp.RE,Symbol},Action}()
    action_order = 1

    function make_action_list(re, names)
        list = ActionList()
        for name in names
            if haskey(predefined_actions, name)  # pick up a predefined action
                action = predefined_actions[name]
            elseif haskey(actions, (re, name))
                action = actions[(re, name)]
            else
                action = Action(name, action_order)
                actions[(re, name)] = action
                action_order += 1
            end
            push!(list, action)
        end
        return list
    end

    # Thompson's construction.
    function rec!(start, re)
        if haskey(re.actions, :enter)
            start_in = NFANode()
            push!(start.edges, (Edge(eps, make_action_list(re, re.actions[:enter])), start_in))
        else
            start_in = start
        end

        re′ = RegExp.shallow_desugar(re)
        head = re′.head
        args = re′.args

        if head == :set
            @assert length(args) == 1
            final_in = NFANode()
            push!(start_in.edges, (Edge(args[1]), final_in))
        elseif head == :cat
            f = start_in
            for arg in args
                f = rec!(f, arg)
            end
            if f == start_in
                final_in = NFANode()
                push!(start_in.edges, (Edge(eps), final_in))
            else
                final_in = f
            end
        elseif head == :alt
            @assert length(args) > 0
            final_in = NFANode()
            for arg in args
                s = NFANode()
                f = rec!(s, arg)
                push!(start_in.edges, (Edge(eps), s))
                push!(       f.edges, (Edge(eps), final_in))
            end
        elseif head == :rep
            @assert length(args) == 1
            s = NFANode()
            f = rec!(s, args[1])
            final_in = NFANode()
            push!(start_in.edges, (Edge(eps), s))
            push!(start_in.edges, (Edge(eps), final_in))
            push!(       f.edges, (Edge(eps), s))
            push!(       f.edges, (Edge(eps), final_in))
        elseif head == :isec || head == :diff
            @assert length(args) == 2
            final_in = NFANode()
            s1 = NFANode()
            f1 = rec!(s1, args[1])
            push!(start_in.edges, (Edge(eps), s1))
            push!(      f1.edges, (Edge(eps), final_in))
            s2 = NFANode()
            f2 = rec!(s2, args[2])
            push!(start_in.edges, (Edge(eps), s2))
            push!(      f2.edges, (Edge(eps), final_in))
            if head == :isec
                revoke = s -> f1 ∉ s.nfanodes || f2 ∉ s.nfanodes
            else  # re.head == :diff
                revoke = s -> f2 ∈ s.nfanodes
            end
            nfa = dfa2nfa(revoke_finals(revoke, nfa2dfa(NFA(start_in, final_in))))
            push!(start_in.edges, (Edge(eps), nfa.start))
            final_in = nfa.final
        else
            error("unsupported operation: $(head)")
        end

        if haskey(re.actions, :all)
            as = make_action_list(re, re.actions[:all])
            for s in traverse(start), (e, _) in s.edges
                if !iseps(e)
                    union!(e.actions, as)
                end
            end
        end

        if haskey(re.actions, :final)
            as = make_action_list(re, re.actions[:final])
            for s in traverse(start), (e, t) in s.edges
                if !iseps(e) && final_in ∈ epsilon_closure(t)
                    union!(e.actions, as)
                end
            end
        end

        if haskey(re.actions, :exit)
            final = NFANode()
            push!(final_in.edges, (Edge(eps, make_action_list(re, re.actions[:exit])), final))
        else
            final = final_in
        end

        if re.when != nothing
            name = re.when
            for s in traverse(start), (e, _) in s.edges
                if !iseps(e)
                    push!(e.precond, (name => TRUE))
                end
            end
        end

        return final
    end

    nfa_start = NFANode()
    nfa_final = rec!(nfa_start, re)
    return NFA(nfa_start, nfa_final)
end

function remove_dead_nodes(nfa::NFA)
    backrefs = Dict(nfa.start => Set{NFANode}())
    for s in traverse(nfa.start), (_, t) in s.edges
        push!(get!(() -> Set{NFANode}(), backrefs, t), s)
    end

    alive = Set{NFANode}()
    unvisited = [nfa.final]
    while !isempty(unvisited)
        s = pop!(unvisited)
        push!(alive, s)
        for t in backrefs[s]
            if t ∉ alive
                push!(unvisited, t)
            end
        end
    end
    @assert nfa.start ∈ alive
    @assert nfa.final ∈ alive

    newnodes = Dict{NFANode,NFANode}()
    new(s) = get!(() -> NFANode(), newnodes, s)
    isvisited(s) = haskey(newnodes, s)
    unvisited = [nfa.start]
    while !isempty(unvisited)
        s = pop!(unvisited)
        s′ = new(s)
        for (e, t) in s.edges
            if t ∈ alive
                if !isvisited(t)
                    push!(unvisited, t)
                end
                push!(s′.edges, (e, new(t)))
            end
        end
    end

    return NFA(new(nfa.start), new(nfa.final))
end