text:
"""
$(SIGNATURES)

Unified gas kinetic scheme (UGKS)

* @args: particle distribution functions and their slopes at left/right sides of interface
* @args: particle velocity quadrature points and weights
* @args: time step
"""
function flux_ugks!(
    fw::T1,
    fh::T2,
    fb::T2,
    wL::T3,
    hL::T4,
    bL::T4,
    wR::T3,
    hR::T4,
    bR::T4,
    u::T5,
    ω::T5,
    inK,
    γ,
    visRef,
    visIdx,
    pr,
    dt,
    dxL,
    dxR,
    shL = zeros(eltype(hL), axes(hL))::T4,
    sbL = zeros(eltype(bL), axes(bL))::T4,
    shR = zeros(eltype(hR), axes(hR))::T4,
    sbR = zeros(eltype(bR), axes(bR))::T4,
) where {T1<:AA{<:FN,1},T2<:AA{<:FN,1},T3<:AA{<:Real,1},T4<:AA{<:FN,1},T5<:AA{<:FN,1}} # 1D2F flux

    #--- reconstruct initial distribution ---#
    δ = heaviside.(u)

    h = @. hL * δ + hR * (1.0 - δ)
    b = @. bL * δ + bR * (1.0 - δ)
    sh = @. shL * δ + shR * (1.0 - δ)
    sb = @. sbL * δ + sbR * (1.0 - δ)

    #--- construct interface variables ---#
    #w = moments_conserve(h, b, u, v, ω)
    #prim = conserve_prim(w, γ)

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)
    Mu1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Muv1 = moments_conserve(MuL1, Mxi1, 0, 0)
    Mu2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)
    Muv2 = moments_conserve(MuR2, Mxi2, 0, 0)

    w = @. primL[1] * Muv1 + primR[1] * Muv2
    prim = conserve_prim(w, γ)
    aL = pdf_slope(prim, (w .- wL) ./ dxL, inK)
    aR = pdf_slope(prim, (wR .- w) ./ dxR, inK)

    Mu, Mxi, MuL, MuR = gauss_moments(prim, inK)
    MauL = moments_conserve_slope(aL, MuL, Mxi, 1)
    MauR = moments_conserve_slope(aR, MuR, Mxi, 1)
    aT = pdf_slope(prim, -prim[1] .* (MauL .+ MauR), inK)

    #--- calculate integral time constants ---#
    tau = vhs_collision_time(prim, visRef, visIdx)

    Mt = zeros(5)
    Mt[4] = tau * (1.0 - exp(-dt / tau)) # f0
    Mt[5] = -tau * dt * exp(-dt / tau) + tau * Mt[4]
    Mt[1] = dt - Mt[4] # M0
    Mt[2] = -tau * Mt[1] + Mt[5]
    Mt[3] = dt^2 / 2.0 - tau * Mt[1]

    #--- calculate flux from M0 ---#
    Muv = moments_conserve(Mu, Mxi, 1, 0)
    MauL = moments_conserve_slope(aL, MuL, Mxi, 2)
    MauR = moments_conserve_slope(aR, MuR, Mxi, 2)
    MauT = moments_conserve_slope(aT, Mu, Mxi, 1)

    @. fw = Mt[1] * prim[1] * Muv + Mt[2] * prim[1] * (MauL + MauR) + Mt[3] * prim[1] * MauT

    #--- calculate flux from f0 ---#
    H = maxwellian(u, prim)
    B = H .* inK ./ (2.0 * prim[end])

    fw[1] += Mt[4] * sum(ω .* u .* h) - Mt[5] * sum(ω .* u .^ 2 .* sh)
    fw[2] += Mt[4] * sum(ω .* u .^ 2 .* h) - Mt[5] * sum(ω .* u .^ 3 .* sh)
    fw[3] +=
        Mt[4] * 0.5 * (sum(ω .* u .^ 3 .* h) + sum(ω .* u .* b)) -
        Mt[5] * 0.5 * (sum(ω .* u .^ 4 .* sh) + sum(ω .* u .^ 2 .* sb))

    @. fh =
        Mt[1] * u * H +
        Mt[2] * u^2 * (aL[1] * H + aL[2] * u * H + 0.5 * aL[3] * (u^2 * H + B)) * δ +
        Mt[2] *
        u^2 *
        (aR[1] * H + aR[2] * u * H + 0.5 * aR[3] * (u^2 * H + B)) *
        (1.0 - δ) +
        Mt[3] * u * (aT[1] * H + aT[2] * u * H + 0.5 * aT[3] * (u^2 * H + B)) +
        Mt[4] * u * h - Mt[5] * u^2 * sh
    @. fb =
        Mt[1] * u * B +
        Mt[2] *
        u^2 *
        (aL[1] * B + aL[2] * u * B + 0.5 * aL[3] * (u^2 * B + Mxi[2] * H)) *
        δ +
        Mt[2] *
        u^2 *
        (aR[1] * B + aR[2] * u * B + 0.5 * aR[3] * (u^2 * B + Mxi[2] * H)) *
        (1.0 - δ) +
        Mt[3] * u * (aT[1] * B + aT[2] * u * B + 0.5 * aT[3] * (u^2 * B + Mxi[2] * H)) +
        Mt[4] * u * b - Mt[5] * u^2 * sb

    return nothing

end

"""
$(SIGNATURES)
"""
function flux_ugks!(
    fw::T1,
    fh::T2,
    fb::T2,
    wL::T3,
    hL::T4,
    bL::T4,
    wR::T3,
    hR::T4,
    bR::T4,
    u::T5,
    v::T5,
    ω::T5,
    inK,
    γ,
    visRef,
    visIdx,
    pr,
    dt,
    dxL,
    dxR,
    len,
    shL = zeros(eltype(hL), axes(hL))::T4,
    sbL = zeros(eltype(bL), axes(bL))::T4,
    shR = zeros(eltype(hR), axes(hR))::T4,
    sbR = zeros(eltype(bR), axes(bR))::T4,
) where {T1<:AA{<:FN,1},T2<:AA{<:FN,2},T3<:AA{<:Real,1},T4<:AA{<:FN,2},T5<:AA{<:FN,2}} # 2D2F flux

    #--- reconstruct initial distribution ---#
    δ = heaviside.(u)

    h = @. hL * δ + hR * (1.0 - δ)
    b = @. bL * δ + bR * (1.0 - δ)
    sh = @. shL * δ + shR * (1.0 - δ)
    sb = @. sbL * δ + sbR * (1.0 - δ)

    #--- construct interface variables ---#
    #w = moments_conserve(h, b, u, v, ω)
    #prim = conserve_prim(w, γ)

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)
    Mu1, Mv1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Muv1 = moments_conserve(MuL1, Mv1, Mxi1, 0, 0, 0)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)
    Muv2 = moments_conserve(MuR2, Mv2, Mxi2, 0, 0, 0)
    w = @. primL[1] * Muv1 + primR[1] * Muv2
    prim = conserve_prim(w, γ)

    aL = pdf_slope(prim, (w .- wL) ./ dxL, inK)
    aR = pdf_slope(prim, (wR .- w) ./ dxR, inK)

    Mu, Mv, Mxi, MuL, MuR = gauss_moments(prim, inK)
    MauL = moments_conserve_slope(aL, MuL, Mv, Mxi, 1, 0)
    MauR = moments_conserve_slope(aR, MuR, Mv, Mxi, 1, 0)
    aT = pdf_slope(prim, -prim[1] .* (MauL .+ MauR), inK)

    #--- calculate integral time constants ---#
    tau = vhs_collision_time(prim, visRef, visIdx)

    Mt = zeros(5)
    Mt[4] = tau * (1.0 - exp(-dt / tau)) # f0
    Mt[5] = -tau * dt * exp(-dt / tau) + tau * Mt[4]
    Mt[1] = dt - Mt[4] # M0
    Mt[2] = -tau * Mt[1] + Mt[5]
    Mt[3] = dt^2 / 2.0 - tau * Mt[1]

    # --- calculate flux from M0 ---#
    Muv = moments_conserve(Mu, Mv, Mxi, 1, 0, 0)
    MauL = moments_conserve_slope(aL, MuL, Mv, Mxi, 2, 0)
    MauR = moments_conserve_slope(aR, MuR, Mv, Mxi, 2, 0)
    MauT = moments_conserve_slope(aT, Mu, Mv, Mxi, 1, 0)

    @. fw = Mt[1] * prim[1] * Muv + Mt[2] * prim[1] * (MauL + MauR) + Mt[3] * prim[1] * MauT

    # --- calculate flux from f0 ---#
    H = maxwellian(u, v, prim)
    B = H .* inK ./ (2.0 * prim[end])

    fw[1] += Mt[4] * sum(ω .* u .* h) - Mt[5] * sum(ω .* u .^ 2 .* sh)
    fw[2] += Mt[4] * sum(ω .* u .^ 2 .* h) - Mt[5] * sum(ω .* u .^ 3 .* sh)
    fw[3] += Mt[4] * sum(ω .* v .* u .* h) - Mt[5] * sum(ω .* v .* u .^ 2 .* sh)
    fw[4] +=
        Mt[4] * 0.5 * (sum(ω .* u .* (u .^ 2 .+ v .^ 2) .* h) + sum(ω .* u .* b)) -
        Mt[5] *
        0.5 *
        (sum(ω .* u .^ 2 .* (u .^ 2 .+ v .^ 2) .* sh) + sum(ω .* u .^ 2 .* sb))

    @. fh =
        Mt[1] * u * H +
        Mt[2] *
        u^2 *
        (aL[1] * H + aL[2] * u * H + aL[3] * v * H + 0.5 * aL[4] * ((u^2 + v^2) * H + B)) *
        δ +
        Mt[2] *
        u^2 *
        (aR[1] * H + aR[2] * u * H + aR[3] * v * H + 0.5 * aR[4] * ((u^2 + v^2) * H + B)) *
        (1.0 - δ) +
        Mt[3] *
        u *
        (aT[1] * H + aT[2] * u * H + aT[3] * v * H + 0.5 * aT[4] * ((u^2 + v^2) * H + B)) +
        Mt[4] * u * h - Mt[5] * u^2 * sh
    @. fb =
        Mt[1] * u * B +
        Mt[2] *
        u^2 *
        (
            aL[1] * B +
            aL[2] * u * B +
            aL[3] * v * B +
            0.5 * aL[4] * ((u^2 + v^2) * B + Mxi[2] * H)
        ) *
        δ +
        Mt[2] *
        u^2 *
        (
            aR[1] * B +
            aR[2] * u * B +
            aR[3] * v * B +
            0.5 * aR[4] * ((u^2 + v^2) * B + Mxi[2] * H)
        ) *
        (1.0 - δ) +
        Mt[3] *
        u *
        (
            aT[1] * B +
            aT[2] * u * B +
            aT[3] * v * B +
            0.5 * aT[4] * ((u^2 + v^2) * B + Mxi[2] * H)
        ) +
        Mt[4] * u * b - Mt[5] * u^2 * sb

    # multiply interface length
    fw .*= len
    fh .*= len
    fb .*= len

    return nothing

end

