text:
# Predicates on rows which can hopefully take advantage of acceleration indexes

abstract type Predicate{names} <: Function; end # Intention to add `!`, `&`, `|` on `Predicate`s.
# TODO: logical not, and, or, etc for predicates?

# Predicates are useful for filtering, etc:

function map(pred::Predicate{names}, t::Table) where {names}
    # First get the indices using the acceleration indices
    t_projected = Project(names)(t)
    index = promote_index(accelerators(t_projected)...)
    return _map(pred, t_projected, index)
end

function _map(pred::Predicate{names}, t::Table{names}, ::AbstractIndex) where {names}
    return map(row -> pred(row), t) # Use Julia's default algorithm
end

function findall(pred::Predicate{names}, t::Table) where {names}
    # First get the indices using the acceleration indices
    t_projected = Project(names)(t)
    index = promote_index(accelerators(t_projected)...)
    return _findall(pred, t_projected, index)
end

function _findall(pred::Predicate{names}, t::Table{names}, ::AbstractIndex) where {names}
    findall(row -> pred(row), t) # Use Julia's default algorithm
end

function filter(pred::Predicate{names}, t::Table) where {names}
    # First get the indices using the acceleration indices
    t_projected = Project(names)(t)
    index = promote_index(accelerators(t_projected)...)
    inds = _filter_indices(pred, t_projected, index) # Uses `map` or `findall` depending on available acceleration indices
    @inbounds return t[inds]
end

function _filter_indices(pred::Predicate{names}, t::Table{names}, index::AbstractIndex) where {names}
    return _map(pred, t, index) # Like Julia, default to map
end

# TODO similarly for findfirst, findlast, findnext, findprev, findmin, findmax

"""
    IsEqual(namedtuple)
    IsEqual(name = value, ...)

Creates an `IsEqual` function, which returns true on any named tuple whose fields of the
given name(s) are `isequal` to those specified.

See also `Equals` for comparing two columns of the same table.
"""
struct IsEqual{names, D <: Tuple} <: Predicate{names}
    data::D
end
IsEqual(;kwargs...) = IsEqual(kwargs.data)
IsEqual(nt::NamedTuple{names}) where {names} = IsEqual{names}(Tuple(nt))
IsEqual{names}(t::T) where {names, T <: Tuple} = IsEqual{names, T}(t)

function (pred::IsEqual{names})(x::NamedTuple{names2}) where {names, names2}
    if names === names2
        return isequal(pred.data, Tuple(x))
    else
        return pred(Project(names)(x))
    end
end

# map for IsEqual

function _map(pred::IsEqual{names}, t::Table{names}, ::UniqueIndex) where {names}
    out = fill(false, length(t))
    i = findfirst(pred, t)
    if i === nothing
        return out
    else
        @inbounds out[i] = true
        return out
    end
end

function _map(pred::IsEqual{names}, t::Table{names}, index::SortIndex{names}) where {names}
    out = fill(false, length(t))

    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = searchsorted(view(t, index.order), searchrow)
    @inbounds out[index.order[range]] = true

    return out
end

function _map(pred::IsEqual{names}, t::Table{names}, index::SortIndex{names2}) where {names, names2}
    out = fill(false, length(t))

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds range = searchsorted(view(t_projected, index.order), searchrow)
    @inbounds out[index.order[range]] = true

    return out
end

function _map(pred::IsEqual{names}, t::Table{names}, index::UniqueSortIndex{names}) where {names}
    n = length(t)
    out = fill(false, n)

    searchrow = NamedTuple{names}(pred.data)
    @inbounds first_greater_or_equal = searchsortedfirst(view(t, index.order), searchrow)
    if first_greater_or_equal <= n
        @inbounds out[index.order[first_greater_or_equal]] = pred(t[index.order[first_greater_or_equal]])
    end

    return out
end

function _map(pred::IsEqual{names}, t::Table{names}, index::UniqueSortIndex{names2}) where {names, names2}
    n = length(t)
    out = fill(false, n)

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds first_greater_or_equal = searchsortedfirst(view(t_projected, index.order), searchrow)
    if first_greater_or_equal <= n
        @inbounds out[index.order[first_greater_or_equal]] = pred(t[index.order[first_greater_or_equal]])
    end

    return out
end

function _map(pred::IsEqual{names}, t::Table{names}, index::HashIndex{names}) where {names}
    out = fill(false, length(t))

    key = pred.data
    if haskey(index.dict, key) # TODO make faster
        inds = index.dict[key]
        @inbounds out[inds] = true
    end

    return out
end

function _map(pred::IsEqual{names}, t::Table{names}, index::HashIndex{names2}) where {names, names2}
    out = fill(false, length(t))

    key = Tuple(Project(names2)(NamedTuple{names}(pred.data)))
    if haskey(index.dict, key)
        inds = index.dict[key] # TODO make faster
        @inbounds for i in inds
            if pred(t[i])
                out[i] = true
            end
        end
    end

    return out
end

function _map(pred::IsEqual{names}, t::Table{names}, index::UniqueHashIndex{names}) where {names}
    out = fill(false, length(t))

    key = pred.data
    if haskey(index.dict, key) # TODO make faster
        i = index.dict[key]
        @inbounds out[i] = true
    end

    return out
end

function _map(pred::IsEqual{names}, t::Table{names}, index::UniqueHashIndex{names2}) where {names, names2}
    out = fill(false, length(t))
    key = Tuple(Project(names2)(NamedTuple{names}(pred.data)))

    if haskey(index.dict, key) # TODO make faster
        i = index.dict[key]
        @inbounds out[i] = pred(t[i])
    end

    return out
end

# findall for IsEqual

function _findall(pred::IsEqual{names}, t::Table{names}, ::UniqueIndex) where {names}
    i = findfirst(pred, t)
    if i === nothing
        return 1:0
    else
        return i:i
    end
end

function _findall(pred::IsEqual{names}, t::Table{names}, index::SortIndex{names}) where {names}
    searchrow = NamedTuple{names}(pred.data)
    return @inbounds index.order[searchsorted(view(t, index.order), searchrow)]
end

function _findall(pred::IsEqual{names}, t::Table{names}, index::SortIndex{names2}) where {names, names2}
    searchrow = Project(names2)(NamedTuple{names}(pred.data))
    t_projected = Project(names2)(t)
    @inbounds range = searchsorted(view(t_projected, index.order), searchrow)
    out = Int[]
    @inbounds for i ∈ range
        j = index.order[i]
        if pred(t[j])
            push!(out, j)
        end
    end

    return out
end

function _findall(pred::IsEqual{names}, t::Table{names}, index::UniqueSortIndex{names}) where {names}
    n = length(t)
    searchrow = NamedTuple{names}(pred.data)
    @inbounds first_greater_or_equal = searchsortedfirst(view(t, index.order), searchrow)
    @inbounds if first_greater_or_equal <= n
        i = index.order[first_greater_or_equal]
        if pred(t[i])
            return i:i
        end
    end
    return 1:0
end

function _findall(pred::IsEqual{names}, t::Table{names}, index::UniqueSortIndex{names2}) where {names, names2}
    n = length(t)
    searchrow = Project(names2)(NamedTuple{names}(pred.data))
    t_projected = Project(names2)(t)
    @inbounds first_greater_or_equal = searchsortedfirst(view(t_projected, index.order), searchrow)
    @inbounds if first_greater_or_equal <= n
        i = index.order[first_greater_or_equal]
        if pred(t[i])
            return i:i
        end
    end
    return 1:0
end

function _findall(pred::IsEqual{names}, t::Table{names}, index::HashIndex{names}) where {names}
    searchrow = pred.data
    return get(() -> Int[], index.dict, searchrow)
end

function _findall(pred::IsEqual{names}, t::Table{names}, index::HashIndex{names2}) where {names, names2}
    searchrow = Tuple(Project(names2)(NamedTuple{names}(pred.data)))
    inds = get(() -> Int[], index.dict, searchrow)
    if length(inds) == 0
        return inds
    end
    return filter(i -> @inbounds(pred(t[i])), inds)
end

function _findall(pred::IsEqual{names}, t::Table{names}, index::UniqueHashIndex{names}) where {names}
    searchrow = pred.data
    i = get(() -> 0, index.dict, searchrow)
    if i > 0
        return i:i
    else
        return 1:0
    end
end

function _findall(pred::IsEqual{names}, t::Table{names}, index::UniqueHashIndex{names2}) where {names, names2}
    searchrow = Tuple(Project(names2)(NamedTuple{names}(pred.data)))
    i = get(() -> 0, index.dict, searchrow)
    if i > 0 && @inbounds(pred(t[i]))
        return i:i
    else
        return 1:0
    end
end

function _filter_indices(pred::IsEqual{names}, t::Table{names}, index::AbstractUniqueIndex) where {names}
    return _findall(pred, t, index)
end

function _filter_indices(pred::IsEqual{names}, t::Table{names}, index::SortIndex{names}) where {names}
    return _findall(pred, t, index)
end

function _filter_indices(pred::IsEqual{names}, t::Table{names}, index::HashIndex{names}) where {names}
    return _findall(pred, t, index)
end

"""
    IsLess(namedtuple)
    IsLess(name = value, ...)

Creates an `IsLess` function, which returns true on any named tuple whose fields of the
given name(s) are `isless` to those specified.
"""
struct IsLess{names, D <: Tuple} <: Predicate{names}
    data::D
end
IsLess(;kwargs...) = IsLess(kwargs.data)
IsLess(nt::NamedTuple{names}) where {names} = IsLess{names}(Tuple(nt))
IsLess{names}(t::T) where {names, T <: Tuple} = IsLess{names, T}(t)

function (pred::IsLess{names})(x::NamedTuple{names2}) where {names, names2}
    if names === names2
        return isless(Tuple(x), pred.data)
    else
        return pred(Project(names)(x))
    end
end

# map for IsLess

function _map(pred::IsLess{names}, t::Table{names}, index::SortIndex{names}) where {names}
    out = fill(false, length(t))

    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = 1:searchsortedlastless(view(t, index.order), searchrow)
    @inbounds out[view(index.order, range)] = true

    return out
end

function _map(pred::IsLess{names}, t::Table{names}, index::SortIndex{names2}) where {names, names2}
    # Be careful about what index you accept - TODO perhaps `promote_index` needs to know the predicate rather than do this here
    names3 = _headsubset(Val(names), Val(names2))
    if names3 === ()
        return _map(pred, t, NoIndex())
    elseif names3 !== names2
        return _map(pred, t, SortIndex{names3}(index.order))
    end

    out = fill(false, length(t))

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds range = searchsorted(view(t_projected, index.order), searchrow)
    @inbounds out[index.order[1:range.start-1]] = true
    @inbounds for i ∈ range
        j = index.order[i]
        out[j] = pred(t[j]) # TODO: Make faster by comparing only columns setdiff(names, names2)
    end

    return out
end

function _map(pred::IsLess{names}, t::Table{names}, index::UniqueSortIndex{names}) where {names}
    out = fill(false, length(t))

    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = 1:searchsortedlastless(view(t, index.order), searchrow)
    @inbounds out[view(index.order, range)] = true

    return out
end

function _map(pred::IsLess{names}, t::Table{names}, index::UniqueSortIndex{names2}) where {names, names2}
    # Be careful about what index you accept - TODO perhaps `promote_index` needs to know the predicate rather than do this here
    names3 = _headsubset(Val(names), Val(names2))
    if names3 === ()
        return _map(pred, t, NoIndex())
    elseif names3 !== names2
        return _map(pred, t, SortIndex{names3}(index.order)) # potentially looses uniqueness
    end

    n = length(t)
    out = fill(false, length(t))

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds i = searchsortedfirst(view(t_projected, index.order), searchrow)
    @inbounds if i > n || !pred(t[index.order[i]]) # TODO: Make faster by comparing only columns setdiff(names, names2)
        i = i - 1
    end
    @inbounds out[index.order[1:i]] = true

    return out
end

# findall for IsLess

function _findall(pred::IsLess{names}, t::Table{names}, index::SortIndex{names}) where {names}
    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = 1:searchsortedlastless(view(t, index.order), searchrow)

    @inbounds return sort(view(index.order, range))
end

function _findall(pred::IsLess{names}, t::Table{names}, index::SortIndex{names2}) where {names, names2}
    # Be careful about what index you accept - TODO perhaps `promote_index` needs to know the predicate rather than do this here
    names3 = _headsubset(Val(names), Val(names2))
    if names3 === ()
        return _findall(pred, t, NoIndex())
    elseif names3 !== names2
        return _findall(pred, t, SortIndex{names3}(index.order))
    end

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds range = searchsorted(view(t_projected, index.order), searchrow)
    @inbounds out = index.order[1:range.start-1]
    @inbounds for i ∈ range
        j = index.order[i]
        if pred(t[j]) # TODO: Make faster by comparing only columns setdiff(names, names2)
            push!(out, j)
        end
    end

    sort!(out)
    return out
end

function _findall(pred::IsLess{names}, t::Table{names}, index::UniqueSortIndex{names}) where {names}
    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = 1:searchsortedlastless(view(t, index.order), searchrow)
    @inbounds return sort(view(index.order, range))
end

