struct RootCondition <: SamplingCondition end
struct NonExtinctCondition <: SamplingCondition end
struct NoCondition <: SamplingCondition end

struct NowhereExtinctCondition{T} <: SamplingCondition
    s::Vector{T}  # signs for inclusion - exclusion
    NowhereExtinctCondition(model, T=Float64) =
        new{T}(treepgf_allbinary_sign(model))
end

condition(wm::WhaleModel) = condition(wm, wm.condition)
condition(wm::WhaleModel, ::NoCondition) = 0.

# log probability of non-extinction
function condition(wm::WhaleModel, ::NonExtinctCondition)
    @unpack 畏 = get胃(wm.rates, root(wm))
    log(1. -geompgf(畏, get系(root(wm))))
end

# log probability of non extinction in both clades stemming from the root
function condition(wm::WhaleModel, ::RootCondition)
    @unpack 畏 = get胃(wm.rates, root(wm))
    f, g = children(root(wm))
    系r = geompgf(畏, get系(root(wm)))
    系f = geompgf(畏, get系(f))
    系g = geompgf(畏, get系(g))
    p = one(畏) - 系f - 系g + 系r
    p > zero(p) ? log(p) : -Inf
end

function condition(wm::WhaleModel, e::NowhereExtinctCondition{T}) where T
    p = treepgf_allbinary(wm)
    # probability of at least one count at each leaf
    c = one(T) - sum((p .* e.s)[1:end-1])
    zero(c) < c < one(c) ? log(c) : -Inf
end

# log probability of non-extinction everywhere? naive implementation
# function condition(wm::WhaleModel{T}, ::NowhereCondition, bound=10) where T
#     ���� = zeros(T, bound, length(wm.order))
#     p = one(T)
#     function walk(n)
#         _pvec!(����, wm, n)
#         for c in children(n) walk(c) end
#         if isleaf(n)
#             p *= sum(����[2:end, id(n)])
#         end
#         return
#     end
#     walk(root(wm))
#     return log(p)
# end
#
# function _pvec!(����, model, n)
#     if isroot(n)
#         @unpack 畏 = get胃(model.rates, n)
#         ����[:,id(n)] = [0. ; pdf.(Geometric(畏), 0:size(����)[1]-2)]
#     else
#         @unpack 位, 渭 = get胃(model.rates, n)
#         t = distance(n)
#         bound = size(����)[1]
#         matrix = [tp(i, j, t, 位, 渭) for i=0:bound-1, j=0:bound-1]
#         p = matrix' * ����[:,id(parent(n))]
#         ����[:,id(n)] .= p /sum(p)
#     end
# end