"""
$(SIGNATURES)

3F2V with AAP model
"""
function flux_ugks!(
    fw::T1,
    fh0::T2,
    fh1::T2,
    fh2::T2,
    wL::T3,
    h0L::T4,
    h1L::T4,
    h2L::T4,
    wR::T3,
    h0R::T4,
    h1R::T4,
    h2R::T4,
    u::T5,
    v::T5,
    ω::T5,
    inK,
    γ,
    mi,
    ni,
    me,
    ne,
    Kn,
    dt,
    dxL,
    dxR,
    len,
    sh0L = zeros(eltype(h0L), axes(h0L))::T4,
    sh1L = zeros(eltype(h1L), axes(h1L))::T4,
    sh2L = zeros(eltype(h2L), axes(h2L))::T4,
    sh0R = zeros(eltype(h0R), axes(h0R))::T4,
    sh1R = zeros(eltype(h1R), axes(h1R))::T4,
    sh2R = zeros(eltype(h2R), axes(h2R))::T4,
) where {T1<:AA{<:FN,2},T2<:AA{<:FN,3},T3<:AA{<:Real,2},T4<:AA{<:FN,3},T5<:AA{<:FN,3}}

    #--- reconstruct initial distribution ---#
    δ = heaviside.(u)

    h0 = @. h0L * δ + h0R * (1.0 - δ)
    h1 = @. h1L * δ + h1R * (1.0 - δ)
    h2 = @. h2L * δ + h2R * (1.0 - δ)

    sh0 = @. sh0L * δ + sh0R * (1.0 - δ)
    sh1 = @. sh1L * δ + sh1R * (1.0 - δ)
    sh2 = @. sh2L * δ + sh2R * (1.0 - δ)

    primL = mixture_conserve_prim(wL, γ)
    primR = mixture_conserve_prim(wR, γ)

    #--- construct interface distribution ---#
    Mu1, Mv1, Mxi1, MuL1, MuR1 = mixture_gauss_moments(primL, inK)
    Muv1 = mixture_moments_conserve(MuL1, Mv1, Mxi1, 0, 0, 0)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = mixture_gauss_moments(primR, inK)
    Muv2 = mixture_moments_conserve(MuR2, Mv2, Mxi2, 0, 0, 0)

    w = similar(wL)
    for j in axes(w, 2)
        @. w[:, j] = primL[1, j] * Muv1[:, j] + primR[1, j] * Muv2[:, j]
    end
    prim = mixture_conserve_prim(w, γ)
    tau = aap_hs_collision_time(prim, mi, ni, me, ne, Kn)
    prim = aap_hs_prim(prim, tau, mi, ni, me, ne, Kn)

    a = mixture_pdf_slope(prim, (wR .- wL) ./ (dxR + dxL), inK)
    Mu, Mv, Mxi, MuL, MuR = mixture_gauss_moments(prim, inK)
    Mau = mixture_moments_conserve_slope(a, Mu, Mv, Mxi, 1, 0, 0)
    aT = mixture_pdf_slope(prim, -prim[1] .* Mau, inK)

    Mt = zeros(2, 2)
    @. Mt[2, :] = tau * (1.0 - exp(-dt / tau)) # f0
    @. Mt[1, :] = dt - Mt[2, :] # M0

    Mt = zeros(5, axes(w, 2))
    for j in axes(Mt, 2)
        Mt[4, j] = tau[j] * (1.0 - exp(-dt / tau[j])) # f0
        Mt[5, j] = -tau[j] * dt * exp(-dt / tau[j]) + tau[j] * Mt[4, j]
        Mt[1, j] = dt - Mt[4, j] # M0
        Mt[2, j] = -tau[j] * Mt[1, j] + Mt[5, j]
        Mt[3, j] = dt^2 / 2.0 - tau[j] * Mt[1, j]
    end

    #--- calculate interface flux ---#
    # flux from M0
    Muv = mixture_moments_conserve(Mu, Mv, Mxi, 1, 0, 0)
    Mau = mixture_moments_conserve_slope(a, Mu, Mv, Mxi, 2, 0, 0)
    MauT = mixture_moments_conserve_slope(aT, Mu, Mv, Mxi, 1, 0, 0)
    for j in axes(fw, 2)
        @. fw[:, j] =
            Mt[1, j] * prim[1, j] * Muv[:, j] +
            Mt[2, j] * prim[1, j] * Mau[:, j] +
            Mt[3, j] * prim[1, j] * MauT[:, j]
    end

    # flux from f0
    H0 = mixture_maxwellian(u, v, prim)
    H1 = similar(H0)
    H2 = similar(H0)
    for j in axes(H0, 3)
        H1[:, :, j] = H0[:, :, j] .* prim[4, j]
        H2[:, :, j] .= H0[:, :, j] .* (prim[4, j]^2 + 1.0 / (2.0 * prim[5, j]))
    end

    for j in axes(fw, 2)
        fw[1, j] +=
            Mt[4, j] * sum(ω[:, :, j] .* u[:, :, j] .* h0[:, :, j]) -
            Mt[5, j] * sum(ω[:, :, j] .* u[:, :, j] .^ 2 .* sh0[:, :, j])
        fw[2, j] +=
            Mt[4, j] * sum(ω[:, :, j] .* u[:, :, j] .^ 2 .* h0[:, :, j]) -
            Mt[5, j] * sum(ω[:, :, j] .* u[:, :, j] .^ 3 .* sh0[:, :, j])
        fw[3, j] +=
            Mt[4, j] * sum(ω[:, :, j] .* v[:, :, j] .* u[:, :, j] .* h0[:, :, j]) -
            Mt[5, j] * sum(ω[:, :, j] .* v[:, :, j] .* u[:, :, j] .^ 2 .* sh0[:, :, j])
        fw[4, j] +=
            Mt[4, j] * sum(ω[:, :, j] .* u[:, :, j] .* h1[:, :, j]) -
            Mt[5, j] * sum(ω[:, :, j] .* u[:, :, j] .^ 2 .* sh1[:, :, j])
        fw[5, j] +=
            Mt[4, j] *
            0.5 *
            (
                sum(
                    ω[:, :, j] .* u[:, :, j] .* (u[:, :, j] .^ 2 .+ v[:, :, j] .^ 2) .*
                    h0[:, :, j],
                ) + sum(ω[:, :, j] .* u[:, :, j] .* h2[:, :, j])
            ) -
            Mt[5, j] *
            0.5 *
            (
                sum(
                    ω[:, :, j] .* u[:, :, j] .^ 2 .* (u[:, :, j] .^ 2 .+ v[:, :, j] .^ 2) .*
                    sh0[:, :, j],
                ) + sum(ω[:, :, j] .* u[:, :, j] .^ 2 .* sh2[:, :, j])
            )

        @. fh0[:, :, j] =
            Mt[1, j] * u[:, :, j] * H0[:, :, j] +
            Mt[2, j] *
            u[:, :, j]^2 *
            (
                a[1, j] * H0[:, :, j] +
                a[2, j] * u[:, :, j] * H0[:, :, j] +
                a[3, j] * v[:, :, j] * H0[:, :, j] +
                a[4, j] * u[:, :, j] * H1[:, :, j] +
                0.5 * a[5, j] * ((u[:, :, j]^2 + v[:, :, j]^2) * H0[:, :, j] + H2[:, :, j])
            ) +
            Mt[3, j] *
            u[:, :, j] *
            (
                aT[1, j] * H0[:, :, j] +
                aT[2, j] * u[:, :, j] * H0[:, :, j] +
                aT[3, j] * v[:, :, j] * H0[:, :, j] +
                aT[4, j] * u[:, :, j] * H1[:, :, j] +
                0.5 * aT[5, j] * ((u[:, :, j]^2 + v[:, :, j]^2) * H0[:, :, j] + H2[:, :, j])
            ) +
            Mt[4, j] * u[:, :, j] * h0[:, :, j] - Mt[5, j] * u[:, :, j]^2 * sh0[:, :, j]
        @. fh1[:, :, j] =
            Mt[1, j] * u[:, :, j] * H1[:, :, j] +
            Mt[2, j] *
            u[:, :, j]^2 *
            (
                a[1, j] * H1[:, :, j] +
                a[2, j] * u[:, :, j] * H1[:, :, j] +
                a[3, j] * v[:, :, j] * H1[:, :, j] +
                a[4, j] * u[:, :, j] * H2[:, :, j] +
                0.5 *
                a[5, j] *
                ((u[:, :, j]^2 + v[:, :, j]^2) * H1[:, :, j] + Mxi[3, j] * H0[:, :, j])
            ) +
            Mt[3, j] *
            u[:, :, j] *
            (
                aT[1, j] * H1[:, :, j] +
                aT[2, j] * u[:, :, j] * H1[:, :, j] +
                aT[3, j] * v[:, :, j] * H1[:, :, j] +
                aT[4, j] * u[:, :, j] * H2[:, :, j] +
                0.5 *
                aT[5, j] *
                ((u[:, :, j]^2 + v[:, :, j]^2) * H1[:, :, j] + Mxi[3, j] * H0[:, :, j])
            ) +
            Mt[4, j] * u[:, :, j] * h1[:, :, j] - Mt[5, j] * u[:, :, j]^2 * sh1[:, :, j]
        @. fh2[:, :, j] =
            Mt[1, j] * u[:, :, j] * H2[:, :, j] +
            Mt[2, j] *
            u[:, :, j]^2 *
            (
                a[1, j] * H2[:, :, j] +
                a[2, j] * u[:, :, j] * H2[:, :, j] +
                a[3, j] * v[:, :, j] * H2[:, :, j] +
                a[4, j] * u[:, :, j] * Mxi[3, j] * H0[:, :, j] +
                0.5 *
                a[5, j] *
                ((u[:, :, j]^2 + v[:, :, j]^2) * H2[:, :, j] + Mxi[4, j] * H0[:, :, j])
            ) +
            Mt[3, j] *
            u[:, :, j] *
            (
                aT[1, j] * H2[:, :, j] +
                aT[2, j] * u[:, :, j] * H2[:, :, j] +
                aT[3, j] * v[:, :, j] * H2[:, :, j] +
                aT[4, j] * u[:, :, j] * Mxi[3, j] * H0[:, :, j] +
                0.5 *
                aT[5, j] *
                ((u[:, :, j]^2 + v[:, :, j]^2) * H2[:, :, j] + Mxi[4, j] * H0[:, :, j])
            ) +
            Mt[4, j] * u[:, :, j] * h2[:, :, j] - Mt[5, j] * u[:, :, j]^2 * sh2[:, :, j]
    end

    @. fw *= len
    @. fh0 *= len
    @. fh1 *= len
    @. fh2 *= len

    return nothing

end


text:
@static if VERSION < v"1.1"
    isnothing(::Any) = false
    isnothing(::Nothing) = true
end

@static if VERSION < v"1.2"
    function map!(f, iter::Base.ValueIterator)
        dict = iter.dict
        for (key, val) in pairs(dict)
            dict[key] = f(val)
        end
        return iter
    end
end

"Return all vars in a term."
get_vars(t::Term) = error("Not implemented.")
get_vars(t::Const) = Set{Var}()
get_vars(t::Var) = Set{Var}([t])
get_vars(t::Compound) =
    length(t.args) > 0 ? union((get_vars(a) for a in t.args)...) : Set{Var}()

"Check if a term is ground (contains no variables)."
is_ground(t::Term) = error("Not implemented.")
is_ground(t::Const) = true
is_ground(t::Var) = false
is_ground(t::Compound) = all(is_ground(a) for a in t.args)

"Check whether a variable appears in a term."
occurs_in(v::Var, t::Term) = error("Not implemented.")
occurs_in(v::Var, t::Const) = false
occurs_in(v::Var, t::Var) = (v.name == t.name)
occurs_in(v::Var, t::Compound) = any(occurs_in(v, a) for a in t.args)

"Performs variable substitution of var by val in a term."
substitute(term::Term, var::Var, val::Term) = error("Not implemented.")
substitute(term::Const, var::Var, val::Term) = term
substitute(term::Var, var::Var, val::Term) = term.name == var.name ? val : term
function substitute(term::Compound, var::Var, val::Term)
    args, ident = nothing, true
    for (i, a) in enumerate(term.args)
        b = substitute(a, var, val)
        if ident && a !== b
            args, ident = collect(term.args[1:i-1]), false
        end
        if !ident push!(args, b) end
    end
    return ident ? term : Compound(term.name, args)
end