function _findall(pred::IsLess{names}, t::Table{names}, index::UniqueSortIndex{names2}) where {names, names2}
    # Be careful about what index you accept - TODO perhaps `promote_index` needs to know the predicate rather than do this here
    names3 = _headsubset(Val(names), Val(names2))
    if names3 === ()
        return _findall(pred, t, NoIndex())
    elseif names3 !== names2
        return _findall(pred, t, SortIndex{names3}(index.order)) # potentially looses uniqueness
    end

    n = length(t)

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds i = searchsortedfirst(view(t_projected, index.order), searchrow)
    @inbounds if i > n || !pred(t[index.order[i]]) # TODO: Make faster by comparing only columns setdiff(names, names2)
        i = i - 1
    end
    @inbounds return sort(view(index.order, 1:i))
end

# filter for IsLess - use default `map` method

"""
    IsLessEqual(namedtuple)
    IsLessEqual(name = value, ...)

Creates an `IsLessEqual` function, which returns true on any named tuple whose fields of the
given name(s) are `isless` or `isequal` to those specified.
"""
struct IsLessEqual{names, D <: Tuple} <: Predicate{names}
    data::D
end
IsLessEqual(;kwargs...) = IsLessEqual(kwargs.data)
IsLessEqual(nt::NamedTuple{names}) where {names} = IsLessEqual{names}(Tuple(nt))
IsLessEqual{names}(t::T) where {names, T <: Tuple} = IsLessEqual{names, T}(t)

function (pred::IsLessEqual{names})(x::NamedTuple{names2}) where {names, names2}
    if names === names2
        return !isless(pred.data, Tuple(x))
    else
        return pred(Project(names)(x))
    end
end

# map for IsLessEqual

function _map(pred::IsLessEqual{names}, t::Table{names}, index::SortIndex{names}) where {names}
    out = fill(false, length(t))

    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = 1:searchsortedlast(view(t, index.order), searchrow)
    @inbounds out[view(index.order, range)] = true

    return out
end

function _map(pred::IsLessEqual{names}, t::Table{names}, index::SortIndex{names2}) where {names, names2}
    # Be careful about what index you accept - TODO perhaps `promote_index` needs to know the predicate rather than do this here
    names3 = _headsubset(Val(names), Val(names2))
    if names3 === ()
        return _map(pred, t, NoIndex())
    elseif names3 !== names2
        return _map(pred, t, SortIndex{names3}(index.order))
    end

    out = fill(false, length(t))

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds range = searchsorted(view(t_projected, index.order), searchrow)
    @inbounds out[index.order[1:range.start-1]] = true
    @inbounds for i ∈ range
        j = index.order[i]
        out[j] = pred(t[j]) # TODO: Make faster by comparing only columns setdiff(names, names2)
    end

    return out
end

function _map(pred::IsLessEqual{names}, t::Table{names}, index::UniqueSortIndex{names}) where {names}
    out = fill(false, length(t))

    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = 1:searchsortedlast(view(t, index.order), searchrow)
    @inbounds out[view(index.order, range)] = true

    return out
end

function _map(pred::IsLessEqual{names}, t::Table{names}, index::UniqueSortIndex{names2}) where {names, names2}
    # Be careful about what index you accept - TODO perhaps `promote_index` needs to know the predicate rather than do this here
    names3 = _headsubset(Val(names), Val(names2))
    if names3 === ()
        return _map(pred, t, NoIndex())
    elseif names3 !== names2
        return _map(pred, t, SortIndex{names3}(index.order)) # potentially looses uniqueness
    end

    n = length(t)
    out = fill(false, length(t))

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds i = searchsortedfirst(view(t_projected, index.order), searchrow)
    @inbounds if i > n || !pred(t[index.order[i]]) # TODO: Make faster by comparing only columns setdiff(names, names2)
        i = i - 1
    end
    @inbounds out[index.order[1:i]] = true

    return out
end

# findall for IsLessEqual

function _findall(pred::IsLessEqual{names}, t::Table{names}, index::SortIndex{names}) where {names}
    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = 1:searchsortedlast(view(t, index.order), searchrow)

    @inbounds return sort(view(index.order, range))
end

function _findall(pred::IsLessEqual{names}, t::Table{names}, index::SortIndex{names2}) where {names, names2}
    # Be careful about what index you accept - TODO perhaps `promote_index` needs to know the predicate rather than do this here
    names3 = _headsubset(Val(names), Val(names2))
    if names3 === ()
        return _findall(pred, t, NoIndex())
    elseif names3 !== names2
        return _findall(pred, t, SortIndex{names3}(index.order))
    end

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds range = searchsorted(view(t_projected, index.order), searchrow)
    @inbounds out = index.order[1:range.start-1]
    @inbounds for i ∈ range
        j = index.order[i]
        if pred(t[j]) # TODO: Make faster by comparing only columns setdiff(names, names2)
            push!(out, j)
        end
    end

    sort!(out)
    return out
end

function _findall(pred::IsLessEqual{names}, t::Table{names}, index::UniqueSortIndex{names}) where {names}
    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = 1:searchsortedlast(view(t, index.order), searchrow)
    @inbounds return sort(view(index.order, range))
end

function _findall(pred::IsLessEqual{names}, t::Table{names}, index::UniqueSortIndex{names2}) where {names, names2}
    # Be careful about what index you accept - TODO perhaps `promote_index` needs to know the predicate rather than do this here
    names3 = _headsubset(Val(names), Val(names2))
    if names3 === ()
        return _findall(pred, t, NoIndex())
    elseif names3 !== names2
        return _findall(pred, t, SortIndex{names3}(index.order)) # potentially looses uniqueness
    end

    n = length(t)

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds i = searchsortedfirst(view(t_projected, index.order), searchrow)
    @inbounds if i > n || !pred(t[index.order[i]]) # TODO: Make faster by comparing only columns setdiff(names, names2)
        i = i - 1
    end
    @inbounds return sort(view(index.order, 1:i))
end

"""
    IsGreater(namedtuple)
    IsGreater(name = value, ...)

Creates an `IsGreater` function, which returns true on any named tuple whose fields of the
given name(s) are not `isless` or `isequal` to those specified.
"""
struct IsGreater{names, D <: Tuple} <: Predicate{names}
    data::D
end
IsGreater(;kwargs...) = IsGreater(kwargs.data)
IsGreater(nt::NamedTuple{names}) where {names} = IsGreater{names}(Tuple(nt))
IsGreater{names}(t::T) where {names, T <: Tuple} = IsGreater{names, T}(t)

function (pred::IsGreater{names})(x::NamedTuple{names2}) where {names, names2}
    if names === names2
        return isless(pred.data, Tuple(x))
    else
        return pred(Project(names)(x))
    end
end

# map for IsGreater

function _map(pred::IsGreater{names}, t::Table{names}, index::SortIndex{names}) where {names}
    out = fill(false, length(t))

    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = searchsortedfirstgreater(view(t, index.order), searchrow):length(t)
    @inbounds out[view(index.order, range)] = true

    return out
end

function _map(pred::IsGreater{names}, t::Table{names}, index::SortIndex{names2}) where {names, names2}
    # Be careful about what index you accept - TODO perhaps `promote_index` needs to know the predicate rather than do this here
    names3 = _headsubset(Val(names), Val(names2))
    if names3 === ()
        return _map(pred, t, NoIndex())
    elseif names3 !== names2
        return _map(pred, t, SortIndex{names3}(index.order))
    end

    out = fill(false, length(t))

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds range = searchsorted(view(t_projected, index.order), searchrow)
    @inbounds for i ∈ range
        j = index.order[i]
        out[j] = pred(t[j]) # TODO: Make faster by comparing only columns setdiff(names, names2)
    end
    @inbounds out[index.order[range.stop+1:length(t)]] = true

    return out
end

function _map(pred::IsGreater{names}, t::Table{names}, index::UniqueSortIndex{names}) where {names}
    out = fill(false, length(t))

    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = searchsortedfirstgreater(view(t, index.order), searchrow):length(t)
    @inbounds out[view(index.order, range)] = true

    return out
end

function _map(pred::IsGreater{names}, t::Table{names}, index::UniqueSortIndex{names2}) where {names, names2}
    # Be careful about what index you accept - TODO perhaps `promote_index` needs to know the predicate rather than do this here
    names3 = _headsubset(Val(names), Val(names2))
    if names3 === ()
        return _map(pred, t, NoIndex())
    elseif names3 !== names2
        return _map(pred, t, SortIndex{names3}(index.order)) # potentially looses uniqueness
    end

    n = length(t)
    out = fill(false, length(t))

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds i = searchsortedlast(view(t_projected, index.order), searchrow)
    @inbounds if i == 0 || !pred(t[index.order[i]]) # TODO: Make faster by comparing only columns setdiff(names, names2)
        i = i + 1
    end
    @inbounds out[index.order[i:n]] = true

    return out
end

# findall for IsGreater

function _findall(pred::IsGreater{names}, t::Table{names}, index::SortIndex{names}) where {names}
    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = searchsortedfirstgreater(view(t, index.order), searchrow):length(t)

    @inbounds return sort(view(index.order, range))
end

function _findall(pred::IsGreater{names}, t::Table{names}, index::SortIndex{names2}) where {names, names2}
    # Be careful about what index you accept - TODO perhaps `promote_index` needs to know the predicate rather than do this here
    names3 = _headsubset(Val(names), Val(names2))
    if names3 === ()
        return _findall(pred, t, NoIndex())
    elseif names3 !== names2
        return _findall(pred, t, SortIndex{names3}(index.order))
    end

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds range = searchsorted(view(t_projected, index.order), searchrow)
    out = Int[]
    @inbounds for i ∈ range
        j = index.order[i]
        if pred(t[j]) # TODO: Make faster by comparing only columns setdiff(names, names2)
            push!(out, j)
        end
    end
    @inbounds append!(out, index.order[range.stop+1:length(t)])

    sort!(out)
    return out
end

function _findall(pred::IsGreater{names}, t::Table{names}, index::UniqueSortIndex{names}) where {names}
    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = searchsortedfirstgreater(view(t, index.order), searchrow):length(t)
    @inbounds return sort(view(index.order, range))
end

function _findall(pred::IsGreater{names}, t::Table{names}, index::UniqueSortIndex{names2}) where {names, names2}
    # Be careful about what index you accept - TODO perhaps `promote_index` needs to know the predicate rather than do this here
    names3 = _headsubset(Val(names), Val(names2))
    if names3 === ()
        return _findall(pred, t, NoIndex())
    elseif names3 !== names2
        return _findall(pred, t, SortIndex{names3}(index.order)) # potentially looses uniqueness
    end

    n = length(t)

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds i = searchsortedlast(view(t_projected, index.order), searchrow)
    @inbounds if i == 0 || !pred(t[index.order[i]]) # TODO: Make faster by comparing only columns setdiff(names, names2)
        i = i + 1
    end
    @inbounds return sort(view(index.order, i:n))
end

"""
    IsGreaterEqual(namedtuple)
    IsGreaterEqual(name = value, ...)

Creates an `IsGreaterEqual` function, which returns true on any named tuple whose fields of the
given name(s) are not `isless` to those specified.
"""
struct IsGreaterEqual{names, D <: Tuple} <: Predicate{names}
    data::D
end
IsGreaterEqual(;kwargs...) = IsGreaterEqual(kwargs.data)
IsGreaterEqual(nt::NamedTuple{names}) where {names} = IsGreaterEqual{names}(Tuple(nt))
IsGreaterEqual{names}(t::T) where {names, T <: Tuple} = IsGreaterEqual{names, T}(t)

function (pred::IsGreaterEqual{names})(x::NamedTuple{names2}) where {names, names2}
    if names === names2
        return !isless(Tuple(x), pred.data)
    else
        return pred(Project(names)(x))
    end
end

# map for IsGreaterEqual

function _map(pred::IsGreaterEqual{names}, t::Table{names}, index::SortIndex{names}) where {names}
    out = fill(false, length(t))

    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = searchsortedfirst(view(t, index.order), searchrow):length(t)
    @inbounds out[view(index.order, range)] = true

    return out
end

function _map(pred::IsGreaterEqual{names}, t::Table{names}, index::SortIndex{names2}) where {names, names2}
    # Be careful about what index you accept - TODO perhaps `promote_index` needs to know the predicate rather than do this here
    names3 = _headsubset(Val(names), Val(names2))
    if names3 === ()
        return _map(pred, t, NoIndex())
    elseif names3 !== names2
        return _map(pred, t, SortIndex{names3}(index.order))
    end

    out = fill(false, length(t))

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds range = searchsorted(view(t_projected, index.order), searchrow)
    @inbounds for i ∈ range
        j = index.order[i]
        out[j] = pred(t[j]) # TODO: Make faster by comparing only columns setdiff(names, names2)
    end
    @inbounds out[index.order[range.stop+1:length(t)]] = true

    return out
end

function _map(pred::IsGreaterEqual{names}, t::Table{names}, index::UniqueSortIndex{names}) where {names}
    out = fill(false, length(t))

    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = searchsortedfirst(view(t, index.order), searchrow):length(t)
    @inbounds out[view(index.order, range)] = true

    return out
end