"Apply substitution to a term."
substitute(term::Term, subst::Subst) = error("Not implemented.")
substitute(term::Const, subst::Subst) = term
substitute(term::Var, subst::Subst) = get(subst, term, term)
function substitute(term::Compound, subst::Subst)
    args, ident = nothing, true
    for (i, a) in enumerate(term.args)
        b = substitute(a, subst)
        if ident && a !== b
            args, ident = collect(term.args[1:i-1]), false
        end
        if !ident push!(args, b) end
    end
    return ident ? term : Compound(term.name, args)
end

"Compose two substitutions (s2 after s1)."
function compose(s1::Subst, s2::Subst)
    subst = Subst(var => substitute(val, s2) for (var, val) in s1)
    return merge(s2, subst)
end

"Compose two substitutions (s2 after s1), modifying s1 in place."
function compose!(s1::Subst, s2::Subst)
    map!(v -> substitute(v, s2), values(s1))
    for (var, val) in s2 get!(s1, var, val) end
    return s1
end

"Replace variables in a term or clause with fresh names."
function freshen(term::Term, vars::Set{Var})
    vmap = Subst(v => Var(gensym(v.name)) for v in vars)
    term = substitute(term, vmap)
    return term, vmap
end

function freshen(clause::Clause, vars::Set{Var})
    vmap = Subst(v => Var(gensym(v.name)) for v in vars)
    clause = Clause(substitute(clause.head, vmap),
                    Term[substitute(t, vmap) for t in clause.body])
    return clause, vmap
end

freshen(t::Term) = freshen!(t, Subst())
freshen(c::Clause) = freshen!(c, Subst())

freshen!(t::Const, vmap::Subst) = t
freshen!(t::Var, vmap::Subst) = get!(vmap, t, Var(gensym(t.name)))
freshen!(t::Compound, vmap::Subst) =
    Compound(t.name, Term[freshen!(a, vmap) for a in t.args])
freshen!(c::Clause, vmap::Subst) =
    Clause(freshen!(c.head, vmap), Term[freshen!(t, vmap) for t in c.body])

freshen!(t::Const, vmap::Subst, vcount::Ref{UInt}) = t
freshen!(t::Var, vmap::Subst, vcount::Ref{UInt}) =
    get!(vmap, t, Var(vcount[] += 1))
freshen!(t::Compound, vmap::Subst, vcount::Ref{UInt}) =
    Compound(t.name, Term[freshen!(a, vmap, vcount) for a in t.args])
freshen!(c::Clause, vmap::Subst, vcount::Ref{UInt}) =
    Clause(freshen!(c.head, vmap, vcount),
           Term[freshen!(t, vmap, vcount) for t in c.body])

"Check whether a term has a matching subterm."
function has_subterm(term::Term, subterm::Term)
    if !isnothing(unify(term, subterm)) return true end
    return any(has_subterm(arg, subterm) for arg in get_args(term))
end

"Find all matching subterms in a term."
function find_subterms(term::Term, subterm::Term)
    init = !isnothing(unify(term, subterm)) ? Term[term] : Term[]
    subterms = (find_subterms(a, subterm) for a in get_args(term))
    return reduce(vcat, subterms; init=init)
end

"Convert a vector of Julia objects to a Julog list of constants."
to_const_list(v::Vector) =
    foldr((i, j) -> Compound(:cons, [Const(i), j]), v; init=Compound(:cend, []))

"Convert a vector of Julog terms to a Julog list."
to_term_list(v::Vector{<:Term}) =
    foldr((i, j) -> Compound(:cons, [i, j]), v; init=Compound(:cend, []))

"Convert a list of Julog terms to a vector of Julog terms.."
to_julia_list(list::Term) =
    list.name == :cons ? [list.args[1]; to_julia_list(list.args[2])] : []

"Rewrite implications using and, or and not, subsume true and false."
function simplify(term::Compound)
    if !(term.name in logicals) return term end
    args = simplify.(term.args)
    if term.name in [:imply, :(=>)]
        cond, body = args
        return simplify(@julog or(not(:cond), and(:cond, :body)))
    elseif term.name == :and
        if length(args) == 1 return args[1]
        elseif any(a -> a.name == false, args) return Const(false)
        elseif all(a -> a.name == true, args) return Const(true)
        else term = Compound(term.name, filter!(a -> a.name != true, args)) end
    elseif term.name == :or
        if length(args) == 1 return args[1]
        elseif any(a -> a.name == true, args) return Const(true)
        elseif all(a -> a.name == false, args) return Const(false)
        else term = Compound(term.name, filter!(a -> a.name != false, args)) end
    elseif term.name in [:not, :!]
        if args[1].name == true return Const(false)
        elseif args[1].name == false return Const(true)
        else return Compound(:not, args) end
    end
    return length(term.args) == 1 ? term.args[1] : term
end
simplify(term::Const) = term
simplify(term::Var) = term

"Convert a term to negation normal form."
function to_nnf(term::Compound)
    term = simplify(term)
    if !(term.name in logicals) return term end
    if term.name in [:not, :!]
        inner = term.args[1]
        if inner.name in [:not, :!]
            term = inner.args[1]
        elseif inner.name in [true, false]
            term = Const(!inner.name)
        elseif inner.name in [:and, :or]
            args = to_nnf.(@julog(not(:a)) for a in inner.args)
            term = Compound(inner.name == :and ? :or : :and, args)
        end
    else
        term = Compound(term.name, to_nnf.(term.args))
    end
    return term
end
to_nnf(term::Const) = term
to_nnf(term::Var) = term

"Convert a term to conjunctive normal form."
function to_cnf(term::Compound)
    term = to_nnf(term)
    if !(term.name in [:and, :or]) return @julog and(or(:term)) end
    subterms = to_cnf.(term.args)
    if term.name == :and
        args = foldl(vcat, [a.args for a in subterms]; init=Compound[])
        term = Compound(:and, args)
    elseif term.name == :or
        stack = Compound[@julog(or())]
        for subterm in subterms
            new_stack = Compound[]
            for disj_i in stack
                for disj_j in subterm.args
                    new_disj = Compound(:or, [disj_i.args; disj_j.args])
                    push!(new_stack, new_disj)
                end
            end
            stack = new_stack
        end
        term = Compound(:and, stack)
    end
    return term
end
to_cnf(term::Const) = @julog and(or(:term))
to_cnf(term::Var) = @julog and(or(:term))

"Convert a term to disjunctive normal form."
function to_dnf(term::Compound)
    term = to_nnf(term)
    if !(term.name in [:and, :or]) return @julog or(and(:term)) end
    subterms = to_dnf.(term.args)
    if term.name == :or
        args = foldl(vcat, (a.args for a in subterms); init=Compound[])
        term = Compound(:or, args)
    elseif term.name == :and
        stack = Compound[@julog(and())]
        for subterm in subterms
            new_stack = Compound[]
            for conj_i in stack
                for conj_j in subterm.args
                    new_conj = Compound(:and, [conj_i.args; conj_j.args])
                    push!(new_stack, new_conj)
                end
            end
            stack = new_stack
        end
        term = Compound(:or, stack)
    end
    return term
end
to_dnf(term::Const) = @julog or(and(:term))
to_dnf(term::Var) = @julog or(and(:term))

"Recursively flatten conjunctions in a term to a list."
flatten_conjs(t::Term) = t.name == :and ? flatten_conjs(get_args(t)) : Term[t]
flatten_conjs(t::Vector{<:Term}) = reduce(vcat, flatten_conjs.(t); init=Term[])

"Recursively flatten disjunctions in term to a list."
flatten_disjs(t::Term) = t.name == :or ? flatten_disjs(get_args(t)) : Term[t]
flatten_disjs(t::Vector{<:Term}) = reduce(vcat, flatten_disjs.(t); init=Term[])

"Instantiate universal quantifiers relative to a set of clauses."
function deuniversalize(term::Compound, clauses::Vector{Clause})
    if term.name == :forall
        cond, body = term.args
        sat, subst = resolve(cond, clauses)
        if !sat return Const(true) end
        instantiated = unique!(Term[substitute(body, s) for s in subst])
        return Compound(:and, instantiated)
    elseif term.name in logicals
        args = Term[deuniversalize(a, clauses) for a in term.args]
        return Compound(term.name, args)
    else
        return term
    end
end
deuniversalize(term::Const, ::Vector{Clause}) = term
deuniversalize(term::Var, ::Vector{Clause}) = term
deuniversalize(c::Clause, clauses::Vector{Clause}) =
    Clause(c.head, [deuniversalize(t, clauses) for t in c.body])

"Convert any clauses with disjunctions in their bodies into a set of clauses."
function regularize_clauses(clauses::Vector{Clause}; instantiate::Bool=false)
    regularized = Clause[]
    for c in clauses
        if length(c.body) == 0
            push!(regularized, c)
        else
            body = Compound(:and, c.body)
            if (instantiate == true) body = deuniversalize(body, clauses) end
            for conj in to_dnf(body).args
                push!(regularized, Clause(c.head, conj.args))
            end
        end
    end
    return regularized
end

ext: jl


text:
module Newtonsmethod

using ForwardDiff, LinearAlgebra

function newtonroot(f, f′; x₀, tolerance = 1E-7, maxiter = 1000)
    # setup the algorithm
    x_old = x₀
    normdiff = Inf
    iter = 1
    while normdiff > tolerance && iter <= maxiter
        x_new = x_old - f(x_old)/f′(x_old)
        normdiff = norm(x_new - x_old)
        x_old = x_new
        iter = iter + 1
    end
    if iter < maxiter
        return (value = x_old, normdiff = normdiff, iter = iter)
    else
        return nothing
    end
end

D(f) = x -> ForwardDiff.derivative(f, x)

newtonroot(f; x₀, tolerance = 1E-7, maxiter = 1000) = newtonroot(f, D(f); x₀ = x₀, tolerance=tolerance, maxiter=maxiter)

export newtonroot

end # module


text:
"""
    MLJModelInterface.fit(model, verbosity, data...) -> fitresult, cache, report

All models must implement a `fit` method. Here `data` is the
output of `reformat` on user-provided data, or some some resampling
thereof. The fallback of `reformat` returns the user-provided data
(eg, a table).

"""
function fit end

# fallback for static transformations
fit(::Static, ::Integer, data...) = (nothing, nothing, nothing)

# fallbacks for supervised models that don't support sample weights:
fit(m::Supervised, verbosity, X, y, w) = fit(m, verbosity, X, y)

# fallback for unsupervised annotators when labels or weights appear:
# this is useful for evaluation and mixed composite models that combine
# both supervised and unsupervised annotators
fit(m::UnsupervisedAnnotator, verbosity, X, y) = fit(m, verbosity, X)
fit(m::UnsupervisedAnnotator, verbosity, X, y, w) = fit(m, verbosity, X)

"""
    MLJModelInterface.update(model, verbosity, fitresult, cache, data...)

Models may optionally implement an `update` method. The fallback calls
`fit`.

"""
function update(m::Model, verbosity, fitresult, cache, data...)
    return fit(m, verbosity, data...)
end

"""
    MLJModelInterface.training_losses(model::M, report)

If `M` is an iterative model type which calculates training losses,
implement this method to return an `AbstractVector` of the losses
in historical order. If the model calculates scores instead, then the
sign of the scores should be reversed.

The following trait overload is alse required:
`supports_training_losses(::Type{<:M}) = true`

"""
training_losses(model, report) = nothing

# to support online learning in the future:
# https://github.com/alan-turing-institute/MLJ.jl/issues/60 :
function update_data end

"""
    MLJModelInterface.reformat(model, args...) -> data

Models optionally overload `reformat` to define transformations of
user-supplied data into some model-specific representation (e.g., from
a table to a matrix). When implemented, the MLJ user can avoid
repeating such transformations unnecessarily, and can additionally
make use of more efficient row subsampling, which is then based on the
model-specific representation of data, rather than the
user-representation. When `reformat` is overloaded,
`selectrows(::Model, ...)` must be as well (see
[`selectrows`](@ref)). Furthermore, the model `fit` method(s), and
operations, such as `predict` and `transform`, must be refactored to
act on the model-specific representions of the data.

To implement the `reformat` data front-end for a model, refer to
"Implementing a data front-end" in the [MLJ
manual](https://alan-turing-institute.github.io/MLJ.jl/dev/adding_models_for_general_use/).

"""
reformat(model::Model, args...) = args

"""
    MLJModelInterface.selectrows(::Model, I, data...) -> sampled_data

A model overloads `selectrows` whenever it buys into the optional
`reformat` front-end for data preprocessing. See [`reformat`](@ref)
for details. The fallback assumes `data` is a tuple and calls
`selectrows(X, I)` for each `X` in `data`, returning the results in a
new tuple of the same length. This call makes sense when `X` is a
table, abstract vector or abstract matrix. In the last two cases, a
new object and *not* a view is returned.

"""
selectrows(::Model, I, data...) = map(X -> selectrows(X, I), data)

# this operation can be optionally overloaded to provide access to
# fitted parameters (eg, coeficients of linear model):
"""
    fitted_params(model, fitresult) -> human_readable_fitresult # named_tuple

Models may overload `fitted_params`. The fallback returns
`(fitresult=fitresult,)`.

Other training-related outcomes should be returned in the `report`
part of the tuple returned by `fit`.

"""
fitted_params(::Model, fitresult) = (fitresult=fitresult,)

"""

    predict(model, fitresult, new_data...)

`Supervised` and `SupervisedAnnotator` models must implement the
`predict` operation. Here `new_data` is the output of `reformat`
called on user-specified data.

"""
function predict end

"""

Models types `M` for which `prediction_type(M) == :probablisitic` may
overload `predict_mean`.

"""
function predict_mean end

"""

Models types `M` for which `prediction_type(M) == :probablisitic` may
overload `predict_mode`.

"""
function predict_mode end

"""

Models types `M` for which `prediction_type(M) == :probablisitic` may
overload `predict_median`.

"""
function predict_median end

"""
`JointProbabilistic` supervised models MUST overload `predict_joint`.

`Probabilistic` supervised models MAY overload `predict_joint`.
"""
function predict_joint end

"""
`Unsupervised` models must implement the `transform` operation.
"""
function transform end

"""

`Unsupervised` models may implement the `inverse_transform` operation.

"""
function inverse_transform end

# models can optionally overload these for enable serialization in a
# custom format:
function save end
function restore end

"""
some meta-models may choose to implement the `evaluate` operations
"""
function evaluate end


text:
# Gurobi callbacks

type CallbackData
    cbdata::Ptr{Void}
    model::Model
end

function gurobi_callback_wrapper(ptr_model::Ptr{Void}, cbdata::Ptr{Void}, where::Cint, userdata::Ptr{Void})

    callback,model = unsafe_pointer_to_objref(userdata)::(Function,Model)
    callback(CallbackData(cbdata,model), where)
    return convert(Cint,0)
end

# User callback function should be of the form:
# callback(cbdata::CallbackData, where::Cint)

function set_callback_func!(model::Model, callback::Function)

    grbcallback = cfunction(gurobi_callback_wrapper, Cint, (Ptr{Void}, Ptr{Void}, Cint, Ptr{Void}))
    usrdata = (callback,model)
    ret = @grb_ccall(setcallbackfunc, Cint, (Ptr{Void}, Ptr{Void}, Any), model.ptr_model, grbcallback, usrdata)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    # we need to keep a reference to the callback function
    # so that it isn't garbage collected
    model.callback = usrdata
    nothing
end

export CallbackData, set_callback_func!

for f in (:cbcut, :cblazy)
    @eval function ($f)(cbdata::CallbackData, ind::Vector{Cint}, val::Vector{Float64},
        sense::Char, rhs::Float64)

        len = length(ind)
        @assert length(val) == len

        ret = @grb_ccall($f, Cint, (Ptr{Void},Cint,Ptr{Cint},Ptr{Float64},
            Char,Float64), cbdata.cbdata, len, ind.-1, val, sense, rhs)
        if ret != 0
            throw(GurobiError(cbdata.model.env, ret))
        end
    end
end

export cbcut, cblazy

function cbsolution(cbdata::CallbackData, sol::Vector{Float64})
    nvar = num_vars(cbdata.model)
    @assert length(sol) >= nvar

    ret = @grb_ccall(cbsolution, Cint, (Ptr{Void},Ptr{Float64}),
        cbdata.cbdata, sol)
    if ret != 0
        throw(GurobiError(cbdata.model.env, ret))
    end
end

function cbget{T}(::Type{T},cbdata::CallbackData, where::Cint, what::Integer)

    out = Array(T,1)
    ret = @grb_ccall(cbget, Cint, (Ptr{Void}, Cint, Cint, Ptr{T}),
        cbdata.cbdata, where, convert(Cint,what), out)
    if ret != 0
        throw(GurobiError(cbdata.model.env, ret))
    end
    return out[1]
end

# Callback constants
# grep GRB_CB gurobi_c.h | awk '{ print "const " substr($2,5) " = " $3; }'
const CB_POLLING = 0
const CB_PRESOLVE = 1
const CB_SIMPLEX = 2
const CB_MIP = 3
const CB_MIPSOL = 4
const CB_MIPNODE = 5
const CB_MESSAGE = 6
const CB_BARRIER = 7

export CB_POLLING, CB_PRESOLVE, CB_SIMPLEX, CB_MIP,
       CB_MIPSOL, CB_MIPNODE, CB_MESSAGE, CB_BARRIER

# grep GRB_CB gurobi_c.h | awk '{ print "(\"" tolower(substr($2,8)) "\"," $3 ")"; }'
const cbconstants = [
("pre_coldel",1000,Cint),
("pre_rowdel",1001,Cint),
("pre_senchg",1002,Cint),
("pre_bndchg",1003,Cint),
("pre_coechg",1004,Cint),
("spx_itrcnt",2000,Float64),
("spx_objval",2001,Float64),
("spx_priminf",2002,Float64),
("spx_dualinf",2003,Float64),
("spx_ispert",2004,Float64),
("mip_objbst",3000,Float64),
("mip_objbnd",3001,Float64),
("mip_nodcnt",3002,Float64),
("mip_solcnt",3003,Cint),
("mip_cutcnt",3004,Cint),
("mip_nodlft",3005,Float64),
("mip_itrcnt",3006,Float64),
###("mipsol_sol",4001),
("mipsol_obj",4002,Float64),
("mipsol_objbst",4003,Float64),
("mipsol_objbnd",4004,Float64),
("mipsol_nodcnt",4005,Float64),
("mipsol_solcnt",4006,Cint),
("mipnode_status",5001,Cint),
###("mipnode_rel",5002),
("mipnode_objbst",5003,Float64),
("mipnode_objbnd",5004,Float64),
("mipnode_nodcnt",5005,Float64),
("mipnode_solcnt",5006,Cint),
##("mipnode_brvar",5007), -- undocumented
##("msg_string",6001), -- not yet implemented:
### documentation is unclear on output type
("runtime",6002, Float64),
("barrier_itrcnt",7001,Cint),
("barrier_primobj",7002,Float64),
("barrier_dualobj",7003,Float64),
("barrier_priminf",7004,Float64),
("barrier_dualinf",7005,Float64),
("barrier_compl",7006,Float64)]

for (cname,what,T) in cbconstants
    fname = Symbol("cbget_$cname")
    @eval ($fname)(cbdata::CallbackData, where::Cint) = cbget($T, cbdata, where, $what)
    eval(Expr(:export,fname))
end

for (fname, what) in ((:cbget_mipsol_sol, 4001), (:cbget_mipnode_rel, 5002))
    @eval function ($fname)(cbdata::CallbackData, where::Cint, out::Vector{Float64})
        nvar = num_vars(cbdata.model)
        @assert length(out) >= nvar
        ret = @grb_ccall(cbget, Cint, (Ptr{Void}, Cint, Cint, Ptr{Float64}),
                         cbdata.cbdata, where, $what, out)
        if ret != 0
            throw(GurobiError(cbdata.model.env, ret))
        end
    end
    @eval function ($fname)(cbdata::CallbackData, where::Cint)
        nvar = num_vars(cbdata.model)
        out = Array(Float64,nvar)
        ($fname)(cbdata, where, out)
        return out
    end
    eval(Expr(:export,fname))
end


text:
##################################
## map for columnwise iterators ##
##################################

import Base.map
function map(f::Function, tdci::TdColumnIterator)
    ## function works on da values only (second part of tuple)!!
    ## function returns array or da
    df = DataFrame()
    for (nam, col) in tdci
        df[nam] = f(col)
    end
    td = Timedata(df, idx(tdci.td))
    return td
end

## getColData
##-----------

## transform output data into a format that can be assigned to
## DataFrame column
function getColData(x::Vector)
    return x
end
function getColData(x::DataArray)
    return x
end
function getColData(x::DataFrame)
    return x[:, 1]
end
function getColData(x::AbstractTimedata)
    return x.vals[:, 1]
end

function map(f::Function, tdvi::TdVariableIterator)
    ## function must not work on DataArray, but on TimeData
    ## function returns array, da or TimeData
    df = DataFrame()
    for tdCol in tdvi
        tdRes = f(tdCol)
        df[names(tdCol)[1]] = getColData(tdRes)
    end
    td = Timedata(df, idx(tdvi.td))
    return td
end

###############################
## map for rowwise iterators ##
###############################

## nthRowElem
##-----------

## extract nth value from given output and transform it into a format
## that does not contain metadata anymore

## the following output formats are allowed for f:
function nthRowElem(x::Vector, n::Int)
    return x[n]
end
function nthRowElem(x::Matrix, n::Int)
    return x[1, n]
end
function nthRowElem(x::DataFrame, n::Int)
    return x[1, n]
end
function nthRowElem(x::AbstractTimedata, n::Int)
    return get(x, 1, n)
end

function map(f::Function, tdri::TdRowIterator)
    nObs, nVars = size(tdri.td)
    vals = Array(Any, nObs, nVars)
    rowCounter = 1;
    for row in tdri # DataFrame
        rowRes = f(row)
        for jj=1:nVars
            vals[rowCounter, jj] = nthRowElem(rowRes, jj)
        end
        rowCounter += 1
    end
    df = anyToDf(vals, names(tdri.td))
    td = Timedata(df, idx(tdri.td))
    return td
end

function map(f::Function, tddi::TdDateIterator)
    nObs, nVars = size(tddi.td)
    vals = Array(Any, nObs, nVars)
    rowCounter = 1;
    for tmRow in tddi # TimeData object
        rowRes = f(tmRow)
        for jj=1:nVars
            vals[rowCounter, jj] = nthRowElem(rowRes, jj)
        end
        rowCounter += 1
    end
    df = anyToDf(vals, names(tddi.td))
    td = Timedata(df, idx(tddi.td))
    return td
end

###################################
## map for elementwise iterators ##
###################################

function map(f::Function, tdei::TdEntryIterator)
    nObs, nVars = size(tdei.td)
    vals = Array(Any, nObs, nVars)
    indRow = 1
    indCol = 1
    for entry in tdei
        vals[indRow, indCol] = f(entry)
        if indRow == nObs # move to next column
            indRow = 1
            indCol += 1
        else
            indRow += 1
        end
    end
    df = anyToDf(vals, names(tdei.td))
    td = Timedata(df, idx(tdei.td))
end

function map(f::Function, tdoi::TdObsIterator)
    nObs, nVars = size(tdoi.td)
    vals = Array(Any, nObs, nVars)
    indRow = 1
    indCol = 1
    for obs in tdoi
        vals[indRow, indCol] = get(f(obs), :)
        if indRow == nObs # move to next column
            indRow = 1
            indCol += 1
        else
            indRow += 1
        end
    end
    df = anyToDf(vals, names(tdoi.td))
    td = Timedata(df, idx(tdoi.td))