function _map(pred::IsGreaterEqual{names}, t::Table{names}, index::UniqueSortIndex{names2}) where {names, names2}
    # Be careful about what index you accept - TODO perhaps `promote_index` needs to know the predicate rather than do this here
    names3 = _headsubset(Val(names), Val(names2))
    if names3 === ()
        return _map(pred, t, NoIndex())
    elseif names3 !== names2
        return _map(pred, t, SortIndex{names3}(index.order)) # potentially looses uniqueness
    end

    n = length(t)
    out = fill(false, length(t))

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds i = searchsortedlast(view(t_projected, index.order), searchrow)
    @inbounds if i == 0 || !pred(t[index.order[i]]) # TODO: Make faster by comparing only columns setdiff(names, names2)
        i = i + 1
    end
    @inbounds out[index.order[i:n]] = true

    return out
end

# findall for IsGreaterEqual

function _findall(pred::IsGreaterEqual{names}, t::Table{names}, index::SortIndex{names}) where {names}
    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = searchsortedfirst(view(t, index.order), searchrow):length(t)

    @inbounds return sort(view(index.order, range))
end

function _findall(pred::IsGreaterEqual{names}, t::Table{names}, index::SortIndex{names2}) where {names, names2}
    # Be careful about what index you accept - TODO perhaps `promote_index` needs to know the predicate rather than do this here
    names3 = _headsubset(Val(names), Val(names2))
    if names3 === ()
        return _findall(pred, t, NoIndex())
    elseif names3 !== names2
        return _findall(pred, t, SortIndex{names3}(index.order))
    end

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds range = searchsorted(view(t_projected, index.order), searchrow)
    out = Int[]
    @inbounds for i ∈ range
        j = index.order[i]
        if pred(t[j]) # TODO: Make faster by comparing only columns setdiff(names, names2)
            push!(out, j)
        end
    end
    @inbounds append!(out, index.order[range.stop+1:length(t)])

    sort!(out)
    return out
end

function _findall(pred::IsGreaterEqual{names}, t::Table{names}, index::UniqueSortIndex{names}) where {names}
    searchrow = NamedTuple{names}(pred.data)
    @inbounds range = searchsortedfirst(view(t, index.order), searchrow):length(t)
    @inbounds return sort(view(index.order, range))
end

function _findall(pred::IsGreaterEqual{names}, t::Table{names}, index::UniqueSortIndex{names2}) where {names, names2}
    # Be careful about what index you accept - TODO perhaps `promote_index` needs to know the predicate rather than do this here
    names3 = _headsubset(Val(names), Val(names2))
    if names3 === ()
        return _findall(pred, t, NoIndex())
    elseif names3 !== names2
        return _findall(pred, t, SortIndex{names3}(index.order)) # potentially looses uniqueness
    end

    n = length(t)

    searchrow = Project{names2}()(NamedTuple{names}(pred.data))
    t_projected = Project{names2}()(t)
    @inbounds i = searchsortedlast(view(t_projected, index.order), searchrow)
    @inbounds if i == 0 || !pred(t[index.order[i]]) # TODO: Make faster by comparing only columns setdiff(names, names2)
        i = i + 1
    end
    @inbounds return sort(view(index.order, i:n))
end

"""
    In(namedtuple)
    In(name = collection, ...)

Creates an `In` function, which returns true on any named tuple whose field(s) of the
given name(s) is `in` the collection(s) specified.
"""
struct In{names, T <: Tuple} <: Predicate{names}
    data::T
end
In(;kwargs...) = In(kwargs.data)
In(nt::NamedTuple{names}) where {names} = In{names}(Tuple(nt))
In{names}(t::T) where {names, T <: Tuple} = In{names, T}(t)

function (pred::In{names})(x::NamedTuple{names2}) where {names, names2}
    if names === names2
        return _all(in, Tuple(x), pred.data)
    else
        return pred(Project(names)(x))
    end
end

# map for `In` of `Interval`

function _map(pred::In{names, Tuple{<:Interval}}, t::Table{names}, index::SortIndex{names}) where {names}
    # TODO Currently expects a single column here - generalize
    # TODO The generalized solution could use a generalized search tree like k-d tree
    out = fill(false, length(t))

    startrow = NamedTuple{names}(pred.data[1].start)
    @inbounds start = searchsortedfirst(view(t, index.order), startrow)

    stoprow = NamedTuple{names}(pred.data[1].start)
    @inbounds stop = searchsortedlast(view(t, index.order), stoprow)

    @inbounds out[view(index.order, start:stop)] = true

    return out
end

function _map(pred::In{names, Tuple{<:Interval}}, t::Table{names}, index::UniqueSortIndex{names}) where {names}
    # TODO Currently expects a single column here - generalize
    # TODO The generalized solution could use a generalized search tree like k-d tree
    out = fill(false, length(t))

    startrow = NamedTuple{names}(pred.data[1].start)
    @inbounds start = searchsortedfirst(view(t, index.order), startrow)

    stoprow = NamedTuple{names}(pred.data[1].start)
    @inbounds stop = searchsortedlast(view(t, index.order), stoprow)

    @inbounds out[view(index.order, start:stop)] = true

    return out
end

# findall for `In` of `Interval`

function _findall(pred::In{names, Tuple{<:Interval}}, t::Table{names}, index::SortIndex{names}) where {names}
    # TODO Currently expects a single column here - generalize
    # TODO The generalized solution could use a generalized search tree like k-d tree
    startrow = NamedTuple{names}(pred.data[1].start)
    @inbounds start = searchsortedfirst(view(t, index.order), startrow)

    stoprow = NamedTuple{names}(pred.data[1].start)
    @inbounds stop = searchsortedlast(view(t, index.order), stoprow)

    @inbounds return index.order[start:stop]
end

function _findall(pred::In{names, Tuple{<:Interval}}, t::Table{names}, index::UniqueSortIndex{names}) where {names}
    # TODO Currently expects a single column here - generalize
    # TODO The generalized solution could use a generalized search tree like k-d tree
    startrow = NamedTuple{names}(pred.data[1].start)
    @inbounds start = searchsortedfirst(view(t, index.order), startrow)

    stoprow = NamedTuple{names}(pred.data[1].start)
    @inbounds stop = searchsortedlast(view(t, index.order), stoprow)

    @inbounds return index.order[start:stop]
end

# filter for `In` of `Interval`

function _filter_indices(pred::IsEqual{names, Tuple{<:Interval}}, t::Table{names}, index::SortIndex{names}) where {names}
    return _findall(pred, t, index)
end

function _filter_indices(pred::IsEqual{names, Tuple{<:Interval}}, t::Table{names}, index::UniqueSortIndex{names}) where {names}
    return _findall(pred, t, index)
end

"""
    Equals(name1, name2)

Creates an `Equals` function, which returns true on any named tuple whose fields of the
two given names are `isequal` to each other.

See also `IsEqual` for comparing one or more columns of a table to a fixed value.
"""
struct Equals{names} <: Predicate{names}
end
@inline Equals(n1::Symbol, n2::Symbol) = Equals{(n1, n2)}()
@inline Equals(names::Tuple{Symbol,Symbol}) = Equals{names}()

function (pred::Equals{names})(x::NamedTuple{names2}) where {names, names2}
    if names === names2
        return isequal(getproperty(x, names[1]), getproperty(x, names[2]))
    else
        return pred(Project(names)(x))
    end
end

function Predicate(pred::Equals{names}, x::NamedTuple{names2}) where {names, names2}
    if names[1] === names2[1]
        return IsEqual(NamedTuple{(names[2],)}(Tuple(x)))
    elseif names[2] === names2[1]
        return IsEqual(NamedTuple{(names[1],)}(Tuple(x)))
    else
        error("Can't create conditional predicate from $pred and $x")
    end
end

# function _map(pred::Equals{names}, t::Table{names}, index::SortIndex{names}) where {names}
#     # TODO sort-merge join algorithm could benefit from `searchsortednext`

#     # TODO finish this!
#     return _map(pred, t, SortIndex{(names[2],)}())

#     # out = fill(false, length(t))
#     # if length(t) == 0
#     #     return out
#     # end

#     # cols = columns(t)
#     # v1 = getproperty(cols, names[1])
#     # v2 = getproperty(cols, names[2])

#     # # The sort-merge join algorithm
#     # i1 = 1
#     # i2 = 0
#     # n = length(t)
#     # @inbounds x = v1[1]
#     # while i1 < n && i2 < n
#     #     if i2 == i1
#     #         @inbounds out[i1] = true
#     #         i1 += 1
#     #     elseif i2 < i1
#     #         i2 = searchsortedfirst(v2, x)
#     #         @inbounds x = v2[i2]
#     #     else

#     #     end

#     #     @inbounds x = v2[1]

#     # end

#     # return out
# end

# function _map(pred::Equals{names}, t::Table{names}, index::SortIndex{names2}) where {names, names2}
#     if length(names2) !== 0 && length(names2) !== 1
#         if names2[1] === names[2] && names2[2] === names[1]
#             p = Project((names[2], names[1]))
#             _map(Equals(names[2], names[1]), p(t), index)
#         end
#     end

#     return _map(pred, t, NoIndex())
# end

"""
    LessThan(name1, name2)

Creates an `LessThan` function, which returns true on any named tuple whose field `name1`
is `isless` to the field `name2`.

See also `IsLess` for comparing one or more columns of a table to a fixed value.
"""
struct LessThan{names} <: Predicate{names}
end
@inline LessThan(n1::Symbol, n2::Symbol) = LessThan{(n1, n2)}()
@inline LessThan(names::Tuple{Symbol,Symbol}) = LessThan{names}()

function (pred::LessThan{names})(x::NamedTuple{names2}) where {names, names2}
    if names === names2
        return isless(getproperty(x, names[1]), getproperty(x, names[2]))
    else
        return pred(Project(names)(x))
    end
end

function Predicate(pred::LessThan{names}, x::NamedTuple{names2}) where {names, names2}
    if names[1] === names2[1]
        return IsGreaterEqual(NamedTuple{(names[2],)}(Tuple(x)))
    elseif names[2] === names2[1]
        return IsLess(NamedTuple{(names[1],)}(Tuple(x)))
    else
        error("Can't create conditional predicate from $pred and $x")
    end
end

"""
    LessEqualThan(name1, name2)

Creates an `LessEqualThan` function, which returns true on any named tuple whose field `name1`
is `isless` or `isequal` to the field `name2`.

See also `IsLessEqual` for comparing one or more columns of a table to a fixed value.
"""
struct LessEqualThan{names} <: Predicate{names}
end
@inline LessEqualThan(n1::Symbol, n2::Symbol) = LessEqualThan{(n1, n2)}()
@inline LessEqualThan(names::Tuple{Symbol,Symbol}) = LessEqualThan{names}()

function (pred::LessEqualThan{names})(x::NamedTuple{names2}) where {names, names2}
    if names === names2
        return !isless(getproperty(x, names[2]), getproperty(x, names[1]))
    else
        return pred(Project(names)(x))
    end
end

function Predicate(pred::LessEqualThan{names}, x::NamedTuple{names2}) where {names, names2}
    if names[1] === names2[1]
        return IsGreater(NamedTuple{(names[2],)}(Tuple(x)))
    elseif names[2] === names2[1]
        return IsLessEqual(NamedTuple{(names[1],)}(Tuple(x)))
    else
        error("Can't create conditional predicate from $pred and $x")
    end
end

"""
    GreaterThan(name1, name2)

Creates an `GreaterThan` function, which returns true on any named tuple whose field `name2`
is `isless` to the field `name1`.

See also `IsGreater` for comparing one or more columns of a table to a fixed value.
"""
struct GreaterThan{names} <: Predicate{names}
end
@inline GreaterThan(n1::Symbol, n2::Symbol) = GreaterThan{(n1, n2)}()
@inline GreaterThan(names::Tuple{Symbol,Symbol}) = GreaterThan{names}()

function (pred::GreaterThan{names})(x::NamedTuple{names2}) where {names, names2}
    if names === names2
        return isless(getproperty(x, names[2]), getproperty(x, names[1]))
    else
        return pred(Project(names)(x))
    end
end

function Predicate(pred::GreaterThan{names}, x::NamedTuple{names2}) where {names, names2}
    if names[1] === names2[1]
        return IsLessEqual(NamedTuple{(names[2],)}(Tuple(x)))
    elseif names[2] === names2[1]
        return IsGreater(NamedTuple{(names[1],)}(Tuple(x)))
    else
        error("Can't create conditional predicate from $pred and $x")
    end
end

"""
    GreaterEqualThan(name1, name2)

Creates an `GreaterEqualThan` function, which returns true on any named tuple whose field `name2`
is `isless` or `isequal` to the field `name1`.

See also `IsGreaterEqual` for comparing one or more columns of a table to a fixed value.
"""
struct GreaterEqualThan{names} <: Predicate{names}
end
@inline GreaterEqualThan(n1::Symbol, n2::Symbol) = GreaterEqualThan{(n1, n2)}()
@inline GreaterEqualThan(names::Tuple{Symbol,Symbol}) = GreaterEqualThan{names}()

function (pred::GreaterEqualThan{names})(x::NamedTuple{names2}) where {names, names2}
    if names === names2
        return !isless(getproperty(x, names[1]), getproperty(x, names[2]))
    else
        return pred(Project(names)(x))
    end
end

function Predicate(pred::GreaterEqualThan{names}, x::NamedTuple{names2}) where {names, names2}
    if names[1] === names2[1]
        return IsLess(NamedTuple{(names[2],)}(Tuple(x)))
    elseif names[2] === names2[1]
        return IsGreaterEqual(NamedTuple{(names[1],)}(Tuple(x)))
    else
        error("Can't create conditional predicate from $pred and $x")
    end
end

"""
    Within(name1, name2, distance)

Creates an `Within` function, which returns true on any named tuple whose field `name1`
is within `±distance` of `name2` (inclusive).

See also `In` for determining if one or more columns of a table is within a given `Interval`.
"""
struct Within{names, T} <: Predicate{names}
    distance::T
end
@inline Within(n1::Symbol, n2::Symbol, distance::T) where {T} = Within{(n1, n2), T}(distance)
@inline Within(names::Tuple{Symbol,Symbol}, distance::T) where {T} = Within{names, T}(distance)

function (pred::Within{names})(x::NamedTuple{names2}) where {names, names2}
    if names === names2
        x0 = getproperty(x, names[1])
        interval = (x0 - pred.distance) .. (x0 + pred.distance)
        return getproperty(x, names[2]) ∈ interval
    else
        return pred(Project(names)(x))
    end
end

function Predicate(pred::Within{names}, x::NamedTuple{names2}) where {names, names2}
    val = Tuple(x)[1]
    if names[1] === names2[1]
        return In(NamedTuple{(names[2],)}(((val - pred.distance) .. (val + pred.distance),)))
    elseif names[2] === names2[1]
        return In(NamedTuple{(names[1],)}(((val - pred.distance) .. (val + pred.distance),)))
    else
        error("Can't create conditional predicate from $pred and $x")
    end
end

# TODO: Some way of dealing with multiple indexes on one table. E.g. two sort indexes and sort-merge filter.


text:
function mc_class_prob(x::Tf, data::Td) where {Tf <: Fit, Td <: Data}
    z = get_samples(x)
    N = length(z[1].cind)
    G = size(z[1].mu, 1)

    res = zeros(Int, N, G)

    for par in x.theta
        res += class_prob(par.xi, data.Xp)
    end

    res / length(z)
end

function mc_class_prob(x::T) where T <: Fit
    z = get_samples(x)
    N = length(z[1].cind)
    G = size(z[1].mu)[1]

    res = zeros(Int, N, G)

    for par in x.theta
        for i in 1:N
            res[i, par.cind[i]] += 1
        end
    end

    res / length(z)
end

## Bayes estimator of class membership
predict_class(x::T) where T <: Fit = last.(map(findmax, eachrow(mc_class_prob(x))))

ext: jl


text:
@load_symbol _transfer_entropy :inform_transfer_entropy _libinform
@load_symbol _transfer_entropy_ensemble :inform_transfer_entropy_ensemble _libinform

function transferentropy(ys :: Vector{Int}, xs :: Vector{Int}, k;
                         b=max(2, max(maximum(xs),maximum(ys))+1))
    if length(xs) != length(ys)
        error("series must have the same shape")
    end
    ccall(_transfer_entropy, Float64, (Ptr{Int}, Ptr{Int}, Int, Int, Int), ys, xs, length(xs), b, k)
end

function transferentropy(ys :: Array{Int,2}, xs :: Array{Int,2}, k;
                         b=max(2, max(maximum(xs),maximum(ys))+1))
    if size(xs) != size(ys)
        error("series must have the same shape")
    end
    ccall(_transfer_entropy_ensemble, Float64, (Ptr{Int}, Ptr{Int}, Int, Int, Int, Int),
        ys, xs, size(xs,2), size(xs,1), b, k)
end


text:
"""model definition of dubins2d

Current model only considers the Dubins shortest path

ref:
- Dubins, Lester E (1957).
  On Curves of Minimal Length with a Constraint on Average Curvature,
  and with Prescribed Initial and Terminal Positions and Tangents.
  American Journal of Mathematics.
"""

using Dubins

struct StateDubins <: AbsState
    x::Float64
    y::Float64
    θ::Float64
end

function get_mid_status(p::StateDubins, q::StateDubins)::StateDubins
    return StateDubins((p.x + q.x) / 2, (p.y + q.y) / 2, diff_angles(q.θ, p.θ) / 2 + p.θ)
end

function dist(q1::StateDubins, q2::StateDubins)::Float64
    return norm([q1.x - q2.x, q1.y - q2.y, diff_angles(q1.θ, q2.θ) / π])
end

function get_dubins_points(
    q_i_from::State,
    q_i_to::State,
    rad::Float64;
    step_dist::Float64 = STEP_DIST,
    n_dividing::Union{Nothing,Int64} = nothing,
)::Union{Nothing,Vector{Vector{Float64}}} where {State<:AbsState}
    # same -> return one point
    q_i_from == q_i_to && return [[q_i_from.x, q_i_from.y, q_i_from.θ]]

    # obtain Dubins shortest path
    (errcode, path_i) = dubins_shortest_path(
        [q_i_from.x, q_i_from.y, q_i_from.θ],
        [q_i_to.x, q_i_to.y, q_i_to.θ],
        DUBINS_TURN_RADIUS * rad,
    )
    errcode != 0 && return nothing

    # sampling
    samples = dubins_path_sample_many(
        path_i,
        isnothing(n_dividing) ? step_dist : dubins_path_length(path_i) / n_dividing,
    )[end]

    if !isnothing(samples)
        # add goal location
        samples = vcat(samples, [[q_i_to.x, q_i_to.y, q_i_to.θ]])
    else
        # when failing to get samples -> add points of q_i_from and q_i_to
        samples = [[q_i_from.x, q_i_from.y, q_i_from.θ], [q_i_to.x, q_i_to.y, q_i_to.θ]]
    end

    return samples
end

function gen_connect(
    q::StateDubins,  # to identify type
    obstacles::Vector{CircleObstacle2D},
    rads::Vector{Float64};
    step_dist::Float64 = STEP_DIST,
    max_dist::Union{Nothing,Float64} = nothing,
)::Function

    # check: q \in C_free
    f(q::StateDubins, i::Int64)::Bool = begin
        # within field
        any(x -> (x < rads[i] || 1 - rads[i] < x), [q.x, q.y]) && return false

        # obstacles
        any([dist([q.x, q.y], o) < o.r + rads[i] for o in obstacles]) && return false

        return true
    end

    f(q_from::StateDubins, q_to::StateDubins, i::Int64)::Bool = begin
        D = dist(q_from, q_to)
        !isnothing(max_dist) && D > max_dist && return false

        P = get_dubins_points(q_from, q_to, rads[i]; step_dist = STEP_DIST)
        isnothing(P) && return false
        for p in P
            # outside
            any(e -> e < rads[i] || 1 - rads[i] < e, p[1:2]) && return false

            # obstacles
            any(o -> dist(p[1:2], o) < o.r + rads[i], obstacles) && return false
        end

        return true
    end

    return f
end

function gen_collide(
    q::StateDubins,
    rads::Vector{Float64};
    step_dist::Float64 = STEP_DIST,
)::Function

    N = length(rads)

    f(
        q_i_from::StateDubins,
        q_i_to::StateDubins,
        q_j_from::StateDubins,
        q_j_to::StateDubins,
        i::Int64,
        j::Int64,
    ) = begin
        P_i = get_dubins_points(q_i_from, q_i_to, rads[i]; step_dist = step_dist)
        P_j = get_dubins_points(q_j_from, q_j_to, rads[j]; step_dist = step_dist)
        for p_i in P_i, p_j in P_j
            dist(p_i[1:2], p_j[1:2]) < rads[i] + rads[j] && return true
        end

        return false
    end

    f(q_i::StateDubins, q_j::StateDubins, i::Int64, j::Int64) = begin
        return dist([q_i.x, q_i.y], [q_j.x, q_j.y]) < rads[i] + rads[j]
    end

    f(Q_from::Vector{Node{StateDubins}}, Q_to::Vector{Node{StateDubins}}) = begin
        for i = 1:N, j = i+1:N
            if f(Q_from[i].q, Q_to[i].q, Q_from[j].q, Q_to[j].q, i, j)
                return true
            end
        end
        return false
    end

    f(Q::Vector{Node{StateDubins}}, q_i_to::StateDubins, i::Int64) = begin
        P_i = get_dubins_points(Q[i].q, q_i_to, rads[i]; step_dist = step_dist)
        for j = 1:N
            j == i && continue
            p_j = [Q[j].q.x, Q[j].q.y]
            for p_i in P_i
                dist(p_i[1:2], p_j) < rads[i] + rads[j] && return true
            end
        end
        return false
    end

    return f
end

function gen_uniform_sampling(q::StateDubins)::Function
    () -> begin
        return StateDubins(rand(2)..., rand() * 2π)
    end
end

function plot_motion!(
    q_from::StateDubins,
    q_to::StateDubins,
    rad::Float64,
    params;
    step_dist::Float64 = STEP_DIST,
)
    q_from == q_to && return
    points = get_dubins_points(q_from, q_to, rad; step_dist = step_dist)
    p = copy(params)
    if haskey(p, :markershape)
        !isnothing(points) && scatter!([points[end][1]], [points[end][2]]; params...)
        delete!(p, :markershape)
    end
    !isnothing(points) && plot!(map(e -> e[1], points), map(e -> e[2], points); p...)
end

function plot_start_goal!(q_init::StateDubins, q_goal::StateDubins, rad::Float64, params)
    plot!([q_init.x], [q_init.y]; markershape = :hex, params...)
    plot!([q_goal.x], [q_goal.y]; markershape = :star, params...)
end

function plot_agent!(q::StateDubins, rad::Float64, color::String)
    plot_circle!(q.x, q.y, rad, color, 3.0, fillalpha = 0.1)
    p_from_x = rad * cos(q.θ) + q.x
    p_from_y = rad * sin(q.θ) + q.y
    plot!([q.x, p_from_x], [q.y, p_from_y], lw = 5, color = color, legend = nothing)
end

function gen_random_instance_StateDubins(;
    N_min::Int64 = 2,
    N_max::Int64 = 8,
    N::Int64 = rand(N_min:N_max),
    num_obs_min::Int64 = 0,
    num_obs_max::Int64 = 10,
    num_obs::Int64 = rand(num_obs_min:num_obs_max),
    rad::Float64 = 0.025,
    rad_obs::Float64 = 0.05,
    rad_min::Float64 = rad,
    rad_max::Float64 = rad,
    rad_obs_min::Float64 = rad_obs,
    rad_obs_max::Float64 = rad_obs,
    TIME_LIMIT::Float64 = 0.5,
)::Tuple{Vector{StateDubins},Vector{StateDubins},Vector{CircleObstacle2D},Vector{Float64}}
    """This model do not consider 'back' operation,
    states too close to walls/obstacles cannot be reached.
    Hence this func inserts 'additional buffer'."""

    _q = StateDubins(0, 0, 0)
    while true
        t_s = now()
        timeover() = elapsed_sec(t_s) > TIME_LIMIT

        # generate obstacles
        obstacles = gen_obstacles(2, num_obs, rad_obs_min, rad_obs_max)

        # determine rads
        rads = map(e -> rand() * (rad_max - rad_min) + rad_min, 1:N)

        connect(q::StateDubins, i::Int64) = begin
            r = rads[i] + rads[i] * DUBINS_TURN_RADIUS
            p = [q.x, q.y]
            return all(x -> r < x < 1 - r, p) && all(o -> dist(p, o) > r + o.r, obstacles)
        end

        # avoiding deadlock
        collide = gen_collide(_q, rads + rads * DUBINS_TURN_RADIUS)

        # generate starts & goals
        config_init, config_goal = gen_config_init_goal(_q, N, connect, collide, timeover)

        timeover() && continue
        return (config_init, config_goal, obstacles, rads)
    end
end

function plot_anim!(
    config_init::Vector{StateDubins},
    config_goal::Vector{StateDubins},
    obstacles::Vector{Obs} where {Obs<:Obstacle},
    ins_params...;
    solution::Union{Nothing,Vector{Vector{Node{StateDubins}}}} = nothing,
    filename::String = "tmp.gif",
    fps::Int64 = 10,
    interpolate_depth::Union{Nothing,Int64} = nothing,
    VERBOSE::Int64 = 0,
)

    if isnothing(solution)
        @warn "solution has not computed yet!"
        return
    end

    T = length(solution)
    N = length(config_init)
    anim = @animate for (t, Q) in enumerate(vcat(solution, [solution[end]]))
        VERBOSE > 0 && @printf("\rplotting t = %d / %d", t, T)
        plot_init!(StateDubins)
        plot_obs!(obstacles)
        plot_traj!(solution, ins_params...; lw = 1.0)
        plot_start_goal!(config_init, config_goal, ins_params...)

        if !isnothing(interpolate_depth) && interpolate_depth > 0 && 1 < t <= T
            depth = sum(map(k -> 2^k, 0:interpolate_depth-1)) + 2
            # compute intermediate states
            for i = 1:N
                P = get_dubins_points(
                    solution[t-1][i].q,
                    solution[t][i].q,
                    ins_params[1][i];
                    n_dividing = depth,
                )
                for (x, y, θ) in P
                    plot_agent!(
                        StateDubins(x, y, θ),
                        map(arr -> arr[i], ins_params)...,
                        get_color(i),
                    )
                end
            end
        else
            for (i, v) in enumerate(Q)
                plot_agent!(v.q, map(arr -> arr[i], ins_params)..., get_color(i))
            end
        end
    end
    VERBOSE > 0 && println()

    dirname = join(split(filename, "/")[1:end-1], "/")
    !isdir(dirname) && mkpath(dirname)

    return gif(anim, filename, fps = fps)