end

## ####################################
## ## map with two collections - old ##
## ####################################

## ## map function for column iterator always returns Timedata object of
## ## same size! Indices remain unchanged!

## ## map for two collections:
## function map(f::Function, tdci::TdColumnIterator, x)
##     # note: `f` must return a consistent length
##     nIter1 = length(tdci)
##     nIter2 = length(x)
##     if nIter1 != nIter2
##         error("length of iterators must be identical")
##     end

##     res = DataFrame()
##     nams = names(tdci.td)
##     state = start(x)
##     for ii=1:nIter1
##         val, state = next(x, state)
##         res[nams[ii]] = f(tdci.td.vals[:, ii], val)
##     end

##     display(res)
##     td = Timedata(res, idx(tdci.td))
##     if isa(tdci.td, AbstractTimenum)
##         try
##             td = convert(Timenum, td)
##         catch
##         end
##     elseif isa(tdci.td, AbstractTimenum)
##         try
##             td = convert(Timematr, td)
##         catch
##         end
##     end
##     return td
## end

############
## chkElw ##
############

function chkElw(f::Function, tdei::TdEntryIterator)
    nObs, nVars = size(tdei.td)
    res = preallocateBoolDf(names(tdei.td), nObs)
    indRow = 1
    indCol = 1
    for entry in tdei
        res[indRow, indCol] = get(f(entry), :)
        if indRow == nObs
            indRow = 1
            indCol += 1
        else
            indRow += 1
        end
    end
    return Timedata(res, idx(tdei.td))
end

function chkElw(f::Function, tdoi::TdEntryIterator)
    nObs, nVars = size(tdoi.td)
    res = preallocateBoolDf(names(tdoi.td), nObs)
    indRow = 1
    indCol = 1
    for obs in tdoi
        res[indRow, indCol] = get(f(obs), :)
        if indRow == nObs
            indRow = 1
            indCol += 1
        else
            indRow += 1
        end
    end
    return Timedata(res, idx(tdoi.td))
end

## chkVars
##--------

function chkVars(f::Function, tdci::TdColumnIterator)
    res = preallocateBoolDf(names(tdci.td), 1)
    for (nam, col) in tdci
        res[nam] = get(f(col), :)
    end
    return res
end

function chkVars(f::Function, tdvi::TdVariableIterator)
    res = preallocateBoolDf(names(tdvi.td), 1)
    for tdCol in tdvi
        res[names(tdCol)] = get(f(tdCol), :)
    end
    return res
end

## chkDates
##---------

function chkDates(f::Function, tdri::TdRowIterator)
    nObs = size(tdri.td, 1)
    res = DataArray(Bool, nObs)
    indCounter = 1
    for df in tdri
        res[indCounter] = get(f(df), :)
        indCounter += 1
    end
    return Timedata(DataFrame(test_is = res), idx(tdri.td))
end

function chkDates(f::Function, tddi::TdDateIterator)
    nObs = size(tddi.td, 1)
    res = DataArray(Bool, nObs)
    indCounter = 1
    for tdRow in tddi
        res[indCounter] = get(f(tdRow), :)
        indCounter += 1
    end
    return Timedata(DataFrame(test_is = res), idx(tddi.td))
end

###################
## collapseDates ##
###################

## reduces columns to single value
## output: DataFrame

function collapseDates(f::Function, tdci::TdColumnIterator)
    df = DataFrame()
    for (nam, col) in tdci
        df[nam] = f(col)
    end
    return df
end

function collapseDates(f::Function, tdvi::TdVariableIterator)
    df = DataFrame()
    for tdCol in tdvi
        tdRes = f(tdCol)
        df[names(tdCol)[1]] = get(tdRes, :) # get single value
    end
    return df
end

##################
## collapseVars ##
##################

## reduces row to single value
## output: TimeData object
function collapseVars(f::Function, tdri::TdRowIterator)
    nObs = size(tdri.td, 1)
    vals = Array(Any, nObs)
    rowCounter = 1
    for row in tdri # DataFrame
        rowRes = f(row)
        vals[rowCounter, 1] = get(rowRes, :) # get single value
        rowCounter += 1
    end
    df = DataFrame(funcVal = anyToDa(vals))
    return Timedata(df, idx(tdri.td))
end

function collapseVars(f::Function, tddi::TdDateIterator)
    nObs = size(tddi.td, 1)
    vals = Array(Any, nObs)
    rowCounter = 1
    for tmRow in tddi # TimeData object
        rowRes = f(tmRow)
        vals[rowCounter, 1] = get(rowRes, :) # get single value
        rowCounter += 1
    end
    df = DataFrame(funcVal = anyToDa(vals))
    return Timedata(df, idx(tddi.td))
end


text:
function initFind(ed::Editor)
    ed.params[:find] = Dict{Symbol, Any}()

    # 1 indexed, 0 none
    ed.params[:find][:last_match] = 0
    ed.params[:find][:direction] = 1
end

function findCallback(ed::Editor, query::String, key::Char)
    # If the params have not been created, init them
    if !(:find in keys(ed.params))
        initFind(ed)
    end

    last_match = ed.params[:find][:last_match]
    direction = ed.params[:find][:direction]

    if key == '\r' || key == '\x1b'
        last_match = 0
        direction = 1
        return
    elseif key == ARROW_RIGHT || key == ARROW_DOWN
        direction = 1
    elseif key == ARROW_LEFT || key == ARROW_UP
        direction = -1
    else
        last_match = 0
        direction = 1
    end

    last_match == 0 && (direction = 1)
    current = last_match
    for i = 1:length(ed.rows)
        current += direction

        # Bounds check
        if current == 0
            # At begenning of document? Go to end
            current = length(ed.rows)
        elseif current == length(ed.rows)+1
            # At end of doc? Got to start
            current = 1
        end

        row = ed.rows[current]
        loc = search(row.chars, query)
        if loc != 0:-1
            last_match = current
            ed.csr.y = current
            ed.csr.x = first(loc)#charX(row, first(loc))
            ed.rowoff = length(ed.rows)
            break
        end
    end

    # Update params
    ed.params[:find][:last_match] = last_match
    ed.params[:find][:direction] = direction

end

function editorFind(ed::Editor, str::String)
    saved_cx, saved_cy = ed.csr.x, ed.csr.y
    saved_coloff = ed.coloff
    saved_rowoff = ed.rowoff

    findCallback(ed, str, ' ')
    query = editorPrompt(ed, "Search (arrow keys for next/prev): ",
                         callback=findCallback,
                         buf=str,
                         showcursor=false)

    # If the user cancels the search, restore view
    if query == ""
        ed.csr.x = saved_cx
        ed.csr.y = saved_cy
        ed.coloff = saved_coloff
        ed.rowoff = saved_rowoff
    end
end

addCommand(:find, editorFind,
    help="type to find, <enter> to get cursor, arrow keys to next/prev")


text:
include("NfRel/NfRel.jl")
include("NfRel/NfRelNS.jl")
include("NfRel/NfRelNSMor.jl")
include("NfRel/NEQ.jl")
include("NfRel/Conjugates.jl")
include("NfRel/ComplexEmbeddings.jl")
include("NfRel/Simplify.jl")
include("NfRel/absolute_field.jl")
include("NfRel/relative_extension.jl")


text:
module Parsers

export PosLen

using Dates

include("utils.jl")

struct Format
    tokens::Vector{Dates.AbstractDateToken}
    locale::Dates.DateLocale
end

const PtrLen = Tuple{Ptr{UInt8}, Int}

const SupportedFloats = Union{Float16, Float32, Float64, BigFloat}
const SupportedTypes = Union{Integer, SupportedFloats, Dates.TimeType, Bool}

"""
    Parsers.Result{T}(code::Parsers.ReturnCode, tlen[, val])

Struct for holding the results of a parsing operations, returned from `Parsers.xparse`.
Contains 3 fields:
  * `code`: parsing bitmask with various flags set based on parsing (see [ReturnCode](@ref))
  * `tlen`: the total number of bytes consumed while parsing
  * `val`: the parsed value; note that `val` is optional when constructing; users *MUST* check `!Parsers.invalid(result.code)` before accessing `result.val`; if parsing fails, `result.val` is undefined
"""
struct Result{T}
    code::ReturnCode
    tlen::Int64
    val::T
    Result{T}(code::ReturnCode, tlen::Integer) where {T} = new{T}(code, tlen)
    Result{T}(code::ReturnCode, tlen::Integer, val) where {T} = new{T}(code, tlen, val)
end

# Result(code::ReturnCode, tlen::Integer, val::T) where {T} = Result{T}(code, tlen, val)

"""
    `Parsers.Options` is a structure for holding various parsing settings when calling `Parsers.parse`, `Parsers.tryparse`, and `Parsers.xparse`. They include:

  * `sentinel=nothing`: valid values include: `nothing` meaning don't check for sentinel values; `missing` meaning an "empty field" should be considered a sentinel value; or a `Vector{String}` of the various string values that should each be checked as a sentinel value. Note that sentinels will always be checked longest to shortest, with the longest valid match taking precedence.
  * `wh1=' '`: the first ascii character to be considered when ignoring leading/trailing whitespace in value parsing
  * `wh2='\t'`: the second ascii character to be considered when ignoring leading/trailing whitespace in value parsing
  * `openquotechar='"'`: the ascii character that signals a "quoted" field while parsing; subsequent characters will be treated as non-significant until a valid `closequotechar` is detected
  * `closequotechar='"'`: the ascii character that signals the end of a quoted field
  * `escapechar='"'`: an ascii character used to "escape" a `closequotechar` within a quoted field
  * `delim=nothing`: if `nothing`, no delimiter will be checked for; if a `Char` or `String`, a delimiter will be checked for directly after parsing a value or `closequotechar`; a newline (`\n`), return (`\r`), or CRLF (`"\r\n"`) are always considered "delimiters", in addition to EOF
  * `decimal='.'`: an ascii character to be used when parsing float values that separates a decimal value
  * `trues=nothing`: if `nothing`, `Bool` parsing will only check for the string `true` or an `Integer` value of `1` as valid values for `true`; as a `Vector{String}`, each string value will be checked to indicate a valid `true` value
  * `falses=nothing`: if `nothing`, `Bool` parsing will only check for the string `false` or an `Integer` value of `0` as valid values for `false`; as a `Vector{String}`, each string value will be checked to indicate a valid `false` value
  * `dateformat=nothing`: if `nothing`, `Date`, `DateTime`, and `Time` parsing will use a default `Dates.DateFormat` object while parsing; a `String` or `Dates.DateFormat` object can be provided for custom format parsing
  * `ignorerepeated=false`: if `true`, consecutive delimiter characters or strings will be consumed until a non-delimiter is encountered; if `false`, only a single delimiter character/string will be consumed. Useful for fixed-width delimited files where fields are padded with delimiters
  * `quoted=false`: whether parsing should check for `openquotechar` and `closequotechar` characters to signal quoted fields
  * `comment=nothing`: a string which, if matched at the start of a line, will make parsing consume the rest of the line
  * `ignoreemptylines=false`: after parsing a value, if a newline is detected, another immediately proceeding newline will be checked for and consumed
  * `debug=false`: if `true`, various debug logging statements will be printed while parsing; useful when diagnosing why parsing returns certain `Parsers.ReturnCode` values
"""
struct Options
    refs::Vector{String} # for holding references to sentinel, trues, falses, cmt strings
    sentinel::Union{Nothing, Missing, Vector{PtrLen}}
    ignorerepeated::Bool
    ignoreemptylines::Bool
    wh1::UInt8
    wh2::UInt8
    quoted::Bool
    oq::UInt8
    cq::UInt8
    e::UInt8
    delim::Union{Nothing, UInt8, PtrLen}
    decimal::UInt8
    trues::Union{Nothing, Vector{PtrLen}}
    falses::Union{Nothing, Vector{PtrLen}}
    dateformat::Union{Nothing, Format}
    cmt::Union{Nothing, PtrLen}