end


text:
import Term: load_code_and_highlight, highlight_syntax, highlight, theme

@testset "\e[34mHIGHLIGHT" begin

    @test_nothrow highlight("test 1 123 33.4 44,5 +1 -2 12 0.5, ,, ...") == "test[#90CAF9][#90CAF9] 1 [/#90CAF9][/#90CAF9][#90CAF9]123[/#90CAF9] [#90CAF9]33.4[/#90CAF9] [#90CAF9]44,5[/#90CAF9] [#90CAF9]+1 [/#90CAF9][#90CAF9]-2 [/#90CAF9][#90CAF9]12[/#90CAF9] [#90CAF9]0.5[/#90CAF9][#90CAF9], [/#90CAF9],[#90CAF9], [/#90CAF9]..."

    @test_nothrow highlight("this is ::Int64") == "this is [#d880e7]::Int64[/#d880e7]"

    @test_nothrow highlight("print", :func) == "[#fff126]print[/#fff126]"

    @test_nothrow highlight("1 + 2", :code) == "[#ffd77a]1 + 2[/#ffd77a]"

    @test_nothrow highlight("this 1 + `test`") == "this 1 + [#ffd77a]`test`[/#ffd77a]"

    @test_nothrow highlight(1)

    @test_nothrow highlight([1, 2, 3])

    @test_nothrow highlight(Int)

    @test_nothrow highlight(print)

    @test_nothrow highlight_syntax("""
This is ::Int64 my style
""") == "This is \e[1m\e[38;2;252;98;98m::\e[22m\e[39mInt64 my style\n"

    @test_nothrow highlight_syntax("""
This is ::Int64 my style
print(x + 2)
""") == "This is \e[1m\e[38;2;252;98;98m::\e[22m\e[39mInt64 my style\n\e[1m\e[38;2;255;245;157mprint\e[22m\e[39m\e[38;2;252;116;116m(\e[39mx \e[1m\e[38;2;252;98;98m+\e[22m\e[39m \e[38;2;144;202;249m2\e[39m\e[38;2;252;116;116m)\e[39m\n"

    @test_nothrow load_code_and_highlight("14_test_highlight.jl", 7)

end


text:
using Memoize
using BenchmarkTools

cd(dirname(@__FILE__))
lines = open("ten.txt", "r") do f
    readlines(f)
end
values = map(x -> parse(Int, x), lines)

push!(values, 0)
push!(values, maximum(values) + 3)

sort!(values)
d = diff(values)

j1s = count(==(1),d)
j3s = count(==(3),d)
println(j1s)
println(j3s)
println(j1s*j3s)
println()

println(values)

# using Memoize feels like cheating, but who cares ...
@memoize function chain1(values, i)
    len = length(values)
    if i == len
        return 1
    else
        value = values[i]
        j = i + 1
        acc = 0
        while (j <= len) && ((values[j] - value) <= 3)
            acc += chain1(values, j)
            j += 1
        end
        return acc
    end
end
n1 = 0
@btime n1 = chain1(values, 1)
println(n1)

# OK, let's also do it with manually memoizing the results
cache = Dict{Int64,Int64}()

function chain2(values, i)
    len = length(values)
    if i == len
        return 1
    else
        value = values[i]
        j = i + 1
        acc = 0
        while (j <= len) && ((values[j] - value) <= 3)
            if !haskey(cache, j)
                c = chain2(values, j)
                cache[j] = c
            else
                c = cache[j]
            end
            acc += c
            j += 1
        end
        return acc
    end
end

n2 = 0
@btime n2 = chain2(values,1)
println(n2)


text:
module ParameterInference
using Catch22
using NonstationaryProcesses

# windows:      x -> Array of windows
# features:     TS array -> feature array
# projections:  feature array -> low dim projection array
# estimates:    projection array --> parameter estimates

Base.@kwdef struct Inference
    # These you can set
    timeseries
    windows = slidingwindow
    features = catch24
    baseline = identity
    filter = identity #nonanrows∘noconstantrows∘noinfrows
    normalisation = identity#standardise
    dimensionalityReduction = principalcomponents
    parameters = timeseries .+ NaN
    # These you should leave to calculate
    windowedTimeseries = (windows)(timeseries)
    windowTimes = (windows)(1:length(timeseries))
    windowEdges = [windowTimes[1, :][:]..., windowTimes[end]]
    windowCentres = mean(windowTimes, dims=1)[:]
    F = (!(features isa AbstractFeatureSet) && !(features isa Vector{Function}) ? features : features(windowedTimeseries)) # Maybe you supplied some pre-computed features?
    F̂ = (filter∘baseline∘filter∘normalisation∘filter)(F) #! Too many filters?
    model = project(Array(F̂), dimensionalityReduction)
    F′ = embed(model, Array(F̂))
    estimates = embed(model, Array(F̂), [1])
    corrtype = corspearman
    ρ = corrtype(mean((windows)(parameters), dims=1)[:], estimates)
end
export Inference

include("Baseline.jl")
include("FeatureClustering.jl")
include("FeatureRepresentations.jl")
include("LowDimensionalProjections.jl")
include("Normalisation.jl")
include("Plotting.jl")
include("Windows.jl")
include("DistributionTests.jl")
include("Information.jl")
include("BaselineFilter.jl")
include("BoxCox.jl")
include("TestBaselines.jl")

Base.identity(x, y) = identity

function infer(x::AbstractVector; kwargs...)
    Inference(timeseries=x; kwargs...)
end

function infer(P::NonstationaryProcesses.Process, dim::Int=1; parameters::Int=1, kwargs...)
    infer(timeseries(P, dim), parameters=NonstationaryProcesses.parameterseries(P, p=parameters); kwargs...)
end
export infer

end


text:
@testset "How to stop NLP" begin
  ###############################################################################
  #
  # The Stopping structure eases the implementation of algorithms and the
  # stopping criterion.
  #
  # We illustrate here the basic features of NLPStopping, which is a
  # specialized version of the Stopping to the case where:
  # pb is an AbstractNLPModel
  # state shares the structure of the NLPAtX.
  #
  # NLPModels is a package to handle non-linear (constrained) optimization.
  # NLPStopping is following this approach.
  #
  ###############################################################################
  #using Test, NLPModels, Stopping

  #We first create a toy problem
  f(x) = sum(x .^ 2)
  x0 = zeros(5)
  nlp = ADNLPModel(f, x0)
  #Warning: see https://github.com/JuliaSmoothOptimizers/NLPModels.jl/blob/master/src/autodiff_model.jl
  #for the proper way of defining an ADNLPModel
  meta = NLPModelMeta(5, x0 = x0, lvar = zeros(5), uvar = Inf * ones(5))
  nlp2 = ADNLPModel(meta, Counters(), ADNLPModels.ForwardDiffAD(5, f, x0), f, x -> [])
  nlp_at_x = NLPAtX(x0)
  x1 = ones(5)

  ###############################################################################
  #1) Initialize the NLPStopping.
  #The specificity here is that the NLPStopping requires another mandatory input:
  #optimality_check which is the function later used to compute the score.
  #Recall that the score is then tested at 0, to declare optimality.
  #Stopping provides a default KKT function.
  stop_nlp = NLPStopping(nlp, nlp_at_x, optimality_check = (x, y) -> KKT(x, y))

  #Another approach is to use the lazy way:
  stop_nlp_lazy = NLPStopping(nlp2) #use nlp.meta.x0 as initial point
  @test stop_nlp_lazy.current_state.x == nlp2.meta.x0

  ###############################################################################
  #2) Fill in
  #Before calling start! and stop! one should fill in current information in the
  #State. -> the optimality_check then exploits this knowledge.
  #As seen before, we by hand use update_and_start and update_and_stop.
  #Another way is to call the fill_in! function:
  fill_in!(stop_nlp, x1, matrix_info = false)
  @test stop_nlp.current_state.x == x1
  @test stop_nlp.current_state.fx == 5.0
  @test stop_nlp.current_state.Hx == zeros(0, 0)
  #Note that since there are no constraints, c(x) and J(x) are not called:
  @test stop_nlp.current_state.Jx == zeros(0, 0)
  @test stop_nlp.current_state.cx == zeros(0)
  #Since there are no bounds on x, the Lagrange multiplier is not updated:
  @test stop_nlp.current_state.mu == zeros(0)

  #would give Hx if matrix_info = true
  fill_in!(stop_nlp_lazy, x1)
  @test stop_nlp_lazy.current_state.Hx != zeros(0, 0)
  #stop_nlp_lazy.pb has bounds, so mu is a vector of size x
  @test size(x0) == size(stop_nlp_lazy.current_state.mu)

  ###############################################################################
  #3) Evaluations
  #Another particularity is that the NLPModels has a counter keeping track of
  #the evaluations of each function.
  #Similarly the NLPStopping has a dictionary keeping all the maximum number of
  #evaluations:
  @test typeof(stop_nlp.meta.max_cntrs) <: Dict
  #For instance the limit in evaluations of objective and gradient:
  @test stop_nlp.meta.max_cntrs[:neval_obj] == typemax(Int)
  @test stop_nlp.meta.max_cntrs[:neval_grad] == typemax(Int)

  #Limit can be set using init_max_counters function:
  stop_nlp.meta.max_cntrs = init_max_counters(obj = 3, grad = 0, hess = 0)
  @test stop_nlp.meta.max_cntrs[:neval_obj] == 3
  @test stop_nlp.meta.max_cntrs[:neval_grad] == 0

  OK = update_and_stop!(stop_nlp, evals = stop_nlp.pb.counters)
  @test OK == true
  @test stop_nlp.meta.resources == true
  @test status(stop_nlp) == :EvaluationLimit

  ###############################################################################
  #4) Unbounded problem
  #An additional feature of the NLPStopping is to provide an _unbounded_problem_check
  #whenever \|c(x)\| or -f(x) become too large.
  stop_nlp.meta.unbounded_threshold = -6.0 #by default 1.0e50
  stop!(stop_nlp)
  @test stop_nlp.meta.unbounded_pb == true
  @test stop_nlp.current_state.fx > stop_nlp.meta.unbounded_threshold
  @test stop_nlp.meta.resources == true #still true as the state has not changed

  ###############################################################################
  #An advanced feature is the possibility to send keywords to optimality_check:
  optimality_fct_test = (x, y; a = 1.0) -> a

  #In this case, the optimality_check function used to compute the score may
  #depend on a parameter (algorithm-dependent for instance)
  stop_nlp_2 = NLPStopping(nlp, nlp_at_x, optimality_check = optimality_fct_test)
  fill_in!(stop_nlp_2, x0)

  OK = stop!(stop_nlp_2, a = 0.0)
  @test OK == true
  @test stop_nlp_2.meta.optimal == true

  #However, note that the same cannot be achieved with update_and_stop!:
  reinit!(stop_nlp_2)
  OK = update_and_stop!(stop_nlp_2, a = 0.0)
  @test OK == false
end


text:
using .Empirikos
using Zygote
using LinearFractional
using StatsFuns

abstract type AbstractRegressionDiscontinuityTarget end
abstract type TargetedRegressionDiscontinuityTarget <: AbstractRegressionDiscontinuityTarget end

struct ConstantTarget <: AbstractRegressionDiscontinuityTarget end

Base.@kwdef struct SharpRegressionDiscontinuityTarget{T,S} <: TargetedRegressionDiscontinuityTarget
    cutoff::T
    denom_est::S = 1.0
end

function (sharp_target::SharpRegressionDiscontinuityTarget)(u)
    likelihood( sharp_target.cutoff, u) / sharp_target.denom_est
end

function update_target(sharp_target::SharpRegressionDiscontinuityTarget, prior)
    denom_est = pdf(prior, sharp_target.cutoff)
    @set sharp_target.denom_est = denom_est
end

function Base.getproperty(obj::RunningVariable{<:EBayesSample}, sym::Symbol)
    if sym === :ZsC
        return response.(obj.Zs) .- obj.cutoff
    else
        return getfield(obj, sym)
    end
end

function StatsDiscretizations.unwrap(Z::EBayesSample)
    response(Z)
end

function StatsDiscretizations.wrap(Z::EBayesSample, val)
    @set Z.Z = val
end

export NoiseInducedRandomization

struct NoiseInducedRandomization{RDT, T,D,F,G,P,V,S} <: RegressionDiscontinuity.SharpRD
    response_lower_bound::T
    response_upper_bound::T
    solver
    convexclass::G
    plugin_G::P
    flocalization::F
    target::RDT
    discretizer::D
    bias_opt_multiplier::S
    σ_squared::V
    τ_range::S
    skip_bias::Bool
    maxbias_opt_bisection::Int
    α::S
end