end

prepare(x::Vector{String}) = sort!(map(ptrlen, x), by=x->x[2], rev=true)
asciival(c::Char) = isascii(c)
asciival(b::UInt8) = b < 0x80

function Options(
            sentinel::Union{Nothing, Missing, Vector{String}},
            wh1::Union{UInt8, Char},
            wh2::Union{UInt8, Char},
            oq::Union{UInt8, Char},
            cq::Union{UInt8, Char},
            e::Union{UInt8, Char},
            delim::Union{Nothing, UInt8, Char, String},
            decimal::Union{UInt8, Char},
            trues::Union{Nothing, Vector{String}},
            falses::Union{Nothing, Vector{String}},
            dateformat::Union{Nothing, String, Dates.DateFormat, Format},
            ignorerepeated, ignoreemptylines, comment, quoted, debug)
    asciival(wh1) && asciival(wh2) || throw(ArgumentError("whitespace characters must be ASCII"))
    asciival(oq) && asciival(cq) && asciival(e) || throw(ArgumentError("openquotechar, closequotechar, and escapechar must be ASCII characters"))
    (oq == delim) || (cq == delim) || (e == delim) && throw(ArgumentError("delim argument must be different than openquotechar, closequotechar, and escapechar arguments"))
    if sentinel isa Vector{String}
        for sent in sentinel
            if startswith(sent, string(Char(wh1))) || startswith(sent, string(Char(wh2)))
                throw(ArgumentError("sentinel value isn't allowed to start with wh1 or wh2 characters"))
            end
            if startswith(sent, string(Char(oq))) || startswith(sent, string(Char(cq)))
                throw(ArgumentError("sentinel value isn't allowed to start with openquotechar, closequotechar, or escapechar characters"))
            end
            if (delim isa UInt8 || delim isa Char) && startswith(sent, string(Char(delim)))
                throw(ArgumentError("sentinel value isn't allowed to start with a delimiter character"))
            elseif delim isa String && startswith(sent, delim)
                throw(ArgumentError("sentinel value isn't allowed to start with a delimiter string"))
            end
        end
        refs = copy(sentinel)
    else
        refs = [""]
    end
    sent = sentinel === nothing || sentinel === missing ? sentinel : prepare(sentinel)
    del = delim === nothing ? nothing : delim isa String ? ptrlen(delim) : delim % UInt8
    if del isa UInt8
        ((wh1 % UInt8) == del || (wh2 % UInt8) == del) && throw(ArgumentError("whitespace characters (`wh1=' '` and `wh2='\\t'` default keyword arguments) must be different than delim argument"))
    end
    if trues !== nothing
        append!(refs, trues)
        trues = prepare(trues)
    end
    if falses !== nothing
        append!(refs, falses)
        falses = prepare(falses)
    end
    if comment === nothing
        cmt = nothing
    else
        push!(refs, comment)
        cmt = ptrlen(comment)
    end
    df = dateformat === nothing ? nothing : dateformat isa String ? Format(dateformat) : dateformat isa Dates.DateFormat ? Format(dateformat) : dateformat
    return Options(refs, sent, ignorerepeated, ignoreemptylines, wh1 % UInt8, wh2 % UInt8, quoted, oq % UInt8, cq % UInt8, e % UInt8, del, decimal % UInt8, trues, falses, df, cmt)
end

Options(;
    sentinel::Union{Nothing, Missing, Vector{String}}=nothing,
    wh1::Union{UInt8, Char}=UInt8(' '),
    wh2::Union{UInt8, Char}=UInt8('\t'),
    openquotechar::Union{UInt8, Char}=UInt8('"'),
    closequotechar::Union{UInt8, Char}=UInt8('"'),
    escapechar::Union{UInt8, Char}=UInt8('"'),
    delim::Union{Nothing, UInt8, Char, String}=nothing,
    decimal::Union{UInt8, Char}=UInt8('.'),
    trues::Union{Nothing, Vector{String}}=nothing,
    falses::Union{Nothing, Vector{String}}=nothing,
    dateformat::Union{Nothing, String, Dates.DateFormat, Format}=nothing,
    ignorerepeated::Bool=false,
    ignoreemptylines::Bool=false,
    comment::Union{Nothing, String}=nothing,
    quoted::Bool=false,
    debug::Bool=false,
) = Options(sentinel, wh1, wh2, openquotechar, closequotechar, escapechar, delim, decimal, trues, falses, dateformat, ignorerepeated, ignoreemptylines, comment, quoted, debug)

const OPTIONS = Options(nothing, UInt8(' '), UInt8('\t'), UInt8('"'), UInt8('"'), UInt8('"'), nothing, UInt8('.'), nothing, nothing, nothing, false, false, nothing, false, false)
const XOPTIONS = Options(missing, UInt8(' '), UInt8('\t'), UInt8('"'), UInt8('"'), UInt8('"'), UInt8(','), UInt8('.'), nothing, nothing, nothing, false, false, nothing, true, false)

# high-level convenience functions like in Base
"""
    Parsers.parse(T, source[, options, pos, len]) => T

Parse a value of type `T` from `source`, which may be a byte buffer (`AbstractVector{UInt8}`), string, or `IO`. Optional arguments include `options`, a [`Parsers.Options`](@ref)
struct, `pos` which indicates the byte position where parsing should begin in a byte buffer `source`, and `len` which is the byte position where parsing should stop in a byte
buffer `source`. If parsing fails for any reason, either invalid value or non-value characters encountered before/after a value, an error will be thrown. To return `nothing`
instead of throwing an error, use [`Parsers.tryparse`](@ref).
"""
function parse(::Type{T}, buf::Union{AbstractVector{UInt8}, AbstractString, IO}, options=OPTIONS, pos::Integer=1, len::Integer=buf isa IO ? 0 : sizeof(buf)) where {T}
    res = xparse2(T, buf isa AbstractString ? codeunits(buf) : buf, pos, len, options)
    code = res.code
    tlen = res.tlen
    fin = buf isa IO || (tlen == (len - pos + 1))
    return ok(code) && fin ? (res.val::T) : throw(Error(buf, T, code, pos, tlen))
end

"""
    Parsers.tryparse(T, source[, options, pos, len]) => Union{T, Nothing}

Parse a value of type `T` from `source`, which may be a byte buffer (`AbstractVector{UInt8}`), string, or `IO`. Optional arguments include `options`, a [`Parsers.Options`](@ref)
struct, `pos` which indicates the byte position where parsing should begin in a byte buffer `source`, and `len` which is the byte position where parsing should stop in a byte
buffer `source`. If parsing fails for any reason, either invalid value or non-value characters encountered before/after a value, `nothing` will be returned. To instead throw
an error, use [`Parsers.parse`](@ref).
"""
function tryparse(::Type{T}, buf::Union{AbstractVector{UInt8}, AbstractString, IO}, options=OPTIONS, pos::Integer=1, len::Integer=buf isa IO ? 0 : sizeof(buf)) where {T}
    res = xparse2(T, buf isa AbstractString ? codeunits(buf) : buf, pos, len, options)
    fin = buf isa IO || (res.tlen == (len - pos + 1))
    return ok(res.code) && fin ? (res.val::T) : nothing
end

"""
    Parsers.xparse(T, buf, pos, len, options) => Parsers.Result{T}

The core parsing function for any type `T`. Takes a `buf`, which can be an `AbstractVector{UInt8}`, `AbstractString`,
or an `IO`. `pos` is the byte position to begin parsing at. `len` is the total # of bytes in `buf` (signaling eof).
`options` is an instance of `Parsers.Options`.

A [`Parsers.Result`](@ref) struct is returned, with the following fields:
      * `res.val` is a value of type `T`, but only if parsing succeeded; for parsed `String`s, no string is returned to avoid excess allocating; if you'd like the actual parsed string value, you can call [`Parsers.getstring`](@ref)
      * `res.code` is a bitmask of parsing codes, use `Parsers.codes(code)` or `Parsers.text(code)` to see the various bit values set. See [`Parsers.ReturnCode`](@ref) for additional details on the various parsing codes
      * `res.tlen`: the total # of bytes consumed while parsing a value, including any quote or delimiter characters; this can be added to the starting `pos` to allow calling `Parsers.xparse` again for a subsequent field/value
"""
function xparse end

# for testing purposes only, it's much too slow to dynamically create Options for every xparse call
function xparse(::Type{T}, buf::Union{AbstractVector{UInt8}, AbstractString, IO}; pos::Integer=1, len::Integer=buf isa IO ? 0 : sizeof(buf), sentinel=nothing, wh1::Union{UInt8, Char}=UInt8(' '), wh2::Union{UInt8, Char}=UInt8('\t'), quoted::Bool=true, openquotechar::Union{UInt8, Char}=UInt8('"'), closequotechar::Union{UInt8, Char}=UInt8('"'), escapechar::Union{UInt8, Char}=UInt8('"'), ignorerepeated::Bool=false, ignoreemptylines::Bool=false, delim::Union{UInt8, Char, PtrLen, AbstractString, Nothing}=UInt8(','), decimal::Union{UInt8, Char}=UInt8('.'), comment=nothing, trues=nothing, falses=nothing, dateformat::Union{Nothing, String, Dates.DateFormat}=nothing, debug::Bool=false) where {T}
    options = Options(sentinel, wh1, wh2, openquotechar, closequotechar, escapechar, delim, decimal, trues, falses, dateformat, ignorerepeated, ignoreemptylines, comment, quoted, debug)
    return xparse(T, buf, pos, len, options)
end

xparse(::Type{T}, buf::AbstractString, pos, len, options=Parsers.XOPTIONS) where {T} =
    xparse(T, codeunits(buf), pos, len, options)

# generic fallback calls Base.tryparse
function xparse(::Type{T}, source, pos, len, options, ::Type{S}=T) where {T, S}
    res = xparse(String, source, pos, len, options)
    code = res.code
    poslen = res.val
    if !Parsers.invalid(code) && !Parsers.sentinel(code)
        str = getstring(source, poslen, options.e)
        x = Base.tryparse(T, str)
        if x === nothing
            return Result{S}(code | INVALID, res.tlen)
        else
            return Result{S}(code, res.tlen, x)
        end
    else
        return Result{S}(code, res.tlen)
    end
end

function xparse(::Type{Char}, source, pos, len, options, ::Type{S}=Char) where {S}
    res = xparse(String, source, pos, len, options)
    code = res.code
    poslen = res.val
    if !Parsers.invalid(code) && !Parsers.sentinel(code)
        return Result{S}(code, res.tlen, first(getstring(source, poslen, options.e)))
    else
        return Result{S}(code, res.tlen)
    end
end

function xparse(::Type{Symbol}, source, pos, len, options, ::Type{S}=Symbol) where {S}
    res = xparse(String, source, pos, len, options)
    code = res.code
    poslen = res.val
    if !Parsers.invalid(code) && !Parsers.sentinel(code)
        if source isa AbstractVector{UInt8}
            sym = ccall(:jl_symbol_n, Ref{Symbol}, (Ptr{UInt8}, Int), pointer(source, poslen.pos), poslen.len)
        else
            sym = Symbol(getstring(source, poslen, options.e))
        end
        return Result{S}(code, res.tlen, sym)
    else
        return Result{S}(code, res.tlen)
    end
end