function NoiseInducedRandomization(;
    response_lower_bound = 0.0,
    response_upper_bound = 1.0,
    convexclass = :default,
    solver,
    plugin_G = :default,
    flocalization = :default,
    target = ConstantTarget(),
    discretizer = :default,
    bias_opt_multiplier = 1.0,
    σ_squared = abs2(response_upper_bound - response_lower_bound),
    τ_range = 2*(response_upper_bound - response_lower_bound),
    skip_bias=false,
    maxbias_opt_bisection = 50,
    α = 0.05
    )
    NoiseInducedRandomization(
        response_lower_bound,
        response_upper_bound,
        solver,
        convexclass,
        plugin_G,
        flocalization,
        target,
        discretizer,
        bias_opt_multiplier,
        σ_squared,
        τ_range,
        skip_bias,
        maxbias_opt_bisection,
        α
    )
end

Base.@kwdef struct NoiseInducedRandomizationWeights{EB, Γ₊, Γ₋, H₊, H₋, M, T}
    γ₊::Γ₊
    γ₋::Γ₋
    h₊::H₊
    h₋::H₋
    eb_sample::EB
    model::M
    marginal_probs_treated::T
    marginal_probs_untreated::T
end

function nir_default_convexclass(Zs::AbstractVector{<:BinomialSample})
    DiscretePriorClass(range(1e-4, stop=1-1e-4, length=400))
end

function nir_default_discretizer(ZsR::RunningVariable{<:BinomialSample})
    cutoff = Int(ZsR.cutoff)
    Zs = ZsR.Zs
    skedasticity(Zs) === Empirikos.Homoskedastic() || throw("Only homoskedastic case supported.")
    K = ntrials(Zs[1])
    ℓ = 0
    u = K

    discr_all = FiniteSupportDiscretizer(0:1:K)
    if ZsR.treated === :≥
        discr_untreated = FiniteSupportDiscretizer(ℓ:1:(cutoff-1))
        discr_treated = FiniteSupportDiscretizer(cutoff:1:u)
    elseif ZsR.treated === :>
        discr_untreated = FiniteSupportDiscretizer(ℓ:1:cutoff)
        discr_treated = FiniteSupportDiscretizer((cutoff+1):1:u)
    elseif ZsR.treated === :≤
        discr_treated = FiniteSupportDiscretizer(ℓ:1:cutoff)
        discr_untreated = FiniteSupportDiscretizer((cutoff+1):1:u)
    elseif ZsR.treated === :<
        discr_treated = FiniteSupportDiscretizer(ℓ:1:(cutoff-1))
        discr_untreated = FiniteSupportDiscretizer(cutoff:1:u)
    else
        throw(ArgumentError("RunningVariable:treated can only be one of :>, :<, :≥, :≤"))
    end

    (untreated = discr_untreated,
     treated = discr_treated,
     all = discr_all)
end

function nir_default_convexclass(Zs::AbstractVector{<:Empirikos.AbstractNormalSample})
    z_min, z_max = extrema(response.(Zs))
    DiscretePriorClass(range(z_min, stop=z_max, length=500))
end

function nir_default_discretizer(ZsR::RunningVariable{<:Empirikos.AbstractNormalSample})
    cutoff = ZsR.cutoff
    Zs = ZsR.Zs
    skedasticity(Zs) === Empirikos.Homoskedastic() || throw("Only homoskedastic case supported.")
    ν = Empirikos.std(Zs[1])
    # TODO: add cases depending on >= who is treated
    ℓ, u = extrema(response.(Zs))
    ℓ = ℓ - ν/50
    u = u + ν/50

    _range_left = collect(RangeHelpers.range(start=below(ℓ), stop=cutoff, step = ν/20 ))
    _range_right = collect(RangeHelpers.range(start=cutoff, stop=above(u), step = ν/20))

    _range_all = [_range_left[1:(end-1)]; cutoff; _range_right[2:end]]

    if ZsR.treated === :≥
        discr_untreated = BoundedIntervalDiscretizer{:closed,:open}(_range_left)
        discr_treated = BoundedIntervalDiscretizer{:closed,:open}(_range_right)
        discr_all = RealLineDiscretizer{:closed,:open}(_range_all)
    elseif ZsR.treated === :>
        discr_untreated = BoundedIntervalDiscretizer{:open,:closed}(_range_left)
        discr_treated = BoundedIntervalDiscretizer{:open,:closed}(_range_right)
        discr_all = RealLineDiscretizer{:open,:closed}(_range_all)
    elseif ZsR.treated === :≤
        discr_treated = BoundedIntervalDiscretizer{:open,:closed}(_range_left)
        discr_untreated = BoundedIntervalDiscretizer{:open,:closed}(_range_right)
        discr_all = RealLineDiscretizer{:open,:closed}(_range_all)
    elseif ZsR.treated === :<
        discr_treated = BoundedIntervalDiscretizer{:closed,:open}(_range_left)
        discr_untreated = BoundedIntervalDiscretizer{:closed,:open}(_range_right)
        discr_all = RealLineDiscretizer{:closed,:open}(_range_all)
    else
        throw(ArgumentError("RunningVariable:treated can only be one of :>, :<, :≥, :≤"))
    end

    (untreated = discr_untreated,
     treated = discr_treated,
     all = discr_all)
end

function initialize(nir::NoiseInducedRandomization, ZsR, Ys)
    Zs = ZsR.Zs
    cutoff = ZsR.cutoff

    n = nobs(Zs)
    if nir.flocalization === :default
        α_min = min(n^(-1/4),nir.α)
        nir = @set nir.flocalization = DvoretzkyKieferWolfowitz(α_min)
    end

    if nir.convexclass === :default
        #TODO: customize by sample type
        nir = @set nir.convexclass = nir_default_convexclass(Zs)
    end

    if nir.plugin_G === :default
        nir = @set nir.plugin_G = NPMLE(; solver = nir.solver, convexclass = nir.convexclass)
    end

    if nir.discretizer === :default
        # check for Homoskedastic
        nir = @set nir.discretizer = nir_default_discretizer(ZsR)
    end

    if nir.σ_squared === :lm
        lm_fit = fit(LinearModel, @formula(Ys~response(Zs)*Ws), RDData(Ys, ZsR))
        σ_squared = mean(abs2, residuals(lm_fit))
        nir = @set nir.σ_squared = σ_squared
    end
    nir
end

function nir_h(γ, Zs_levels, disc = Zs_levels)
    function h(u)
        dot(likelihood.(Zs_levels, u), γ.(disc))
    end
end

function add_bias_constraint!(model, nir::NoiseInducedRandomization{<:ConstantTarget}, h₊, h₋, t)
    M = nir.bias_opt_multiplier * (nir.response_upper_bound - nir.response_lower_bound)
    for u in Distributions.support(nir.convexclass)
        @constraint(model,  M*(h₊(u) - h₋(u)) <=  t)
        @constraint(model,  M*(h₊(u) - h₋(u)) >= -t)
    end
end

function add_bias_constraint!(model, nir::NoiseInducedRandomization{<:TargetedRegressionDiscontinuityTarget}, h₊, h₋, t)
    target = nir.target
    M = nir.bias_opt_multiplier * (nir.response_upper_bound - nir.response_lower_bound)
    τ_range = nir.τ_range/2

    @variable(model, t1)
    @variable(model, t2)

    for u in Distributions.support(nir.convexclass)
        @constraint(model,  M*(h₊(u) - h₋(u)) <=  t1)
        @constraint(model,  M*(h₊(u) - h₋(u)) >= -t1)
        @constraint(model,  τ_range*(h₊(u) - target(u)) <=  t2)
        @constraint(model,  τ_range*(h₊(u) - target(u)) >= -t2)
        @constraint(model,  τ_range*(h₋(u) - target(u)) <=  t2)
        @constraint(model,  τ_range*(h₋(u) - target(u)) >= -t2)
    end
    @constraint(model, t1 + t2 <= t)
end

function nir_weights_quadprog(nir::NoiseInducedRandomization, Zs)
    model = Model(nir.solver)
    σ_squared = nir.σ_squared
    n = nobs(Zs)
    eb_sample = skedasticity(Zs) === Empirikos.Homoskedastic() ? Zs[1] : throw("Only implemented for homoskedastic samples.")

    γ₊ = StatsDiscretizations.add_discretized_function!(model, nir.discretizer.treated)
    γ₋ = StatsDiscretizations.add_discretized_function!(model, nir.discretizer.untreated)

    Zs_levels_treated = Empirikos.set_response.(eb_sample, nir.discretizer.treated)
    Zs_levels_untreated = Empirikos.set_response.(eb_sample, nir.discretizer.untreated)

    h₊ = nir_h(γ₊, Zs_levels_treated, nir.discretizer.treated)
    h₋ = nir_h(γ₋, Zs_levels_untreated, nir.discretizer.untreated)

    @variable(model, t)

    add_bias_constraint!(model, nir, h₊, h₋, t)

    marginal_probs_treated = pdf.(nir.plugin_G, Zs_levels_treated)
    marginal_probs_untreated = pdf.(nir.plugin_G, Zs_levels_untreated)

    #return (marginal_probs_treated,marginal_probs_untreated, nir.plugin_G, Zs_levels_treated, Zs_levels_untreated)
    @constraint(model, dot(marginal_probs_treated, γ₊.(Zs_levels_treated)) ==  1)
    @constraint(model, dot(marginal_probs_untreated, γ₋.(Zs_levels_untreated)) ==  1)

    @variable(model, s)
    @constraint(model, [s;
                        t;
                        sqrt(σ_squared/n) .* γ₊.(nir.discretizer.treated) .* sqrt.(marginal_probs_treated);
                        sqrt(σ_squared/n) .* γ₋.(nir.discretizer.untreated) .* sqrt.(marginal_probs_untreated)
                    ]  ∈ SecondOrderCone())

    @objective(model, Min, s)
    optimize!(model)

    γ₊_fun = JuMP.value(γ₊, z->Empirikos.set_response(eb_sample, z))
    γ₊_fun = @set γ₊_fun.discretizer = nir.discretizer.all
    γ₋_fun = JuMP.value(γ₋, z->Empirikos.set_response(eb_sample, z))
    γ₋_fun = @set γ₋_fun.discretizer = nir.discretizer.all

    #return γ₊_fun, Zs_levels_treated
    h₊ = nir_h(γ₊_fun, Zs_levels_treated)
    h₋ = nir_h(γ₋_fun, Zs_levels_untreated)

    NoiseInducedRandomizationWeights(;γ₊=γ₊_fun, γ₋=γ₋_fun,
        h₊ = h₊, h₋ = h₋,
        model=model,
        eb_sample = eb_sample,
        marginal_probs_treated=marginal_probs_treated,
        marginal_probs_untreated=marginal_probs_untreated
    )
end

function nir_maxbias(nir::NoiseInducedRandomization{<:ConstantTarget}, γs, Zs)
    M = nir.response_upper_bound - nir.response_lower_bound

    bias_model = LinearFractionalModel(nir.solver)
    convexclass = nir.convexclass
    us = support(convexclass)

    hplus  = γs.h₊.(us)
    hminus = γs.h₋.(us)

    π = Empirikos.prior_variable!(bias_model, convexclass)
    _nparams = Empirikos.nparams(convexclass)

    @variable(bias_model, α[1:_nparams] >= 0)

    dkw_band = StatsBase.fit(nir.flocalization, Zs)
    Empirikos.flocalization_constraint!(bias_model, dkw_band, π)

    @constraint(bias_model.model, dot(π.finite_param, hplus) == 1)

    @objective(bias_model.model, JuMP.MOI.MAX_SENSE, dot(π.finite_param, hminus))
    optimize!(bias_model)

    _ξ_max = JuMP.objective_value(bias_model)

    @objective(bias_model.model, JuMP.MOI.MIN_SENSE, dot(π.finite_param, hminus))
    optimize!(bias_model)
    _ξ_min = JuMP.objective_value(bias_model)

    @constraint(bias_model.model, α .<= M.*π.finite_param)

    ξ = _ξ_min
    ξs = range(_ξ_min; stop=_ξ_max, length=nir.maxbias_opt_bisection)
    max_vals = zeros(length(ξs))
    #min_vals = zeros(length(ξs))

    @constraint(bias_model.model, parametric_constraint, dot(π.finite_param, hminus) == _ξ_min)

    for (i, ξ) in enumerate(ξs)
        JuMP.set_normalized_rhs(parametric_constraint, ξ)
        @objective(bias_model.model, JuMP.MOI.MAX_SENSE, dot(α, hplus) - dot(α, hminus)/ξ)
        optimize!(bias_model)
        max_vals[i] = JuMP.objective_value(bias_model)

        #@objective(bias_model.model, JuMP.MOI.MIN_SENSE, dot(α, hplus) - dot(α, hminus)/ξ)
        #optimize!(bias_model)
        #min_vals[i] = JuMP.objective_value(bias_model)
    end

    ξ = ξs[argmax(max_vals)]
    JuMP.set_normalized_rhs(parametric_constraint, ξ)
    @objective(bias_model.model, JuMP.MOI.MAX_SENSE, dot(α, hplus) - dot(α, hminus)/ξ)
    optimize!(bias_model)
    maxbias = JuMP.objective_value(bias_model)

    (maxbias=maxbias, bias_model=bias_model, ξ=ξ, ξ_min = _ξ_min, ξ_max = _ξ_max,
     ξs = ξs, grid_vals = max_vals)
end

function nir_maxbias(nir::NoiseInducedRandomization{<:TargetedRegressionDiscontinuityTarget}, γs, Zs)

    M = nir.response_upper_bound - nir.response_lower_bound
    τ_range = nir.τ_range

    bias_model = LinearFractionalModel(nir.solver)
    convexclass = nir.convexclass
    us = support(convexclass)

    hplus  = γs.h₊.(us)
    hminus = γs.h₋.(us)
    ws = nir.target.(us)

    π = Empirikos.prior_variable!(bias_model, convexclass)
    _nparams = Empirikos.nparams(convexclass)

    @variable(bias_model, α[1:_nparams] >= 0)
    @variable(bias_model, dΤ[1:_nparams] >= 0)

    dkw_band = StatsBase.fit(nir.flocalization, Zs)
    Empirikos.flocalization_constraint!(bias_model, dkw_band, π)

    @constraint(bias_model.model, dot(π.finite_param, hplus) == 1)

    @objective(bias_model.model, JuMP.MOI.MAX_SENSE, dot(π.finite_param, hminus))
    optimize!(bias_model)

    _ξ_max = JuMP.objective_value(bias_model)

    @objective(bias_model.model, JuMP.MOI.MIN_SENSE, dot(π.finite_param, hminus))
    optimize!(bias_model)
    _ξ_min = JuMP.objective_value(bias_model)

    @constraint(bias_model.model, α .<= M .*    π.finite_param)
    @constraint(bias_model.model, dΤ .<= τ_range .* π.finite_param)

    ξ = _ξ_min
    ξs = range(_ξ_min; stop=_ξ_max, length=nir.maxbias_opt_bisection)

    @constraint(bias_model.model, parametric_constraint, dot(π.finite_param, hminus) == ξ)

    max_ws = zeros(length(ξs))
    min_ws = zeros(length(ξs))

    for (i, ξ) in enumerate(ξs)
        JuMP.set_normalized_rhs(parametric_constraint, ξ)
        @objective(bias_model.model, JuMP.MOI.MAX_SENSE, dot(π.finite_param, ws))
        optimize!(bias_model)
        max_ws[i] = JuMP.objective_value(bias_model)

        @objective(bias_model.model, JuMP.MOI.MIN_SENSE, dot(π.finite_param, ws))
        optimize!(bias_model)
        min_ws[i] = JuMP.objective_value(bias_model)
    end

    @constraint(bias_model.model, weight_constraint, dot(π.finite_param, ws) ==  max_ws[1])

    argmax_ws = zeros(length(ξs))
    maxvals_ws = zeros(length(ξs))

    #TODOs: last resolve.
    for (i, ξ) in enumerate(ξs)#enumerate(ξs[2:(end-1)])
        JuMP.set_normalized_rhs(parametric_constraint, ξ)
        _w_list = RangeHelpers.range(start = min_ws[i], stop=max(min_ws[i], max_ws[i]), step=below(0.02 * _ξ_min))
        tmp_vals = zeros(length(_w_list))
        for (j,w) in enumerate(_w_list)
            JuMP.set_normalized_rhs(weight_constraint, w)
            @objective(
                bias_model.model,
                JuMP.MOI.MAX_SENSE,
                dot(α, hplus) - dot(α, hminus)/ξ + dot(dΤ, hplus) - dot(dΤ, ws) / w
            )
            optimize!(bias_model)
            tmp_vals[j] = JuMP.objective_value(bias_model)
        end
        argmax_w_idx = argmax(tmp_vals)
        argmax_ws[i] = _w_list[argmax_w_idx]
        maxvals_ws[i] = tmp_vals[argmax_w_idx]
    end

    maxbias = maximum(maxvals_ws)

    (maxbias=maxbias, bias_model=bias_model, ξ=ξ, ξ_min = _ξ_min, ξ_max = _ξ_max,
     ξs = ξs, max_ws=max_ws, min_ws=min_ws, argmax_ws=argmax_ws, grid_vals = maxvals_ws)
end

Base.@kwdef struct FittedNoiseInducedRandomization
    coeftable
    confint
    γs
    maxbias_fit
    nir
    model = nothing
end

function Base.show(io::IO, nir::FittedNoiseInducedRandomization)
    Base.println(io::IO, "RD analysis with Noise Induced Randomization (NIR)")
    Base.show(io, nir.coeftable)
end

function StatsBase.fit(nir::NoiseInducedRandomization, ZsR::RunningVariable, Ys)
    nir = initialize(nir, ZsR, Ys)
    convexclass = nir.convexclass
    Zs = ZsR.Zs
    nir = @set nir.plugin_G = StatsBase.fit(nir.plugin_G, summarize(nir.discretizer.all.(Zs)))

    if !isa(nir.target, ConstantTarget)
        nir = @set nir.target = update_target(nir.target, nir.plugin_G)
    end

    γs = nir_weights_quadprog(nir, Zs)
    nir_inference(nir, ZsR, Ys, γs)
end

function nir_sensitivity(nir_prev::FittedNoiseInducedRandomization, ZsR::RunningVariable, Ys,
                                 τ_range)
    nir = nir_prev.nir
    nir = @set nir.τ_range = τ_range

    nir_inference(nir, ZsR, Ys, nir_prev.γs)
end

function nir_inference(nir::NoiseInducedRandomization, ZsR::RunningVariable, Ys, γs)

    Zs = ZsR.Zs

    #return γs
    γ₊_denom = sum( γs.γ₊.(Zs,0))
    γ₋_denom = sum( γs.γ₋.(Zs,0))

    μ̂γ₊ = dot( γs.γ₊.(Zs,0), Ys) / γ₊_denom
    μ̂γ₋ = dot( γs.γ₋.(Zs,0), Ys) / γ₋_denom
    τ̂ = μ̂γ₊ - μ̂γ₋

    σ̂_squared =   dot( γs.γ₊.(Zs,0).^2, (Ys .- μ̂γ₊).^2)/γ₊_denom^2 +
                  dot( γs.γ₋.(Zs,0).^2, (Ys .- μ̂γ₋).^2)/γ₋_denom^2

    se = sqrt(σ̂_squared)

    if !nir.skip_bias
        maxbias_fit =  nir_maxbias(nir, γs, Zs)
        maxbias = maxbias_fit.maxbias
    else
        maxbias_fit = nothing
        maxbias = 0.0
    end

    level = 1-nir.α
    ci_halfwidth = bias_adjusted_gaussian_ci(se; maxbias = maxbias, level = level)
    ci = [τ̂ - ci_halfwidth, τ̂ + ci_halfwidth]
    levstr = isinteger(level * 100) ? string(Integer(level * 100)) : string(level * 100)
    res = [τ̂ se maxbias first(ci) last(ci) ci_halfwidth]
    colnms = ["τ̂"; "se"; "max bias"; "Lower $levstr%"; "Upper $levstr%"; "CI halfwidth"]
    rownms = ["Weighted RD estimand"]
    coeftbl = CoefTable(res, colnms, rownms, 0, 0)

    ci = Empirikos.BiasVarianceConfidenceInterval(estimate = τ̂, α = nir.α, se = se, maxbias=maxbias)

    FittedNoiseInducedRandomization(coeftable = coeftbl, confint=ci,
        γs = γs, maxbias_fit = maxbias_fit,
        nir = nir)
end

@recipe function f(rdd_weights::RegressionDiscontinuity.NoiseInducedRandomizationWeights{<:BinomialSample})
    γ₋ = rdd_weights.γ₋
    γ₊ = rdd_weights.γ₊
    ys = [.- collect(γ₋.dictionary); collect(γ₊.dictionary)]
    xs = response.([collect(keys(γ₋.dictionary)); collect(keys(γ₊.dictionary))])

    @series begin
        seriestype := :scatter
        seriescolor := :lightblue
        markerstrokewidth := 0
        markersize :=5
        label := ""
        xs, ys
    end

    label --> nothing
    seriestype --> :sticks
    linecolor --> :lightblue
    linewidth --> 2
    yguide --> L"\gamma(z)"
    xguide --> L"z"
    xs, ys
end

### Curvature Bound ###

Base.@kwdef struct NoiseInducedRandomizationCurvature{EB, F, G<:DiscretePriorClass, T, S, TT} <: WorstCaseCurvatureSelector
    response_lower_bound::T = 0.0
    response_upper_bound::T = 1.0
    density_lower_bound::S
    ebayes_sample::EB
    flocalization::F = nothing
    solver
    convexclass::G =  DiscretePriorClass(-2.5:0.01:2.5)
    n_grid_f_prime::Int = 200
    width_grid_f_double_prime::TT = 0.02
end

Base.@kwdef struct FittedNoiseInducedRandomizationCurvature{S,T, C<:NoiseInducedRandomizationCurvature}
    max_curvature::S
    us::T
    πs::T
    α0s::T
    f
    f′
    f′′
    nir_curvature::C
end

function curvature(nir_curvature::FittedNoiseInducedRandomizationCurvature)
    nir_curvature.max_curvature
end

function curvature(nir_curvature::NoiseInducedRandomizationCurvature)
    curvature(fit( nir_curvature ))
end

function fit(nir_curvature::NoiseInducedRandomizationCurvature)
    convexclass = nir_curvature.convexclass
    us = support(convexclass)
    _nparams = Empirikos.nparams(convexclass)
    M = nir_curvature.response_upper_bound - nir_curvature.response_lower_bound
    Z = nir_curvature.ebayes_sample
    model = Model(nir_curvature.solver)

    @variable(model, π[1:_nparams] >= 0)
    @variable(model, t >= 0)
    @constraint(model, sum(π) == t)
    @variable(model, H[1:_nparams] >= 0)
    @constraint(model, H .<=  M.* π)

    if isa(Z, BinomialSample)
        @show "binom"
        n_trials = ntrials(Z)
        function quasibinom_pdf(n, p, m)
            val = betalogpdf(m + 1, n - m + 1, p) - log(n + 1)
            exp(val)
        end
        f_vec = [z -> quasibinom_pdf(n_trials, u, z) for u in us]
    else
        f_vec = [z -> likelihood(Empirikos.set_response(nir_curvature.ebayes_sample, z), u) for u in us]
    end
    cutoff = Float64(response(Z))
    f_vec_0 = [f(cutoff) for f in f_vec]
    f_vec_prime_0 = [gradient(f, cutoff)[1] for f in f_vec]
    f_vec_double_prime_0 = [hessian(f, cutoff) for f in f_vec]

    f = @expression(model, dot(π, f_vec_0))
    f′ = @expression(model, dot(π, f_vec_prime_0))
    f′′ = @expression(model, dot(π, f_vec_double_prime_0))

    α = @expression(model, dot(H, f_vec_0))
    α′ = @expression(model, dot(H, f_vec_prime_0))
    α′′ = @expression(model, dot(H, f_vec_double_prime_0))

    f_min = nir_curvature.density_lower_bound
    @constraint(model, f >= f_min*t)
    @constraint(model, f == 1)

    @objective(model, JuMP.MOI.MAX_SENSE,  f′)
    optimize!(model)
    _ξ_max = JuMP.value(f′)
    ξs = range(0.0; stop = _ξ_max, length = nir_curvature.n_grid_f_prime)

    @constraint(model, parametric_constraint, f′ == ξs[1])

    # range for f' and f''
    max_ws = zeros(length(ξs))
    min_ws = zeros(length(ξs))

    for (i, ξ) in enumerate(ξs)
        JuMP.set_normalized_rhs(parametric_constraint, ξ)
        @objective(model, JuMP.MOI.MAX_SENSE, f′′)
        optimize!(model)
        max_ws[i] = JuMP.value(f′′)

        @objective(model, JuMP.MOI.MIN_SENSE, f′′)
        optimize!(model)
        min_ws[i] = JuMP.value(f′′)
    end

    @constraint(model, weight_constraint, f′′ ==  max_ws[1])

    argmax_ws = zeros(length(ξs))
    maxvals_ws = zeros(length(ξs))

    for (i, ξ) in enumerate(ξs)
        JuMP.set_normalized_rhs(parametric_constraint, ξ)
        _w_list = RangeHelpers.range(start = min_ws[i], stop=max(min_ws[i], max_ws[i]),
                     step=below(nir_curvature.width_grid_f_double_prime))
        tmp_vals = zeros(length(_w_list))
        for (j,w) in enumerate(_w_list)
            JuMP.set_normalized_rhs(weight_constraint, w)
            @objective(
                model,
                JuMP.MOI.MAX_SENSE,
                α′′ - 2*α′*ξ - α*w + 2α*ξ^2
            )
            optimize!(model)
            tmp_vals[j] = JuMP.objective_value(model)
        end
        argmax_w_idx = argmax(tmp_vals)
        argmax_ws[i] = _w_list[argmax_w_idx]
        maxvals_ws[i] = tmp_vals[argmax_w_idx]
    end

    max_curvature = maximum(maxvals_ws)
    _idx_ξ_worst_case = argmax(maxvals_ws)
    _ξ_worst_case= ξs[_idx_ξ_worst_case]
    w_worst_case = argmax_ws[_idx_ξ_worst_case]
    JuMP.set_normalized_rhs(parametric_constraint, _ξ_worst_case)
    JuMP.set_normalized_rhs(weight_constraint, w_worst_case)
    @objective(
        model,
        JuMP.MOI.MAX_SENSE,
        α′′ - 2*α′*_ξ_worst_case - α*w_worst_case + 2α*_ξ_worst_case^2
    )
    optimize!(model)

    π_worst = JuMP.value.(π)
    H_worst = JuMP.value.(H)

    πs = π_worst ./ sum(π_worst)
    α0s = ifelse.( iszero.(π_worst), zero(eltype(π_worst)) , H_worst ./ π_worst)

    myf = z -> dot(H_worst, [f(z) for f in f_vec])/dot(π_worst, [f(z) for f in f_vec])
    myf′ = z -> gradient(myf, z)[1]
    myf′′ = z -> hessian(myf, z)

    FittedNoiseInducedRandomizationCurvature(;
        max_curvature = myf′′(cutoff),
        us=collect(us),
        πs=πs,
        α0s=α0s,
        f=myf,
        f′=myf′,
        f′′=myf′′,
        nir_curvature = nir_curvature
    )