function xparse(::Type{T}, source::Union{AbstractVector{UInt8}, IO}, pos, len, options::Options=XOPTIONS, ::Type{S}=T) where {T <: SupportedTypes, S}
    startpos = vstartpos = pos
    sentinel = options.sentinel
    wh1 = options.wh1; wh2 = options.wh2
    code = SUCCESS
    quoted = false
    sentinelpos = 0
    if eof(source, pos, len)
        code = (sentinel === missing ? SENTINEL : INVALID) | EOF
        @goto donedone
    end
    b = peekbyte(source, pos)
    # strip leading whitespace
    while b == wh1 || b == wh2
        pos += 1
        incr!(source)
        if eof(source, pos, len)
            code = INVALID | EOF
            @goto donedone
        end
        b = peekbyte(source, pos)
    end
    # check for start of quoted field
    if options.quoted
        quoted = b == options.oq
        if quoted
            code = QUOTED
            pos += 1
            vstartpos = pos
            incr!(source)
            if eof(source, pos, len)
                code |= INVALID_QUOTED_FIELD
                @goto donedone
            end
            b = peekbyte(source, pos)
            # ignore whitespace within quoted field
            while b == wh1 || b == wh2
                pos += 1
                incr!(source)
                if eof(source, pos, len)
                    code |= INVALID_QUOTED_FIELD | EOF
                    @goto donedone
                end
                b = peekbyte(source, pos)
            end
        end
    end
    # check for sentinel values if applicable
    if sentinel !== nothing && sentinel !== missing
        sentinelpos = checksentinel(source, pos, len, sentinel)
    end
    x, code, pos = typeparser(T, source, pos, len, b, code, options)
    if sentinel !== nothing && sentinel !== missing && sentinelpos >= pos
        # if we matched a sentinel value that was as long or longer than our type value
        code &= ~(OK | INVALID | OVERFLOW)
        pos = sentinelpos
        fastseek!(source, pos - 1)
        code |= SENTINEL
        if eof(source, pos, len)
            code |= EOF
        end
    elseif sentinel === missing && pos == vstartpos
        code &= ~(OK | INVALID)
        code |= SENTINEL
    end
    if (code & EOF) == EOF
        if quoted
            # if we detected a quote character, it's an invalid quoted field due to eof in the middle
            code |= INVALID_QUOTED_FIELD
        end
        @goto donedone
    end

@label donevalue
    b = peekbyte(source, pos)
    # donevalue means we finished parsing a value or sentinel, but didn't reach len, b is still the current byte
    # strip trailing whitespace
    while b == wh1 || b == wh2
        pos += 1
        incr!(source)
        if eof(source, pos, len)
            code |= EOF
            if quoted
                code |= INVALID_QUOTED_FIELD
            end
            @goto donedone
        end
        b = peekbyte(source, pos)
    end
    if options.quoted
        # for quoted fields, find the closing quote character
        # we should be positioned at the correct place to find the closing quote character if everything is as it should be
        # if we don't find the quote character immediately, something's wrong, so mark INVALID
        if quoted
            same = options.cq == options.e
            first = true
            while true
                pos += 1
                incr!(source)
                if same && b == options.e
                    if eof(source, pos, len)
                        code |= EOF
                        if !first
                            code |= INVALID
                        end
                        @goto donedone
                    elseif peekbyte(source, pos) != options.cq
                        if !first
                            code |= INVALID
                        end
                        break
                    end
                    code |= ESCAPED_STRING
                    pos += 1
                    incr!(source)
                elseif b == options.e
                    if eof(source, pos, len)
                        code |= INVALID_QUOTED_FIELD | EOF
                        @goto donedone
                    end
                    code |= ESCAPED_STRING
                    pos += 1
                    incr!(source)
                elseif b == options.cq
                    if !first
                        code |= INVALID
                    end
                    if eof(source, pos, len)
                        code |= EOF
                        @goto donedone
                    end
                    break
                end
                if eof(source, pos, len)
                    code |= INVALID_QUOTED_FIELD | EOF
                    @goto donedone
                end
                first = false
                b = peekbyte(source, pos)
            end
            b = peekbyte(source, pos)
            # ignore whitespace after quoted field
            while b == wh1 || b == wh2
                pos += 1
                incr!(source)
                if eof(source, pos, len)
                    code |= EOF
                    @goto donedone
                end
                b = peekbyte(source, pos)
            end
        end
    end

    if options.delim !== nothing
        delim = options.delim
        # now we check for a delimiter; if we don't find it, keep parsing until we do
        if !options.ignorerepeated
            # we're checking for a single appearance of a delimiter
            if delim isa UInt8
                if b == delim
                    pos += 1
                    incr!(source)
                    code |= DELIMITED
                    @goto donedone
                end
            elseif delim isa PtrLen
                predelimpos = pos
                pos = checkdelim(source, pos, len, delim)
                if pos > predelimpos
                    # found the delimiter we were looking for
                    code |= DELIMITED
                    @goto donedone
                end
            else
                error()
            end
        else
            # keep parsing as long as we keep matching delims/newlines
            if delim isa UInt8
                matched = false
                matchednewline = false
                while true
                    if b == delim
                        matched = true
                        code |= DELIMITED
                        pos += 1
                        incr!(source)
                    elseif !matchednewline && b == UInt8('\n')
                        matchednewline = matched = true
                        pos += 1
                        incr!(source)
                        pos = checkcmtemptylines(source, pos, len, options)
                        code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
                    elseif !matchednewline && b == UInt8('\r')
                        matchednewline = matched = true
                        pos += 1
                        incr!(source)
                        if !eof(source, pos, len) && peekbyte(source, pos) == UInt8('\n')
                            pos += 1
                            incr!(source)
                        end
                        pos = checkcmtemptylines(source, pos, len, options)
                        code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
                    else
                        break
                    end
                    if eof(source, pos, len)
                        @goto donedone
                    end
                    b = peekbyte(source, pos)
                end
                if matched
                    @goto donedone
                end
            elseif delim isa PtrLen
                matched = false
                matchednewline = false
                while true
                    predelimpos = pos
                    pos = checkdelim(source, pos, len, delim)
                    if pos > predelimpos
                        matched = true
                        code |= DELIMITED
                    elseif !matchednewline && b == UInt8('\n')
                        matchednewline = matched = true
                        pos += 1
                        incr!(source)
                        pos = checkcmtemptylines(source, pos, len, options)
                        code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
                    elseif !matchednewline && b == UInt8('\r')
                        matchednewline = matched = true
                        pos += 1
                        incr!(source)
                        if !eof(source, pos, len) && peekbyte(source, pos) == UInt8('\n')
                            pos += 1
                            incr!(source)
                        end
                        pos = checkcmtemptylines(source, pos, len, options)
                        code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
                    else
                        break
                    end
                    if eof(source, pos, len)
                        @goto donedone
                    end
                    b = peekbyte(source, pos)
                end
                if matched
                    @goto donedone
                end
            else
                error()
            end
        end
        # didn't find delimiter, but let's check for a newline character
        if b == UInt8('\n')
            pos += 1
            incr!(source)
            pos = checkcmtemptylines(source, pos, len, options)
            code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
            @goto donedone
        elseif b == UInt8('\r')
            pos += 1
            incr!(source)
            if !eof(source, pos, len) && peekbyte(source, pos) == UInt8('\n')
                pos += 1
                incr!(source)
            end
            pos = checkcmtemptylines(source, pos, len, options)
            code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
            @goto donedone
        end
        # didn't find delimiter or newline, so we're invalid, keep parsing until we find delimiter, newline, or len
        code |= INVALID_DELIMITER
        while true
            pos += 1
            incr!(source)
            if eof(source, pos, len)
                code |= EOF
                @goto donedone
            end
            b = peekbyte(source, pos)
            if !options.ignorerepeated
                if delim isa UInt8
                    if b == delim
                        pos += 1
                        incr!(source)
                        code |= DELIMITED
                        @goto donedone
                    end
                elseif delim isa PtrLen
                    predelimpos = pos
                    pos = checkdelim(source, pos, len, delim)
                    if pos > predelimpos
                        # found the delimiter we were looking for
                        code |= DELIMITED
                        @goto donedone
                    end
                else
                    error()
                end
            else
                if delim isa UInt8
                    matched = false
                    matchednewline = false
                    while true
                        if b == delim
                            matched = true
                            code |= DELIMITED
                            pos += 1
                            incr!(source)
                        elseif !matchednewline && b == UInt8('\n')
                            matchednewline = matched = true
                            pos += 1
                            incr!(source)
                            pos = checkcmtemptylines(source, pos, len, options)
                            code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
                        elseif !matchednewline && b == UInt8('\r')
                            matchednewline = matched = true
                            pos += 1
                            incr!(source)
                            if !eof(source, pos, len) && peekbyte(source, pos) == UInt8('\n')
                                pos += 1
                                incr!(source)
                            end
                            pos = checkcmtemptylines(source, pos, len, options)
                            code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
                        else
                            break
                        end
                        if eof(source, pos, len)
                            @goto donedone
                        end
                        b = peekbyte(source, pos)
                    end
                    if matched
                        @goto donedone
                    end
                elseif delim isa PtrLen
                    matched = false
                    matchednewline = false
                    while true
                        predelimpos = pos
                        pos = checkdelim(source, pos, len, delim)
                        if pos > predelimpos
                            matched = true
                            code |= DELIMITED
                        elseif !matchednewline && b == UInt8('\n')
                            matchednewline = matched = true
                            pos += 1
                            incr!(source)
                            pos = checkcmtemptylines(source, pos, len, options)
                            code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
                        elseif !matchednewline && b == UInt8('\r')
                            matchednewline = matched = true
                            pos += 1
                            incr!(source)
                            if !eof(source, pos, len) && peekbyte(source, pos) == UInt8('\n')
                                pos += 1
                                incr!(source)
                            end
                            pos = checkcmtemptylines(source, pos, len, options)
                            code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
                        else
                            break
                        end
                        if eof(source, pos, len)
                            @goto donedone
                        end
                        b = peekbyte(source, pos)
                    end
                    if matched
                        @goto donedone
                    end
                else
                    error()
                end
            end
            # didn't find delimiter, but let's check for a newline character
            if b == UInt8('\n')
                pos += 1
                incr!(source)
                pos = checkcmtemptylines(source, pos, len, options)
                code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
                @goto donedone
            elseif b == UInt8('\r')
                pos += 1
                incr!(source)
                if !eof(source, pos, len) && peekbyte(source, pos) == UInt8('\n')
                    pos += 1
                    incr!(source)
                end
                pos = checkcmtemptylines(source, pos, len, options)
                code |= NEWLINE | ifelse(eof(source, pos, len), EOF, SUCCESS)
                @goto donedone
            end
        end
    end

@label donedone
    tlen = pos - startpos
    if valueok(code)
        y::T = x
        return Result{S}(code, tlen, y)
    else
        return Result{S}(code, tlen)
    end
end

# condensed version of xparse that doesn't worry about quoting or delimiters; called from Parsers.parse/Parsers.tryparse
@inline function xparse2(::Type{T}, source::Union{AbstractVector{UInt8}, IO}, pos, len, options::Options=XOPTIONS, ::Type{S}=T) where {T <: SupportedTypes, S}
    startpos = pos
    sentinel = options.sentinel
    wh1 = options.wh1; wh2 = options.wh2
    code = SUCCESS
    quoted = false
    if eof(source, pos, len)
        code = (sentinel === missing ? SENTINEL : INVALID) | EOF
        @goto donedone
    end
    b = peekbyte(source, pos)
    # strip leading whitespace
    while b == wh1 || b == wh2
        pos += 1
        incr!(source)
        if eof(source, pos, len)
            code = INVALID | EOF
            @goto donedone
        end
        b = peekbyte(source, pos)
    end
    x, code, pos = typeparser(T, source, pos, len, b, code, options)
    if (code & EOF) == EOF
        @goto donedone
    end

@label donevalue
    b = peekbyte(source, pos)
    # donevalue means we finished parsing a value or sentinel, but didn't reach len, b is still the current byte
    # strip trailing whitespace
    while b == wh1 || b == wh2
        pos += 1
        incr!(source)
        if eof(source, pos, len)
            code |= EOF
            if quoted
                code |= INVALID_QUOTED_FIELD
            end
            @goto donedone
        end
        b = peekbyte(source, pos)
    end