end


text:
const SetCounter = Dict{Symbol, Vector{Int}}
const SetDeckCounter = Dict{Symbol,Array{Int64,2}}
const SetRarityCounter = Dict{Symbol,Array{Int64,2}}

abstract type RarityGrouping end

function Base.:getindex(x::RarityGrouping, rarity::Int)
    rarity == 1 && return x.commons
    rarity == 2 && return x.rares
    rarity == 3 && return x.epics
    return x.legendaries
end

Base.:start(iter::RarityGrouping) = 0
function Base.:next(iter::RarityGrouping, state)
    state == 0 && return (iter.commons, 1)
    state == 1 && return (iter.rares, 2)
    state == 2 && return (iter.epics, 3)
    return (iter.legendaries, 4)
end
Base.:done(iter::RarityGrouping, state) = state > 3

struct Collection <: RarityGrouping
    commons::Dict{Symbol,Vector{Int}}
    rares::Dict{Symbol,Vector{Int}}
    epics::Dict{Symbol,Vector{Int}}
    legendaries::Dict{Symbol,Vector{Int}}
end

function Collection(sets, ncard_db)::Collection
    Collection(
        generate_collection_table(sets, ncard_db, 1),
        generate_collection_table(sets, ncard_db, 2),
        generate_collection_table(sets, ncard_db, 3),
        generate_collection_table(sets, ncard_db, 4)
    )
end

struct DecklistArray <: RarityGrouping
    commons::Dict{Symbol,Array{Array{Tuple{Int64,Int64},1},1}}
    rares::Dict{Symbol,Array{Array{Tuple{Int64,Int64},1},1}}
    epics::Dict{Symbol,Array{Array{Tuple{Int64,Int64},1},1}}
    legendaries::Dict{Symbol,Array{Array{Tuple{Int64,Int64},1},1}}
end

function Base.:length(x::DecklistArray)
    z = x.commons
    k = first(keys(z))
    length(z[k])
end

struct SingleSetDeckArray <: RarityGrouping
    set::Symbol
    commons::Vector{Tuple{Int,Int}}
    rares::Vector{Tuple{Int,Int}}
    epics::Vector{Tuple{Int,Int}}
    legendaries::Vector{Tuple{Int,Int}}
end

function SingleSetDeckArray(dla::DecklistArray, set::Symbol, deck::Int)
    ans = Vector{Vector{Tuple{Int,Int}}}(undef, 4)
    for (r,x) in enumerate(dla)
        ans[r] = x[set][deck]
    end
    SingleSetDeckArray(set, ans...)
end

function deck_to_int(deck)

    ans = map(i -> Dict(map( x -> (x, Tuple{Int,Int}[]), sets)), 1:4)
    for (name,count) in values(deck)
        set,id,rarity = card_db[name]
        if rarity > 0
            pos = findunique( rarity_db[set][rarity], id)
            push!( ans[rarity][set], (pos, count))
        end
    end
    ans
end

function DecklistArray(decks)::DecklistArray
    ans = map( i -> Dict( map( x-> (x, Vector{Vector{Tuple{Int,Int}}}()), sets)), 1:4)
    for v in values(decks)
        deckarray = deck_to_int(v)
        for i in 1:4
            for set in sets
                push!(ans[i][set], deckarray[i][set])
            end
        end
    end
    DecklistArray(
        ans[1], ans[2], ans[3], ans[4]
    )
end

@with_kw mutable struct AccountState
    dust::Int = 0
    bonus_packs::Dict{Symbol,Int} = Dict(map(x -> (x, 0), sets))
    pitytimer::Dict{Symbol,Vector{Int}} = Dict(map(x -> (x, [0,0,0,30]), sets))
end

@with_kw struct HearthstoneSimParameters
    rarity_probs::Vector{Float64} = [71.4,12.51,2.18,0.53] / 100
    golden_probs::Vector{Float64} = 1 ./ [47.9, 17.1, 18.0, 11.8]
    use_pitytimer::Bool = true
    rarity_by_set::Dict{Symbol,NTuple{4,Int}} = ncard_db
    startfresh::Bool = true
    bonus_dust::Int = 400
    sets_with_bonus::Vector{Symbol} = [:GILNEAS]
    welcome_bundle::Bool = false
end

@with_kw struct SimOutput
    total_packs::SetCounter
    total_packs_for_deck::SetDeckCounter
    cards_from_packs::SetRarityCounter
    pack_card_rarity::Array{Int,2}
end

function SimOutput(sets::Vector{Symbol}, nreps::Int, ndeck::Int)::SimOutput
    SimOutput(
        total_packs = make_counter(sets, nreps),
        total_packs_for_deck = make_counter(sets, (ndeck, nreps)),
        cards_from_packs = make_counter(sets, (nreps, 4)),
        pack_card_rarity = zeros(Int, (nreps, 4))
    )
end

function copy_into_output!(total::SimOutput, partial::SimOutput; startindex::Int, sets = sets)
    endindex = startindex + size(partial.pack_card_rarity,1) - 1
    range = startindex:endindex

    for set in sets
        total.total_packs[set][range] = partial.total_packs[set]
        total.total_packs_for_deck[set][:,range] = partial.total_packs_for_deck[set]
        total.cards_from_packs[set][range,:] = partial.cards_from_packs[set]
    end
    total.pack_card_rarity[range,:] = partial.pack_card_rarity
    endindex + 1
end

@with_kw mutable struct SimState
    rep::Int
    deck::Int
    set::Symbol
    npacks::Int
end

max_stars_repo_name: asieev/HearthstoneSim

max_issues_repo_name: asieev/HearthstoneSim

max_forks_repo_name: asieev/HearthstoneSim


text:
# Copyright (c) 2015: Lyndon White.
# Reference:
# 1. https://github.com/oxinabox/Pipe.jl/blob/master/src/Pipe.jl
# This function is modified from Lyndon White (2015), with the MIT License.

module Pipe
#Reflow piped things replacing the _ in the next section

const PLACEHOLDER = :_

function replace(arg::Any, target)
  arg #Normally do nothing
end

function replace(arg::Symbol, target)
  if arg == PLACEHOLDER
    target
  else
    arg
  end
end

function replace(arg::Expr, target)
  rep = copy(arg)
  rep.args = map(x -> replace(x, target), rep.args)
  rep
end

find_underscores(x::Any) = false
find_underscores(x::Symbol) = x == :_
function find_underscores(x::Expr)
  map(find_underscores, x.args) |> any
end

function rewrite(ff::Expr, target)
  # @show ff.args
  # dump(ff.args)
  # check underscores; if not found, insert to the first position
  _found = find_underscores(ff)

  if !_found
    agrs = ff.args
    if typeof(ff.args[2]) == LineNumberNode
      ff.args = [agrs[1:2]..., :_, agrs[3:end]...]
    else
      ff.args = [agrs[1], :_, agrs[2:end]...]
    end
  end
  # dump(ff.args)

  rep_args = map(x -> replace(x, target), ff.args)
  if ff.args != rep_args
    #_ subsitution
    ff.args = rep_args
    return ff
  end

  #No subsitution was done (no _ found)
  #Apply to a function that is being returned by ff,
  #(ff could be a function call or something more complex)
  rewrite_apply(ff, target)
end

function rewrite_broadcasted(ff::Expr, target)
  temp_var = gensym()
  rep_args = map(x -> replace(x, temp_var), ff.args)
  if ff.args != rep_args
    #_ subsitution
    ff.args = rep_args
    return :($temp_var -> $ff)
  end

  #No subsitution was done (no _ found)
  #Apply to a function that is being returned by ff,
  #(ff could be a function call or something more complex)
  rewrite_apply_broadcasted(ff, target)
end

function rewrite_apply(ff, target)
  :($ff($target)) #function application
end

function rewrite_apply_broadcasted(ff, target)
  temp_var = gensym()
  :($temp_var -> $ff($temp_var))
end

function rewrite(ff::Symbol, target)
  if ff == PLACEHOLDER
    target
  else
    rewrite_apply(ff, target)
  end
end

function rewrite_broadcasted(ff::Symbol, target)
  if ff == PLACEHOLDER
    target
  else
    rewrite_apply_broadcasted(ff, target)
  end
end

function funnel(ee::Any) #Could be a Symbol could be a literal
  ee #first (left most) input
end

function funnel(ee::Expr)
  if (ee.args[1] == :|>)
    target = funnel(ee.args[2]) #Recurse
    rewrite(ee.args[3], target)
  elseif (ee.args[1] == :.|>)
    target = funnel(ee.args[2]) #Recurse
    rewritten = rewrite_broadcasted(ee.args[3], target)
    ee.args[3] = rewritten
    ee
  else
    #Not in a piping situtation
    ee #make no change
  end
end

macro pipe(ee)
  esc(funnel(ee))
end

export @pipe
end


text:
"""
Notes:
 - `rand` will store randomness only when trace type matches TraceR.
 - `randc` never stores randomness.
 - `randr` will store and replay randomness regardless trace type (N.B. Particle Gibbs uses `randr`).
 - `fork1` will perform replaying immediately and fix the particle weight to 1.
 - `fork2` will perform lazy replaying and accumulate likelihoods like a normal particle.
"""

module Traces
using Distributions

include("taskcopy.jl")
include("tarray.jl")

export Trace, TraceR, TraceC, current_trace, fork, fork2, randr, TArray, tzeros

type Trace{T}
  task :: Task
  randomness :: Array{Any, 1}    # elem t is the randomness created by the t’th assume call.
  index :: Int                 # index of current randomness
  num_produce :: Int           # num of produce calls from trace, each produce corresponds to an observe.
  Trace() = (res = new(); res.randomness = Array{Any,1}(); res.index = 0; res.num_produce = 0; res)
end

function call{T}(::Type{Trace{T}}, f::Function)
  res = Trace{T}();
  res.task = Task(()->f());
  if isa(res.task.storage, Void)
    res.task.storage = ObjectIdDict()
  end
  res.task.storage[:turing_trace] = res # create a backward reference in task_local_storage
  res
end

typealias TraceR Trace{:R} # Task Copy
typealias TraceC Trace{:C} # Replay

# generate a new random variable, replay if t.counter < length(t.randomness)
function randr( t::Trace, distr :: Distribution )
  t.index += 1
  if t.index < length(t.randomness)
    res = t.randomness[t.index]
  else # sample, record
    res = Distributions.rand(distr)
    push!(t.randomness, res)
  end
  return res
end

# generate a new random variable, no replay
randc(t::Trace, distr :: Distribution) = Distributions.rand(distr)

Distributions.rand(t::TraceR, distr :: Distribution) = randr(t, distr)
Distributions.rand(t::TraceC, distr :: Distribution) = randc(t, distr)

# step to the next observe statement, return log likelihood
Base.consume(t::Trace) = (t.num_produce += 1; Base.consume(t.task))

# Task copying version of fork for both TraceR and TraceC.
function forkc(trace :: Trace)
  newtrace = typeof(trace)()
  newtrace.task = Base.copy(trace.task)
  n_rand = min(trace.index, length(trace.randomness))
  newtrace.randomness = trace.randomness[1:n_rand]
  newtrace.index = trace.index
  newtrace.num_produce = trace.num_produce
  newtrace.task.storage[:turing_trace] = newtrace
  newtrace
end

# fork s and replay until observation t; drop randomness between y_t:T if keep == false
#  N.B.: PG requires keeping all randomness even we only replay up to observation y_t
function forkr(trace :: TraceR, t :: Int, keep :: Bool)
  # Step 0: create new task and copy randomness
  newtrace = TraceR(trace.task.code)
  newtrace.randomness = deepcopy(trace.randomness)

  # Step 1: Call consume t times to replay randomness
  map((i)->consume(newtrace), 1:t)

  # Step 2: Remove remaining randomness if keep==false
  if !keep
    newtrace.randomness = newtrace.randomness[1:newtrace.index]
  end
  newtrace
end

# Default fork implementation, replay immediately.
fork(s :: TraceR) = forkr(s, s.num_produce, false)
fork(s :: TraceC) = forkc(s)

# Lazily replay on demand, note that:
#  - lazy replay is only possible for TraceR
#  - lazy replay accumulates likelihoods
#  - lazy replay is useful for implementing PG (i.e. ref particle)
fork2(s :: TraceR) = forkr(s, 0, true)

current_trace() = current_task().storage[:turing_trace]

end