@label donedone
    tlen = pos - startpos
    return valueok(code) ? Result{S}(code, tlen, x) : Result{S}(code, tlen)
end

@inline function xparse2(::Type{T}, source, pos, len, options, ::Type{S}=T) where {T, S}
    res = xparse(String, source, pos, len, options)
    code = res.code
    poslen = res.val
    if !Parsers.invalid(code) && !Parsers.sentinel(code)
        str = getstring(source, poslen, options.e)
        x = Base.tryparse(T, str)
        if x === nothing
            return Result{S}(code | INVALID, res.tlen)
        else
            return Result{S}(code, res.tlen, x)
        end
    else
        return Result{S}(code, res.tlen)
    end
end

@inline function xparse2(::Type{Char}, source, pos, len, options, ::Type{S}=Char) where {S}
    res = xparse(String, source, pos, len, options)
    code = res.code
    poslen = res.val
    if !Parsers.invalid(code) && !Parsers.sentinel(code)
        return Result{S}(code, res.tlen, first(getstring(source, poslen, options.e)))
    else
        return Result{S}(code, res.tlen)
    end
end

@inline function xparse2(::Type{Symbol}, source, pos, len, options, ::Type{S}=Symbol) where {S}
    res = xparse(String, source, pos, len, options)
    code = res.code
    poslen = res.val::PosLen
    if !Parsers.invalid(code) && !Parsers.sentinel(code)
        if source isa AbstractVector{UInt8}
            sym = ccall(:jl_symbol_n, Ref{Symbol}, (Ptr{UInt8}, Int), pointer(source, poslen.pos), poslen.len)
        else
            sym = Symbol(getstring(source, poslen, options.e))
        end
        return Result{S}(code, res.tlen, sym)
    else
        return Result{S}(code, res.tlen)
    end
end

function checkdelim!(buf, pos, len, options::Options)
    pos > len && return pos
    delim = options.delim
    @inbounds b = buf[pos]
    valuepos = pos
    if !options.ignorerepeated
        # we're checking for a single appearance of a delimiter
        if delim isa UInt8
            b == delim && return pos + 1
        else
            pos = checkdelim(buf, pos, len, delim)
            pos > valuepos && return pos
        end
    else
        # keep parsing as long as we keep matching delim
        if delim isa UInt8
            matched = false
            while b == delim
                matched = true
                pos += 1
                pos > len && return pos
                @inbounds b = buf[pos]
            end
            matched && return pos
        elseif delim isa PtrLen
            matched = false
            predelimpos = pos
            pos = checkdelim(buf, pos, len, delim)
            while pos > predelimpos
                matched = true
                pos > len && return pos
                predelimpos = pos
                pos = checkdelim(buf, pos, len, delim)
            end
            matched && return pos
        else
            error()
        end
    end
    return pos
end

function checkcmtemptylines(source, pos, len, options::Options)
    cmt = options.cmt
    while !eof(source, pos, len)
        skipped = false
        if options.ignoreemptylines
            b = peekbyte(source, pos)
            if b == UInt8('\n')
                pos += 1
                incr!(source)
                skipped = true
            elseif b == UInt8('\r')
                pos += 1
                incr!(source)
                if !eof(source, pos, len) && peekbyte(source, pos) == UInt8('\n')
                    pos += 1
                    incr!(source)
                end
                skipped = true
            end
        end
        matched = false
        if cmt !== nothing && !eof(source, pos, len)
            precmtpos = pos
            pos = checkdelim(source, pos, len, cmt)
            if pos > precmtpos
                matched = true
                eof(source, pos, len) && break
                b = peekbyte(source, pos)
                while b != UInt8('\n') && b != UInt8('\r')
                    pos += 1
                    incr!(source)
                    eof(source, pos, len) && break
                    b = peekbyte(source, pos)
                end
                pos += 1
                incr!(source)
                if b == UInt8('\r') && !eof(source, pos, len) && peekbyte(source, pos) == UInt8('\n')
                    pos += 1
                    incr!(source)
                end
            end
        end
        (skipped | matched) || break
    end
    return pos
end

include("ints.jl")
include("floats.jl")
include("strings.jl")
include("bools.jl")
include("dates.jl")
include("ryu.jl")

function __init__()
    resize!(empty!(BIGINT), Threads.nthreads())
    resize!(empty!(BIGFLOAT), Threads.nthreads())
    return
end

include("precompile.jl")

end # module


text:
module SSLTestSetIndex
using SimpleSDMLayers
using Test

S = SimpleSDMResponse(zeros(10, 10), 0.0, 1.0, 0.0, 1.0)

# Change a single value
changepoint = Point(0.2, 0.5)
S[changepoint] = 1.0
@test S[changepoint] == 1.0

# Change multiple by multiple values
changepoints = [Point(rand(2)...) for i in 1:8]
changevalues = rand(length(changepoints))
S[changepoints] = changevalues
# NOTE we test that the values in the layer are part of the pool of possible
# values, because it is entirely possible that two points will fall within the
# same cell
@test all(.!isnothing.(indexin(S[changepoints], changevalues)))

# Change multiple by multiple values as a matrix
changepoints = reshape([Point(rand(2)...) for i in 1:8], (2, 4))
changevalues = reshape(rand(length(changepoints)), (2,4))
S[changepoints] = changevalues
@test all(.!isnothing.(indexin(S[changepoints], changevalues)))

#=

# Change multiple values by a single value
changepoints = [Point(rand(2)...) for i in 1:8]
S[changepoints] .= 5.0
@test all(S[changepoints] .== 5.0)

# Change multiple values by a single value as a matrix
changepoints = reshape([Point(rand(2)...) for i in 1:8], (2, 4))
S[changepoints] .= 5.0
@test all(S[changepoints] .== 5.0)

=#

end


text:
function inner_render!(canvas::Canvas, scene::Geometry, camera::Camera, render_func!::Function, args...)
    w, h = width(canvas), height(canvas)
    for x = 1:h
        sx = x / h
        for y = 1:w
            sy = 1 - y / w
            ray = generate_ray(camera, sx, sy)
            result = ray ∩ scene
            if hit(result)
                render_func!(canvas, x, y, scene, result.value, args...)
            end
        end
    end
    canvas
end

function put_color!(canvas::Canvas, x::Int64, y::Int64, color::ColorRGB)
    canvas[1, y, x] = color.r
    canvas[2, y, x] = color.g
    canvas[3, y, x] = color.b
    canvas[4, y, x] = 1.0
end

function ray_trace_color(scene::Geometry, result::IntersectResult{T}, max_reflect::Int64) where T <: GeometryWithMaterial
    reflectiveness = result.geometry.material.reflectiveness
    color = sample(result.geometry.material, result.ray, result.position, result.normal)
    color *= 1 - reflectiveness
    if reflectiveness > 0 && max_reflect > 0
        new_ray = Ray3(result.position, result.ray.direction ⋅ result.normal * (-2) * result.normal + result.ray.direction)
        new_result = new_ray ∩ scene
        if hit(new_result)
            reflected_color = ray_trace_color(scene, new_result.value, max_reflect - 1)
        else
            reflected_color = Black
        end
        color += reflected_color
    end
    return color
end

function ray_trace_render!(canvas::Canvas, x::Int64, y::Int64, scene::Geometry, result::IntersectResult{T},
    max_reflect::Int64) where T <: GeometryWithMaterial
    color = ray_trace_color(scene, result, max_reflect)
    put_color!(canvas, x, y, color)
end

function light_render!(canvas::Canvas, x::Int64, y::Int64, scene::Geometry, result::IntersectResult{T},
    lights::Lights) where T
    color = Black
    for light in lights
        light_sample = sample(light, scene, result.position)
        if !in_shadow(light_sample)
            ndotl = result.normal ⋅ light_sample.light_vector
            if ndotl >= 0
                color += light_sample.irradiance * ndotl
            end
        end
    end
    put_color!(canvas, x, y, color)
end

function render!(canvas::Canvas, scene::Geometry, camera::Camera, max_reflect::Int64=3)
    inner_render!(canvas, scene, camera, ray_trace_render!, max_reflect)
end


text:
x, y = generators(MPoly{Float64}, :x, :y)

p = 4.0 + x^3 + 3.0*x*y^2 + y^3
q =  x^(-2)*y^2
r = x^(-1) * y

@test diff(p, :x, 1) == 3.0x^2 + 3.0y^2
@test diff(p, :x, 2) == 6.0x
@test diff(p, :x, 3) == 6.0 * one(p)
@test diff(p, :y, 1) == 6x*y + 3y^2
@test diff(p, :z) == zero(p)
@test p == 4.0 + x^3 + 3*x*y^2 + y^3 # Check p is left alone

@test diff(q, :x) == -2*x^(-3)*y^2
@test diff(q, :y) == 2*x^(-2)*y

p2 = x
@test diff(p2, :x, 1) == 1.0 * one(p)
@test diff(p2, :y, 1) == zero(p)

p1 = 3.0 + 3x + y^2
@test integrate(p1, :x, 1) == 3.0x + 3*x^2 * 0.5 + x*y^2
@test integrate(p1, :y, 2) == 1.5*y^2 + 1.5x*y^2 + (1/3 * 1/4) * y^4
pz = MPoly{Float64}(:z)
@test integrate(p1, :z) == pz * p1
@test integrate(p1, :z, 2) == pz^2 * p1 / 2.0
@test p1 == 3.0 + 3x + y^2 # Check p is left alone

@test integrate(q, :x) == -x^(-1)*y^2
@test integrate(r, :y) == x^(-1)*y^2/2
@test_throws ArgumentError integrate(q, :x, 2)
@test_throws ArgumentError integrate(r, :x)


text:
using SVGDraw
using Test

@testset "SVGDraw.jl" begin
    # Write your tests here.
end


text:
abstract type AbstractImageNoise <: AbstractImageFilter end
abstract type AbstractImageDenoiseAlgorithm <: AbstractImageFilter end

@filter_api apply_noise AbstractImageNoise
@filter_api reduce_noise AbstractImageDenoiseAlgorithm

"""
    apply_noise([::Type,] img, n::AbstractImageNoise, args...; rng=GLOBAL_RNG, kwargs...)

Add/Apply noise `n` to image `img`.

# Examples

```julia
n = AdditiveWhiteGaussianNoise(0.1)
img = testimage("lena_gray_256")
apply_noise(img, n)

# sometimes we need to pass `rng` keyword argument to
# generate reproducible noise
apply_noise(img, n; rng = MersenneTwister(0))
```

See also: [`apply_noise!`](@ref apply_noise!)
"""
apply_noise

"""
    apply_noise!([out,] img, n::AbstractImageNoise, args...; rng=GLOBAL_RNG, kwargs...)

Add/Apply noise `n` to image `img`.

If `out` is specified, it will be changed in place. Otherwise `img` will be changed in place.

# Examples

```julia
n = AdditiveWhiteGaussianNoise(0.1)
img = testimage("lena_gray_256")
apply_noise!(img, n)

# sometimes we need to pass `rng` keyword argument to
# generate reproducible noise
apply_noise!(img, n; rng = MersenneTwister(0))
```

See also: [`apply_noise`](@ref apply_noise)
"""
apply_noise!

"""
    reduce_noise([::Type,] img, f::AbstractImageDenoiseAlgorithm, args...)

Remove noise of image `img` using algorithm `f`.

See also: [`reduce_noise!`](@ref reduce_noise!)
"""
reduce_noise

"""
    reduce_noise!([out,] img, f::AbstractImageDenoiseAlgorithm, args...)

Remove noise of image `img` using algorithm `f`.

If `out` is specified, it will be changed in place. Otherwise `img` will be changed in place.

See also: [`reduce_noise`](@ref reduce_noise)
"""
reduce_noise!


text:
module PlusMinus

export plusminus, ±

plusminus(x::Number,y::Number) = (x-y,x+y)

±(x,y) = plusminus(x,y)

end # module