text:
function pitches_gaps!(s, q) #Sometimes in the midifile there are small gaps between pitches, this function removes them.
    for v in s
        for i = 1:(size(v)[1]-1)
            m = abs(v[i,2] - v[i+1,1])
            #printlrintln(Voz[i,2],'\t',Voz[i+1,1],'\t',m)
            if m < q
                v[i,2] = v[i+1,1]
                #println(Voz[i,2],'\t', Voz[i+1,1])
            end
        end
      end
end
#########################################################################################################################################################################
function silence_gaps!(s,q) #Sometimes there are small gaps between silences and pitches, this function removes them.
    sm = q
    for v in s
        for i = 1:(size(v)[1]-1)
            m = mod(v[i,2] - v[i+1,1],sm)
            if m > 0
                v[i,2] = v[i,2] + (sm - m)
            end
        end
        x = length(v[:,1])
        n = mod(v[x,2] - v[x,1], sm)
        if n > 0
            v[x,2] = v[x,2] + (sm - n)
        end
    end
end
#################################################################################################################################################

function min_voces(n_vs, dv) #this function returns the minimum time value of all voices.
    ns = size(n_vs)[1]
    mins = Array{Float64}(undef,ns)
    for i = 1:ns
        dif = n_vs[i][:,2] - n_vs[i][:,1]
        difn = filter(x->x >= dv,dif)
        mins[i] = minimum(difn)
    end
    return minimum(mins)
end

################################################################################################################################################
function max_tempo(Voces, ns) #This function is to find what is the maximum time when the piece ends.
    Tfinal = Array{Float64}(undef,ns)
    for i = 1:ns
        Tfinal[i] = maximum(Voces[i][:,2])
    end
    return maximum(Tfinal)
end
################################################################################################################################
#rounding the possible fractional numbers obtained by the subdivision.
function rounding!(voces, q)
    nv = size(voces)[1]
    for i =1:nv
        voces[i][:,1] = map(x-> ceil(x/q), voces[i][:,1])
        voces[i][:,2] = map(x-> ceil(x/q), voces[i][:,2])
    end
end
####################################################################################################
##########################################################################################################
#Next function filter the undef voices in left in the array of voices...
function filter_undef!(voces)
    nv = size(voces)[1]
    tt = Array{Bool}(undef,nv)
    for i=1:nv
        tt[i] = isassigned(voces,i)
    end
    deleteat!(voces,find(x->x==false,tt))
end
####################################################################################################################
function indice(tmax::Int64) #esta funcion solo regresa el arreglo que lleva el indice(numero) de las notas (o el eje x)
    ind = Array{Float64}(undef,tmax)
    for i = 1:tmax
        ind[i] = convert(Float64, i)
    end
    return(ind)
end
##############################################################################################################################
function notas_hertz!(notas::Array{Float64,1}) #esta funcion convierte el arreglo de entrada de numero de nota a su valor en hertz con A4 = 440Hz
    a4 = 440.0
    x = 1.0/12.0
    r = 2 ^ x
    for i = 1:length(notas)
        if notas[i] == 0.0 #si la nota ya es silencio, la omite
            continue
        end
        d = notas[i] - 69 # 69 es el A4
        notas[i] = a4 * r ^ d
    end
end
################################################################################
function serie_notas(Voz, tmax) #constructs the time series of the pithces start - final - voice
    serie = zeros(tmax)
    ini = Voz[:,1]
    dura = Voz[:,2] - Voz[:,1] #duration
    for i = 1:(length(dura)) #aqui se asignan las notas
        j = 0
        x = convert(Int64, ini[i]) + 1
        while dura[i] > 0
            serie[x+j] = Voz[i,3]
            j += 1
            dura[i] -= 1
        end
    end
    return(serie)
end
###########################################################################################
#Next function is to get the subgroups of size sz in a sequence, subgroups are with overlap.
function get_subs(t_s,sz)
    subs = []
    t = length(t_s)
    for i = 1:(t-sz+1)
        sub = []
        for j = i:(i+sz-1)
            push!(sub,t_s[j])
        end
        push!(subs,sub)
    end

    return subs
end
###################################################################################
#New functions to construct time series and pitch sequences
function get_onoff_notes(s)
    ini = s[s[:,3].==" Note_on_c", :] #getting the the starting places
    fin = s[s[:,3].==" Note_off_c",:] #ending
    dn = unique(ini[:,5]) #different notes.
    ndn = length(dn) #number of different notes
    posin_dn = Array{Vector}(undef, ndn) #initialize arrays for getting information of each note
    posfn_dn = Array{Vector}(undef, ndn)
    in_fi = Array{Matrix}(undef,ndn)
    for i = 1:ndn       #for each different note it gets when are they played
        tmp = ini[ini[:,5].==dn[i],:][:,2]
        tmp2 = fin[fin[:,5].==dn[i],:][:,2]
        l1 = length(tmp); l2 = length(tmp2)
        if l1 != l2
            nt = min(l1,l2)
            note = [dn[i] for j=1:nt]
            in_fi[i] = [tmp[1:nt] tmp2[1:nt] note]
        else
            note = [dn[i] for j = 1:length(tmp)]
            in_fi[i] = [tmp tmp2 note]
        end
    end
    nmat = vcat(in_fi...) #getting all notes in the same array
    return nmat[sortperm(nmat[:,1]),:]  #returns the array sorted by appereance.
end
################################################################################
function get_new_voice(v)
    new_v = []
    ids = []
    iids = [i for i=1:size(v)[1]]
    c = 1
    while c <= size(v)[1]-1
        if v[c+1,1] <= v[c,1]
            #println(v[c,:],'\t',v[c+1,:])
            push!(new_v, v[c+1,:])
            push!(ids,c+1)
            c+=1
        else
            #push!(ids,c)
            c+=1
        end
        #c += 1
        #println(c)
    end
    deleteat!(iids, ids)
    old_v = v[iids,:]
    new_v = convert(Array{Any,2},transpose(reshape(vcat(new_v...),3,:)))
    return [old_v, new_v]
end
function find_more_voices(v)
    v_1 =[]
    nvs = get_new_voice(v)
    push!(v_1,nvs[1])
    lv = nvs[2]
    br = false
    while br == false #this is to check if there are more voices
        ns = get_new_voice(lv)
        push!(v_1, ns[1])
        br = isempty(ns[2])
        if br; break; end
        #println(br)
        lv = ns[2]
    end
    return v_1
end
################################################################################
function get_onon_notes(s)
    ini = s[s[:,6].!=0, :] #getting the the starting places
    fin = s[s[:,6].==0,:] #ending
    dn = unique(ini[:,5]) #different notes.
    ndn = length(dn) #number of different notes
    posin_dn = Array{Vector}(undef, ndn) #initialize arrays for getting information of each note
    posfn_dn = Array{Vector}(undef, ndn)
    in_fi = Array{Matrix}(undef,ndn)
    for i = 1:ndn       #for each different note it gets when are they played
        tmp = ini[ini[:,5].==dn[i],:][:,2]
        tmp2 = fin[fin[:,5].==dn[i],:][:,2]
        l1 = length(tmp); l2 = length(tmp2)
        if l1 != l2
            nt = min(l1,l2)
            note = [dn[i] for j=1:nt]
            in_fi[i] = [tmp[1:nt] tmp2[1:nt] note]
        else
            note = [dn[i] for j = 1:length(tmp)]
            in_fi[i] = [tmp tmp2 note]
        end
    end
    nmat = vcat(in_fi...) #getting all notes in the same array
    return nmat[sortperm(nmat[:,1]),:]  #returns the array sorted by appereance.
end
################################################################################
#next function gets the multivariate time series of a pieces from a csv file
function get_time_series(s)
    sd = 16
    #mq / sd
    #sd is the smallest reference duration, it would divide the quarter of a note.
    mq = s[1,6] #quarter of a note
    nv = s[findlast(x-> x==" Note_on_c",s[:,3]),1] #- 1 #estimates how many voices are in the midi
    voces = Array{Matrix}(undef,nv) #initialze an array
    #next function is to get the time in miliseconds when the pitch starts and ends.
    if nv ==1
        #nv = 1
        if findfirst(x-> x==" Note_off_c",s[:,3]) == 0
            b = s[s[:,1].==1,:] #takes the events of the voice i
            mat = b[b[:,3].==" Note_on_c", :]
            mat = mat[find(x->x!="", mat[:,5]),:]
            voces = get_onon_notes(mat) #construct an array of of information of initial time, finish time, pitch and intensity.
        else
            mat = s[s[:,1].==1,:]
            mat = mat[find(x->x!="", mat[:,5]),:]
            voces = get_onoff_notes(mat)
        end
    else
        voces = Array{Matrix}(undef,nv) #initialze an array
        if findfirst(x-> x==" Note_off_c",s[:,3]) == 0
            #checks if the midi has events of note_off
            for i = 1:(nv) #if does not, it construct the series in this way
                b = s[s[:,1].==i,:] #takes the events of the voice i
                if size(b[b[:,3].==" Note_on_c", :])[1] == 0; continue; end #checks if there are notes in the channel
                mat = b[b[:,3].==" Note_on_c", :]
                mat = mat[find(x->x!="", mat[:,5]),:]
                voces[i] =  get_onon_notes(mat)#construct an array of information of initial time, finish time, pitch and intensity.
            end
        else
            for i = 1:nv #if has note_off events , it does this way
                ini = findfirst(x-> x==" Note_on_c", s[s[:,1].==i,3]) #initial time
                fin = findlast(x-> x==" Note_off_c", s[s[:,1].==i,3]) #finish time
                if ini == 0 || fin == 0; continue; end
                mat = s[s[:,1].==i,:]
                mat = mat[findall(x->x!="", mat[:,5]),:]
                voces[i] = get_onoff_notes(mat) #construct an array of information of initial time, finish time, pitch and intensity
            end
        end
    end
    if nv > 1
        filter_undef!(voces)
        filter!(x->length(x)>0,voces)
        n_vs = []
        nv = length(voces)
        for i=1:nv
            m_v = find_more_voices(voces[i])
            for j =1:length(m_v)
                push!(n_vs, m_v[j])
            end
        end
    else
        try
            n_vs = []
            m_v = find_more_voices(voces[1])
            for j =1:length(m_v)
                push!(n_vs, m_v[j])
            end
        catch
            n_vs = []
            m_v = find_more_voices(voces)
            for j =1:length(m_v)
                push!(n_vs, m_v[j])
            end
        end
    end
    filter!(x->length(x)>0,n_vs)
    nv = size(n_vs)[1] #the real number of voices
    q = round(Int, mq / (mq / min_voces(n_vs, mq / sd)))
    #q = mq/ round(Int, mq / min_voces(n_vs, mq / sd)) #defines the unit of time  q
    pitches_gaps!(n_vs,q) #fixes some gaps of miliseconds between the notes
    silence_gaps!(n_vs,q) #fixes some gaps of milisecons of rests
    rounding!(n_vs, q) #rounds up some possible decimals in the time
    tmax = round(Int,max_tempo(n_vs,nv)) #gets the total lenght of the piece
    series = zeros(tmax,nv+1)  #this is the initialization of the output array
    series[:,1] = indice(tmax) #this is just the index
    for i=2:(nv+1) #constructs the output time series
        series[:,i] = serie_notas(n_vs[i-1],tmax)
    end
    return series
end

################################################################################
#get the sequences of pitches
function get_pitch_seq(s)
    sd = 16
    #mq / sd
    #sd is the smallest reference duration, it would divide the quarter of a note.
    mq = s[1,6] #quarter of a note
    nv = s[findlast(s[:,3], " Note_on_c"),1]# - 1 #estimates how many voices are in the midi
    voces = Array{Matrix}(undef,nv) #initialze an array
    #next function is to get the time in miliseconds when the pitch starts and ends.
    if nv ==1
        #nv = 1
        if findfirst(x-> x==" Note_off_c",s[:,3]) == 0
            b = s[s[:,1].==1,:] #takes the events of the voice i
            mat = b[b[:,3].==" Note_on_c", :]
            mat = mat[find(x->x!="", mat[:,5]),:]
            voces = get_onon_notes(mat) #construct an array of of information of initial time, finish time, pitch and intensity.
        else
            mat = s[s[:,1].==1,:]
            mat = mat[find(x->x!="", mat[:,5]),:]
            voces = get_onoff_notes(mat)
        end
    else
        voces = Array{Matrix}(undef,nv) #initialze an array
        if findfirst(x-> x==" Note_off_c",s[:,3]) == 0
            #checks if the midi has events of note_off
            for i = 1:(nv) #if does not, it construct the series in this way
                b = s[s[:,1].==i,:] #takes the events of the voice i
                if size(b[b[:,3].==" Note_on_c", :])[1] == 0; continue; end #checks if there are notes in the channel
                mat = b[b[:,3].==" Note_on_c", :]
                mat = mat[find(x->x!="", mat[:,5]),:]
                voces[i] =  get_onon_notes(mat)#construct an array of information of initial time, finish time, pitch and intensity.
            end
        else
            for i = 1:nv #if has note_off events , it does this way
                ini = findfirst(x-> x==" Note_on_c", s[s[:,1].==i,3]) #initial time
                fin = findlast(x-> x==" Note_off_c", s[s[:,1].==i,3]) #finish time
                if ini == 0 || fin == 0; continue; end
                mat = s[s[:,1].==i,:]
                mat = mat[findall(x->x!="", mat[:,5]),:]
                voces[i] = get_onoff_notes(mat) #construct an array of information of initial time, finish time, pitch and intensity
            end
        end
    end
    if nv > 1
        filter_undef!(voces)
        filter!(x->length(x)>0,voces)
        n_vs = []
        nv = length(voces)
        for i=1:nv
            m_v = find_more_voices(voces[i])
            for j =1:length(m_v)
                push!(n_vs, m_v[j])
            end
        end
    else
        try
            n_vs = []
            m_v = find_more_voices(voces[1])
            for j =1:length(m_v)
                push!(n_vs, m_v[j])
            end
        catch
            n_vs = []
            m_v = find_more_voices(voces)
            for j =1:length(m_v)
                push!(n_vs, m_v[j])
            end
        end
    end
    filter!(x->length(x)>0,n_vs)
    nv = size(n_vs)[1] #the real number of voices
    pitch_seq = Array{Vector,1}(undef, nv)
    for i = 1:nv
        pitch_seq[i] = n_vs[i][:,3]
    end
    return pitch_seq
end


text:
"""
    build a single histogram containing all grads and weight information
"""
function hist_kernel!(h::CuDeviceArray{T,3}, δ::CuDeviceMatrix{T}, xid::CuDeviceMatrix{S}, 𝑖, 𝑗) where {T,S}

    nbins = size(h, 2)

    it = threadIdx().x
    id, jd = blockDim().x, blockDim().y
    ib, j, k = blockIdx().x, blockIdx().y, blockIdx().z
    ig, jg = gridDim().x, gridDim().y

    shared = @cuDynamicSharedMem(T, nbins)
    fill!(shared, 0)
    sync_threads()

    i_tot = length(𝑖)
    iter = 0
    while iter * id * ig < i_tot
        i = it + id * (ib - 1) + iter * id * ig
        if i <= length(𝑖) && j <= length(𝑗)
            # depends on shared to be assigned to a single feature
            @inbounds i_idx = 𝑖[i]
            @inbounds CUDA.atomic_add!(pointer(shared, xid[i_idx, 𝑗[j]]), δ[i_idx, k])
        end
        iter += 1
    end
    sync_threads()
    # loop to cover cases where nbins > nthreads
    for iter in 1:(nbins - 1) ÷ id + 1
        bin_id = it + id * (iter - 1)
        # bin_id = it
        if bin_id <= nbins
            @inbounds δid = Base._to_linear_index(h, k, bin_id, 𝑗[j])
            @inbounds CUDA.atomic_add!(pointer(h, δid), shared[bin_id])
        end
    end
    sync_threads()
    return nothing
end

# base approach - block built along the cols first, the rows (limit collisions)
function update_hist_gpu!(
    h::CuArray{T,3},
    δ::CuMatrix{T},
    X_bin::CuMatrix{UInt8},
    𝑖::CuVector{S},
    𝑗::CuVector{S},
    K; MAX_THREADS=128) where {T,S}

    fill!(h, 0.0)
    thread_i = min(MAX_THREADS, length(𝑖))
    # thread_j = 1
    threads = (thread_i,)
    blocks = (1, length(𝑗), 2 * K + 1)

    @cuda blocks = blocks threads = threads shmem = sizeof(T) * size(h, 2) hist_kernel!(h, δ, X_bin, 𝑖, 𝑗)
    return
end

# update the vector of length 𝑖 pointing to associated node id
function update_set_kernel!(mask, set, best, x_bin)
    it = threadIdx().x
    ibd = blockDim().x
    ibi = blockIdx().x
    i = it + ibd * (ibi - 1)
    @inbounds if i <= length(set)
        @inbounds mask[i] = x_bin[set[i]] <= best
    end
    return nothing
end

function update_set_gpu(set, best, x_bin; MAX_THREADS=1024)
    mask = CUDA.zeros(Bool, length(set))
    thread_i = min(MAX_THREADS, length(set))
    threads = (thread_i,)
    blocks = (length(set) ÷ thread_i + 1,)
    @cuda blocks = blocks threads = threads update_set_kernel!(mask, set, best, x_bin)
    left, right = set[mask], set[.!mask]
    return left, right
end

# operate on hist_gpu
"""
find_split_gpu!
    Find best split over gpu histograms
"""

function find_split_gpu!(hist::AbstractArray{T,3}, edges::Vector{Vector{T}}, 𝑗::AbstractVector{S}, params::EvoTypes) where {T,S}

    hist_cum_L = cumsum(hist, dims=2)
    # hist_cum_R = sum(hist, dims=2) .- hist_cum_L
    hist_cum_R = hist_cum_L[:,end:end,:] .- hist_cum_L

    gains_L = get_hist_gains_gpu(hist_cum_L[:,1:(end - 1),:], 𝑗, params.λ)
    gains_R = get_hist_gains_gpu(hist_cum_R[:,1:(end - 1),:], 𝑗, params.λ)
    gains = gains_L + gains_R

    best = findmax(gains)
    gain, bin, feat = best[1], best[2][1], UInt32(best[2][2])
    cond = edges[feat][bin]
    # gainL, gainR = gains_L[bin, feat], gains_R[bin, feat]
    gainL, gainR = Array(gains_L)[bin, feat], Array(gains_R)[bin, feat]

    # ∑L = hist_cum_L[:, bin, feat]
    # ∑R = hist_cum_R[:, bin, feat]
    ∑L = Array(hist_cum_L[:, bin, feat])
    ∑R = Array(hist_cum_R[:, bin, feat])

    return (gain = gain, bin = bin, feat = feat, cond = cond,
        gainL = gainL, gainR = gainR,
        ∑L = ∑L, ∑R = ∑R)
end

function hist_gains_gpu!(gains::CuDeviceMatrix{T}, h::CuDeviceArray{T,3}, 𝑗::CuDeviceVector{S}, λ::T) where {T,S}

    i = threadIdx().x
    j = 𝑗[blockIdx().x]
    K = (size(h, 1) - 1) ÷ 2

    @inbounds 𝑤 = h[2 * K + 1, i, j]
    if 𝑤 > 1e-5
        @inbounds for k in 1:K
            if k == 1
                gains[i, j] = (h[k, i, j]^2 / (h[2 * K + k - 1, i, j] + λ * 𝑤)) / 2
            else
                gains[i, j] += (h[k, i, j]^2 / (h[2 * K + k - 1, i, j] + λ * 𝑤)) / 2
            end
        end
    end

    return nothing
end

function get_hist_gains_gpu(h::CuArray{T,3}, 𝑗::CuVector{S}, λ::T; MAX_THREADS=1024) where {T,S}

    gains = CUDA.fill(T(-Inf), size(h, 2) - 1, size(h, 3))

    thread_i = min(size(gains, 1), MAX_THREADS)
    threads = thread_i
    blocks = length(𝑗)

    @cuda blocks = blocks threads = threads hist_gains_gpu!(gains, h, 𝑗, λ)
    return gains
end

ext: jl


text:
"""
    VerticallyStretchedCartesianGrid{FT, TX, TY, TZ, R, A} <: AbstractGrid{FT, TX, TY, TZ}

A Cartesian grid with with constant horizontal grid spacings `Δx` and `Δy`, and
non-uniform or stretched vertical grid spacing `Δz` between cell centers and cell faces,
topology `{TX, TY, TZ}`, and coordinate ranges of type `R` (where a range can be used) and
`A` (where an array is needed).
"""
struct VerticallyStretchedCartesianGrid{FT, TX, TY, TZ, R, A} <: AbstractGrid{FT, TX, TY, TZ}

    # Number of grid points in (x,y,z).
     Nx :: Int
     Ny :: Int
     Nz :: Int

    # Halo size in (x,y,z).
     Hx :: Int
     Hy :: Int
     Hz :: Int

    # Domain size [m].
     Lx :: FT
     Ly :: FT
     Lz :: FT

    # Grid spacing [m].
     Δx :: FT
     Δy :: FT
    ΔzF :: A
    ΔzC :: A

    # Range of coordinates at the centers of the cells.
     xC :: R
     yC :: R
     zC :: A

    # Range of grid coordinates at the faces of the cells.
    # Note: there are Nx+1 faces in the x-dimension, Ny+1 in the y, and Nz+1 in the z.
     xF :: R
     yF :: R
     zF :: A
end

function VerticallyStretchedCartesianGrid(FT=Float64, arch=CPU();
                                              size, x, y, zF,
                                              halo = (1, 1, 1),
                                          topology = (Periodic, Periodic, Bounded))

    TX, TY, TZ = validate_topology(topology)
    Lx, Ly, x, y = validate_vertically_stretched_grid_size_and_xy(FT, size, halo, x, y)

    Nx, Ny, Nz = size
    Hx, Hy, Hz = halo

    # Initialize vertically-stretched arrays on CPU
    Lz, zF, zC, ΔzF, ΔzC = generate_stretched_vertical_grid(FT, topology[3], Nz, Hz, zF)

    # Convert to appropriate array type for arch
     zF = convert(array_type(arch), zF)
     zC = convert(array_type(arch), zC)
    ΔzF = convert(array_type(arch), ΔzF)
    ΔzC = convert(array_type(arch), ΔzC)

    # Construct uniform horizontal grid
    Lh, Nh, Hh, X₁ = (Lx, Ly), size[1:2], halo[1:2], (x[1], y[1])
    Δx, Δy = Δh = Lh ./ Nh

    # Face-node limits in x, y, z
    xF₋, yF₋ = XF₋ = @. X₁ - Hh * Δh
    xF₊, yF₊ = XF₊ = @. XF₋ + total_extent(topology[1:2], Hh, Δh, Lh)

    # Cell-node limits in x, y, z
    xC₋, yC₋ = XC₋ = @. XF₋ + Δh / 2
    xC₊, yC₊ = XC₊ = @. XC₋ + Lh + Δh * (2Hh - 1)

    # Total length of Cell and Face quantities
    TFx, TFy, TFz = total_length.(Face, topology, size, halo)
    TCx, TCy, TCz = total_length.(Cell, topology, size, halo)

    # Include halo points in coordinate arrays
    xF = range(xF₋, xF₊; length = TFx)
    yF = range(yF₋, yF₊; length = TFy)

    xC = range(xC₋, xC₊; length = TCx)
    yC = range(yC₋, yC₊; length = TCy)

    # Offset.
     xC = OffsetArray(xC,  -Hx)
     yC = OffsetArray(yC,  -Hy)
     zC = OffsetArray(zC,  -Hz)
    ΔzC = OffsetArray(ΔzC, 1 - Hz)

     xF = OffsetArray(xF,  -Hx)
     yF = OffsetArray(yF,  -Hy)
     zF = OffsetArray(zF,  -Hz)
    ΔzF = OffsetArray(ΔzF, -Hz)

    return VerticallyStretchedCartesianGrid{FT, TX, TY, TZ, typeof(xF), typeof(zF)}(
        Nx, Ny, Nz, Hx, Hy, Hz, Lx, Ly, Lz, Δx, Δy, ΔzF, ΔzC, xC, yC, zC, xF, yF, zF)
end

#####
##### Vertically stretched grid utilities
#####

get_z_face(z::Function, k) = z(k)
get_z_face(z::AbstractVector, k) = z[k]

lower_exterior_ΔzF(z_topo,          zFi, Hz) = [zFi[end - Hz + k] - zFi[end - Hz + k - 1] for k = 1:Hz]
lower_exterior_ΔzF(::Type{Bounded}, zFi, Hz) = [zFi[2]  - zFi[1] for k = 1:Hz]

upper_exterior_ΔzF(z_topo,          zFi, Hz) = [zFi[k + 1] - zFi[k] for k = 1:Hz]
upper_exterior_ΔzF(::Type{Bounded}, zFi, Hz) = [zFi[end]   - zFi[end - 1] for k = 1:Hz]

function generate_stretched_vertical_grid(FT, z_topo, Nz, Hz, zF_generator)

    # Ensure correct type for zF and derived quantities
    interior_zF = zeros(FT, Nz+1)

    for k = 1:Nz+1
        interior_zF[k] = get_z_face(zF_generator, k)
    end

    Lz = interior_zF[Nz+1] - interior_zF[1]

    # Build halo regions
    ΔzF₋ = lower_exterior_ΔzF(z_topo, interior_zF, Hz)
    ΔzF₊ = upper_exterior_ΔzF(z_topo, interior_zF, Hz)

    z¹, zᴺ⁺¹ = interior_zF[1], interior_zF[Nz+1]

    zF₋ = [z¹   - sum(ΔzF₋[k:Hz]) for k = 1:Hz] # locations of faces in lower halo
    zF₊ = [zᴺ⁺¹ + ΔzF₊[k]         for k = 1:Hz] # locations of faces in width of top halo region

    zF = vcat(zF₋, interior_zF, zF₊)

    # Build cell centers, cell center spacings, and cell interface spacings
    TCz = total_length(Cell, z_topo, Nz, Hz)
     zC = [ (zF[k + 1] + zF[k]) / 2 for k = 1:TCz ]
    ΔzC = [  zC[k] - zC[k - 1]      for k = 2:TCz ]

    # Trim face locations for periodic domains
    TFz = total_length(Face, z_topo, Nz, Hz)
    zF = zF[1:TFz]

    ΔzF = [zF[k + 1] - zF[k] for k = 1:TFz-1]

    return Lz, zF, zC, ΔzF, ΔzC
end


text:
module test_h3_lib_regions

using Test
using H3.Lib
using H3.API

Lib.polyfill
Lib.maxPolyfillSize
Lib.h3SetToLinkedGeo
Lib.destroyLinkedPolygon

using .Lib: H3Index, GeoCoord, Geofence, GeoPolygon

# https://github.com/uber/h3/blob/v3.7.2/src/apps/testapps/testPolyfill.c

sfVerts = [GeoCoord(0.659966917655, -2.1364398519396),
           GeoCoord(0.6595011102219, -2.1359434279405),
           GeoCoord(0.6583348114025, -2.1354884206045),
           GeoCoord(0.6581220034068, -2.1382437718946),
           GeoCoord(0.6594479998527, -2.1384597563896),
           GeoCoord(0.6599990002976, -2.1376771158464)]
sfGeofence = Geofence(6, pointer(sfVerts))
sfGeoPolygon = GeoPolygon(sfGeofence, 0, C_NULL)
@test Lib.maxPolyfillSize(Ref(sfGeoPolygon), 9) == 5613

holeVerts = [GeoCoord(0.6595072188743, -2.1371053983433),
             GeoCoord(0.6591482046471, -2.1373141048153),
             GeoCoord(0.6592295020837, -2.1365222838402)]
holeGeofence = Geofence(3, pointer(holeVerts))
holeGeoPolygon = GeoPolygon(sfGeofence, 1, pointer_from_objref(Ref(holeGeofence)))
@test Lib.maxPolyfillSize(Ref(holeGeoPolygon), 9) == 5613

emptyVerts = [GeoCoord(0.659966917655, -2.1364398519394),
              GeoCoord(0.659966917655, -2.1364398519395),
              GeoCoord(0.659966917655, -2.1364398519396)]
emptyGeofence = Geofence(3, pointer(emptyVerts))
emptyGeoPolygon = GeoPolygon(emptyGeofence, 0, C_NULL)
@test Lib.maxPolyfillSize(Ref(emptyGeoPolygon), 9) == 15

# https://h3geo.org/docs/api/regions

to_geocoord((lat, lon)) = GeoCoord(deg2rad(lat), deg2rad(lon))

sfVerts = map(to_geocoord, [
    [37.813318999983238, -122.4089866999972145],
    [37.7198061999978478, -122.3544736999993603],
    [37.8151571999998453, -122.4798767000009008]
])
res = 7

sfGeofence = Geofence(length(sfVerts), pointer(sfVerts))
sfGeoPolygon = GeoPolygon(sfGeofence, 0, C_NULL)
numHexagons = Lib.maxPolyfillSize(Ref(sfGeoPolygon), res)
@test numHexagons == 100

p_hexagons = Libc.calloc(numHexagons, sizeof(H3Index))
Lib.polyfill(Ref(sfGeoPolygon), res, p_hexagons)
p = Base.unsafe_convert(Ptr{H3Index}, p_hexagons)
w = unsafe_wrap(Vector{H3Index}, p, numHexagons)
@test filter(!iszero, w) == [0x087283082bffffff, 0x0872830870ffffff, 0x0872830820ffffff, 0x087283082effffff, 0x0872830828ffffff, 0x087283082affffff, 0x0872830876ffffff]
Libc.free(p_hexagons)

end # module test_h3_lib_regions


text:
"""

    AbstractMenu

The supertype for all Menu types.
See AbstractMenu.jl for descriptions of functions mentioned in this
doc string.

# Functions

The following functions can be called on all <:AbstractMenu types.
Details can be found in

## Exported

  - `request(m::AbstractMenu)`
  - `request(msg::AbstractString, m::AbstractMenu)`

## Hidden

  - `printMenu(m::AbstractMenu, cursor::Int; init::Bool=false)`

# Subtypes

All subtypes must contain the feilds `pagesize::Int` and
`pageoffset::Int`. They must also implement the following functions.

## Necessary Functions

These functions must be implemented for all subtypes of AbstractMenu.

  - `pick(m::AbstractMenu, cursor::Int)`
  - `cancel(m::AbstractMenu)`
  - `options(m::AbstractMenu)`
  - `writeLine(buf::IOBuffer, m::AbstractMenu, idx::Int, cur::Bool)`

## Optional Functions

These functions do not need to be implemented for all AbstractMenu
subtypes.

  - `header(m::AbstractMenu)`
  - `keypress(m::AbstractMenu, i::UInt32)`

"""
abstract type AbstractMenu end

# NECESSARY FUNCTIONS
# These functions must be implemented for all subtypes of AbstractMenu
######################################################################

# This function must be implemented for all menu types. It defines what
#   happens when a user presses the Enter key while the menu is open.
# If this function returns true, `request()` will exit.
pick(m::AbstractMenu, cursor::Int) = error("unimplemented")

# This function must be implemented for all menu types. It defines what
#   happends when a user cancels ('q' or ctrl-c) a menu. `request()` will
#   always exit after calling this function.
cancel(m::AbstractMenu) = error("unimplemented")

# This function must be implemented for all menu types. It should return
#   a list of strings to be displayed as options in the current page.
options(m::AbstractMenu) = error("unimplemented")

# This function must be implemented for all menu types. It should write
#   the option at index `idx` to the buffer. If cursor is `true` it
#   should also display the cursor
function writeLine(buf::IOBuffer, m::AbstractMenu, idx::Int, cur::Bool)
    error("unimplemented")
end

# OPTIONAL FUNCTIONS
# These functions do not need to be implemented for all Menu types
##################################################################

# If `header()` is defined for a specific menu type, display the header
#  above the menu when it is rendered to the screen.
header(m::AbstractMenu) = ""

# If `keypress()` is defined for a specific menu type, send any
#   non-standard keypres event to this function. If the function returns
#   true, `request()` will exit.
keypress(m::AbstractMenu, i::UInt32) = false

"""

    request(m::AbstractMenu)

Display the menu and enter interactive mode. Returns `m.selected` which
varies based on menu type.
"""
request(m::AbstractMenu) = request(terminal, m)

function request(term::Base.Terminals.TTYTerminal, m::AbstractMenu)
    cursor = 1

    menu_header = header(m)
    if !CONFIG[:supress_output] && menu_header != ""
        println(term.out_stream, menu_header)
    end

    printMenu(term.out_stream, m, cursor, init=true)

    raw_mode_enabled = enableRawMode(term)
    raw_mode_enabled && print(term.out_stream, "\x1b[?25l") # hide the cursor
    try
        while true
            c = readKey(term.in_stream)

            if c == Int(ARROW_UP)

                if cursor > 1
                    # move selection up
                    cursor -= 1
                    # scroll the page
                    if cursor < (2+m.pageoffset) && m.pageoffset > 0
                        m.pageoffset -= 1
                    end
                elseif CONFIG[:scroll_wrap]
                    # wrap to bottom
                    cursor = length(options(m))
                    m.pageoffset = length(options(m)) - m.pagesize
                end

            elseif c == Int(ARROW_DOWN)

                if cursor < length(options(m))
                    # move selection up
                    cursor += 1
                    # scroll page
                    if cursor >= m.pagesize + m.pageoffset && m.pagesize + m.pageoffset < length(options(m))
                        m.pageoffset += 1
                    end
                elseif CONFIG[:scroll_wrap]
                    # wrap to top
                    cursor = 1
                    m.pageoffset = 0
                end

            elseif c == Int(PAGE_UP)
                # If we're at the bottom, move the page 1 less to move the cursor up from
                # the bottom entry, since we try to avoid putting the cursor at bounds.
                m.pageoffset -= m.pagesize - (cursor == length(options(m)) ? 1 : 0)
                m.pageoffset = max(m.pageoffset, 0)
                cursor -= m.pagesize
                cursor = max(cursor, 1)

            elseif c == Int(PAGE_DOWN)
                m.pageoffset += m.pagesize - (cursor == 1 ? 1 : 0)
                m.pageoffset = min(m.pageoffset, length(options(m)) - m.pagesize)
                cursor += m.pagesize
                cursor = min(cursor, length(options(m)))

            elseif c == Int(HOME_KEY)
                cursor = 1
                m.pageoffset = 0

            elseif c == Int(END_KEY)
                cursor = length(options(m))
                m.pageoffset = length(options(m)) - m.pagesize

            elseif c == 13 # <enter>
                # will break if pick returns true
                pick(m, cursor) && break
            elseif c == UInt32('q')
                cancel(m)
                break
            elseif c == 3 # ctrl-c
                cancel(m)
                if CONFIG[:ctrl_c_interrupt]
                    throw(InterruptException())
                else
                    break
                end
            else
                # will break if keypress returns true
                keypress(m, c) && break
            end

            printMenu(term.out_stream, m, cursor)
        end
    finally
        # always disable raw mode even even if there is an
        #  exception in the above loop
        if raw_mode_enabled
            print(term.out_stream, "\x1b[?25h") # unhide cursor
            disableRawMode(term)
        end
    end
    println(term.out_stream)

    return m.selected
end

"""

    request([term,] msg::AbstractString, m::AbstractMenu)

Shorthand for `println(msg); request(m)`.
"""
request(msg::AbstractString, m::AbstractMenu) =
    request(terminal, msg, m)

function request(term::Base.Terminals.TTYTerminal,
                 msg::AbstractString, m::AbstractMenu)
    println(term.out_stream, msg)
    request(term, m)
end

# The generic printMenu function is used for displaying the state of a
#   menu to the screen. Menus must implement `writeLine` and `options`
#   and have fields `pagesize::Int` and `pageoffset::Int` as part of
#   their type definition
function printMenu(out, m::AbstractMenu, cursor::Int; init::Bool=false)
    CONFIG[:supress_output] && return

    buf = IOBuffer()

    # Move the cursor to the beginning of where it should print
    # Don't do this on the initial print
    lines = m.pagesize-1
    if init
        m.pageoffset = 0
    else
        print(buf, "\x1b[999D\x1b[$(lines)A")
    end

    for i in (m.pageoffset+1):(m.pageoffset + m.pagesize)
        print(buf, "\x1b[2K")

        if i == m.pageoffset+1 && m.pageoffset > 0
            # first line && scrolled past first entry
            print(buf, CONFIG[:up_arrow])
        elseif i == m.pagesize+m.pageoffset && i != length(options(m))
            # last line && not last option
            print(buf, CONFIG[:down_arrow])
        else
            # non special line
            print(buf, " ")
        end

        writeLine(buf, m, i, i == cursor)

        # dont print an \r\n on the last line
        i != (m.pagesize+m.pageoffset) && print(buf, "\r\n")
    end

    print(out, String(take!(buf)))
end


text:
# Translating a SQL node graph to a SQL statement.

function render(n; dialect = :default)
    actx = AnnotateContext()
    n′ = annotate(convert(SQLNode, n), actx)
    @debug "FunSQL.annotate\n" * sprint(pprint, n′) _group = Symbol("FunSQL.annotate")
    resolve!(actx)
    @debug "FunSQL.resolve!\n" * sprint(pprint, n′) _group = Symbol("FunSQL.resolve!")
    link!(actx)
    @debug "FunSQL.link!\n" * sprint(pprint, n′) _group = Symbol("FunSQL.link!")
    tctx = TranslateContext(dialect, actx.path_map)
    c = translate(n′, tctx)
    @debug "FunSQL.translate\n" * sprint(pprint, c) _group = Symbol("FunSQL.translate")
    sql = render(c, dialect = dialect)
    @debug "FunSQL.render\n" * sql _group = Symbol("FunSQL.render")
    sql
end

struct TranslateContext
    dialect::SQLDialect
    path_map::PathMap
    aliases::Dict{Symbol, Int}
    vars::Dict{Symbol, SQLClause}
    subs::Dict{SQLNode, SQLClause}

    TranslateContext(dialect, path_map::PathMap) =
        new(dialect, path_map, Dict{Symbol, Int}(), Dict{Symbol, SQLClause}(), Dict{SQLNode, SQLClause}())

    function TranslateContext(ctx::TranslateContext; vars = nothing, subs = nothing)
        new(ctx.dialect, ctx.path_map, ctx.aliases, something(vars, ctx.vars), something(subs, ctx.subs))
    end
end

allocate_alias(ctx::TranslateContext, n::SQLNode) =
    allocate_alias(ctx, (n[]::BoxNode).type.name)

function allocate_alias(ctx::TranslateContext, alias::Symbol)
    n = get(ctx.aliases, alias, 0) + 1
    ctx.aliases[alias] = n
    Symbol(alias, '_', n)
end

# Translating scalar nodes.

function translate(n, ctx::TranslateContext, subs::Dict{SQLNode, SQLClause})
    ctx′ = TranslateContext(ctx, subs = subs)
    translate(n, ctx′)
end

function translate(n::SQLNode, ctx)
    c = get(ctx.subs, n, nothing)
    if c === nothing
        c = convert(SQLClause, translate(n[], ctx))
    end
    c
end

translate(ns::Vector{SQLNode}, ctx) =
    SQLClause[translate(n, ctx) for n in ns]

translate(::Nothing, ctx) =
    nothing

translate(n::AggregateNode, ctx) =
    translate(Val(n.name), n, ctx)

translate(@nospecialize(name::Val{N}), n::AggregateNode, ctx) where {N} =
    translate_default(n, ctx)

function translate_default(n::AggregateNode, ctx)
    args = translate(n.args, ctx)
    filter = translate(n.filter, ctx)
    AGG(uppercase(string(n.name)), distinct = n.distinct, args = args, filter = filter)
end

function translate(::Val{:count}, n::AggregateNode, ctx)
    args = !isempty(n.args) ? translate(n.args, ctx) : [OP("*")]
    filter = translate(n.filter, ctx)
    AGG(:COUNT, distinct = n.distinct, args = args, filter = filter)
end

translate(n::Union{AsNode, HighlightNode}, ctx) =
    translate(n.over, ctx)

function translate(n::ExtendedBindNode, ctx)
    vars′ = copy(ctx.vars)
    for (name, i) in n.label_map
        vars′[name] = translate(n.list[i], ctx)
    end
    ctx′ = TranslateContext(ctx, vars = vars′)
    translate(n.over, ctx′)
end

function translate(n::BoxNode, ctx)
    res = assemble(n, ctx)
    complete(res)
end

translate(n::FunctionNode, ctx) =
    translate(Val(n.name), n, ctx)

translate(@nospecialize(name::Val{N}), n::FunctionNode, ctx) where {N} =
    translate_default(n, ctx)

function translate_default(n::FunctionNode, ctx)
    args = translate(n.args, ctx)
    if Base.isidentifier(n.name)
        FUN(uppercase(string(n.name)), args = args)
    else
        OP(n.name, args = args)
    end
end

for (name, op) in (:not => :NOT,
                   :like => :LIKE,
                   :exists => :EXISTS,
                   :(==) => Symbol("="),
                   :(!=) => Symbol("<>"))
    @eval begin
        translate(::Val{$(QuoteNode(name))}, n::FunctionNode, ctx) =
            OP($(QuoteNode(op)),
               args = SQLClause[translate(arg, ctx) for arg in n.args])
    end
end

for (name, op, default) in ((:and, :AND, true), (:or, :OR, false))
    @eval begin
        function translate(::Val{$(QuoteNode(name))}, n::FunctionNode, ctx)
            args = translate(n.args, ctx)
            if isempty(args)
                LIT($default)
            elseif length(args) == 1
                args[1]
            elseif length(args) == 2 && (@dissect args[1] LIT(val = val)) && val == $default
                args[2]
            elseif (@dissect args[1] OP(name = name, args = args′)) && name === $(QuoteNode(op))
                OP($(QuoteNode(op)), args = SQLClause[args′..., args[2:end]...])
            else
                OP($(QuoteNode(op)), args = args)
            end
        end
    end
end

for (name, op, default) in (("in", "IN", false), ("not in", "NOT IN", true))
    @eval begin
        function translate(::Val{Symbol($name)}, n::FunctionNode, ctx)
            if length(n.args) <= 1
                LIT($default)
            else
                args = translate(n.args, ctx)
                if length(args) == 2 && @dissect args[2] (SELECT() || UNION())
                    OP($op, args = args)
                else
                    OP($op, args[1], FUN("", args = args[2:end]))
                end
            end
        end
    end
end

translate(::Val{Symbol("is null")}, n::FunctionNode, ctx) =
    OP(:IS, SQLClause[translate(arg, ctx) for arg in n.args]..., missing)

translate(::Val{Symbol("is not null")}, n::FunctionNode, ctx) =
    OP("IS NOT", SQLClause[translate(arg, ctx) for arg in n.args]..., missing)

translate(::Val{:case}, n::FunctionNode, ctx) =
    CASE(args = SQLClause[translate(arg, ctx) for arg in n.args])

for (name, op) in (("between", "BETWEEN"), ("not between", "NOT BETWEEN"))
    @eval begin
        function translate(::Val{Symbol($name)}, n::FunctionNode, ctx)
            if length(n.args) == 3
                args = SQLClause[translate(arg, ctx) for arg in n.args]
                OP($op, args[1], args[2], args[3] |> KW(:AND))
            else
                translate_default(n, ctx)
            end
        end
    end
end

for (name, op) in (("current_date", "CURRENT_DATE"),
                   ("current_timestamp", "CURRENT_TIMESTAMP"))
    @eval begin
        function translate(::Val{Symbol($name)}, n::FunctionNode, ctx)
            if isempty(n.args)
                OP($op)
            else
                translate_default(n, ctx)
            end
        end
    end
end

translate(n::LiteralNode, ctx) =
    LIT(n.val)

translate(n::SortNode, ctx) =
    SORT(over = translate(n.over, ctx), value = n.value, nulls = n.nulls)

function translate(n::VariableNode, ctx)
    c = get(ctx.vars, n.name, nothing)
    if c === nothing
        c = VAR(n.name)
    end
    c
end

# Translating subquery nodes.

struct Assemblage
    clause::Union{SQLClause, Nothing}
    cols::OrderedDict{Symbol, SQLClause}
    repl::Dict{SQLNode, Symbol}

    Assemblage(clause; cols = OrderedDict{Symbol, SQLClause}(), repl = Dict{SQLNode, Symbol}()) =
        new(clause, cols, repl)
end

function complete(cols::OrderedDict{Symbol, SQLClause})
    list = SQLClause[]
    for (name, c) in cols
        if !((@dissect c ID(name = id_name)) && id_name == name)
            c = AS(over = c, name = name)
        end
        push!(list, c)
    end
    if isempty(list)
        push!(list, LIT(missing))
    end
    list
end

function complete(res::Assemblage)
    clause = res.clause
    if !(@dissect clause SELECT() || UNION())
        list = complete(res.cols)
        clause = SELECT(over = clause, list = list)
    else
        @assert isempty(res.cols)
    end
    @assert clause !== nothing
    clause
end

function make_subs(base::Assemblage, ::Nothing)
    subs = Dict{SQLNode, SQLClause}()
    for (ref, name) in base.repl
        subs[ref] = base.cols[name]
    end
    subs
end

function make_subs(base::Assemblage, base_as::Symbol)
    subs = Dict{SQLNode, SQLClause}()
    base_cache = Dict{Symbol, SQLClause}()
    for (ref, name) in base.repl
        subs[ref] = get(base_cache, name) do
            ID(over = base_as, name = name)
        end
    end
    subs
end

function make_repl(refs::Vector{SQLNode})::Dict{SQLNode, Symbol}
    repl = Dict{SQLNode, Symbol}()
    dups = Dict{Symbol, Int}()
    for ref in refs
        name′ = name = label(ref)
        k = get(dups, name, 0) + 1
        if k > 1
            name′ = Symbol(name, '_', k)
            while name′ in keys(dups)
                k += 1
                name′ = Symbol(name, '_', k)
            end
            dups[name] = k
        end
        repl[ref] = name′
        dups[name′] = 1
    end
    repl
end

function make_repl(trns::Vector{Pair{SQLNode, SQLClause}})::Tuple{Dict{SQLNode, Symbol}, OrderedDict{Symbol, SQLClause}}
    repl = Dict{SQLNode, Symbol}()
    cols = OrderedDict{Symbol, SQLClause}()
    dups = Dict{Symbol, Int}()
    renames = Dict{Tuple{Symbol, SQLClause}, Symbol}()
    for (ref, c) in trns
        name′ = name = label(ref)
        k = get(dups, name, 0) + 1
        if k > 1
            name′ = get(renames, (name, c), nothing)
            if name′ !== nothing
                repl[ref] = name′
                continue
            end
            name′ = Symbol(name, '_', k)
            while name′ in keys(dups)
                k += 1
                name′ = Symbol(name, '_', k)
            end
            dups[name] = k
        end
        cols[name′] = c
        dups[name′] = 1
        renames[name, c] = name′
        repl[ref] = name′
    end
    (repl, cols)
end

assemble(n::SQLNode, ctx) =
    assemble(n[], ctx)

function assemble(n::BoxNode, ctx)
    refs′ = SQLNode[]
    for ref in n.refs
        if (@dissect ref over |> HandleBoundNode(handle = handle)) && handle == n.handle
            push!(refs′, over)
        else
            push!(refs′, ref)
        end
    end
    res = assemble(n.over, refs′, ctx)
    repl′ = Dict{SQLNode, Symbol}()
    for ref in n.refs
        if (@dissect ref over |> HandleBoundNode(handle = handle)) && handle == n.handle
            repl′[ref] = res.repl[over]
        else
            repl′[ref] = res.repl[ref]
        end
    end
    Assemblage(res.clause, cols = res.cols, repl = repl′)
end

assemble(n::SQLNode, refs, ctx) =
    assemble(n[], refs, ctx)

function assemble(::Nothing, refs, ctx)
    @assert isempty(refs)
    Assemblage(nothing)
end

function assemble(n::AppendNode, refs, ctx)
    base = assemble(n.over, ctx)
    results = [n.over => base]
    for l in n.list
        res = assemble(l, ctx)
        push!(results, l => res)
    end
    dups = Dict{SQLNode, SQLNode}()
    seen = Dict{Symbol, SQLNode}()
    for ref in refs
        name = base.repl[ref]
        if name in keys(seen)
            other_ref = seen[name]
            if all(res.repl[ref] === res.repl[other_ref] for (l, res) in results)
                dups[ref] = seen[name]
            end
        else
            seen[name] = ref
        end
    end
    urefs = SQLNode[ref for ref in refs if !(ref in keys(dups))]
    repl = make_repl(urefs)
    for (ref, uref) in dups
        repl[ref] = repl[uref]
    end
    cs = SQLClause[]
    for (l, res) in results
        if !(@dissect res.clause (SELECT() || UNION()))
            as = nothing
            tail = res.clause
        else
            as = allocate_alias(ctx, l)
            tail = FROM(AS(over = complete(res), name = as))
        end
        subs = make_subs(res, as)
        cols = OrderedDict{Symbol, SQLClause}()
        for ref in refs
            !(ref in keys(dups)) || continue
            name = repl[ref]
            cols[name] = subs[ref]
        end
        c = SELECT(over = tail, list = complete(cols))
        push!(cs, c)
    end
    c = UNION(over = cs[1], all = true, list = cs[2:end])
    Assemblage(c, repl = repl)
end

function assemble(n::AsNode, refs, ctx)
    res = assemble(n.over, ctx)
    repl′ = Dict{SQLNode, Symbol}()
    for ref in refs
        if @dissect ref over |> NameBound()
            @assert over !== nothing
            repl′[ref] = res.repl[over]
        else
            repl′[ref] = res.repl[ref]
        end
    end
    Assemblage(res.clause, cols = res.cols, repl = repl′)
end

function assemble(n::ExtendedBindNode, refs, ctx)
    vars′ = copy(ctx.vars)
    for (name, i) in n.label_map
        vars′[name] = translate(n.list[i], ctx)
    end
    ctx′ = TranslateContext(ctx, vars = vars′)
    assemble(n.over, ctx′)
end

function assemble(n::DefineNode, refs, ctx)
    base = assemble(n.over, ctx)
    if !any(ref -> (@dissect ref Get(name = name)) && name in keys(n.label_map), refs)
        return base
    end
    if !(@dissect base.clause (SELECT() || UNION()))
        base_as = nothing
        c = base.clause
    else
        base_as = allocate_alias(ctx, n.over)
        c = FROM(AS(over = complete(base), name = base_as))
    end
    subs = make_subs(base, base_as)
    repl = Dict{SQLNode, Symbol}()
    trns = Pair{SQLNode, SQLClause}[]
    tr_cache = Dict{Symbol, SQLClause}()
    for ref in refs
        if (@dissect ref (nothing |> Get(name = name))) && name in keys(n.label_map)
            col = get!(tr_cache, name) do
                def = n.list[n.label_map[name]]
                translate(def, ctx, subs)
            end
            push!(trns, ref => col)
        else
            push!(trns, ref => subs[ref])
        end
    end
    repl, cols = make_repl(trns)
    Assemblage(c, cols = cols, repl = repl)
end

function assemble(n::ExtendedJoinNode, refs, ctx)
    left = assemble(n.over, ctx)
    if @dissect left.clause (tail := FROM() || JOIN())
        left_as = nothing
    else
        left_as = allocate_alias(ctx, n.over)
        tail = FROM(AS(over = complete(left), name = left_as))
    end
    lateral = !isempty(n.lateral)
    subs = make_subs(left, left_as)
    if lateral
        right = assemble(n.joinee, TranslateContext(ctx, subs = subs))
    else
        right = assemble(n.joinee, ctx)
    end
    if @dissect right.clause ((joinee := nothing |> ID() |> AS(name = right_as)) |> FROM())
        for (ref, name) in right.repl
            subs[ref] = right.cols[name]
        end
    else
        right_as = allocate_alias(ctx, n.joinee)
        joinee = AS(over = complete(right), name = right_as)
        right_cache = Dict{Symbol, SQLClause}()
        for (ref, name) in right.repl
            subs[ref] = get(right_cache, name) do
                ID(over = right_as, name = name)
            end
        end
    end
    on = translate(n.on, ctx, subs)
    trns = Pair{SQLNode, SQLClause}[]
    for ref in refs
        push!(trns, ref => subs[ref])
    end
    repl, cols = make_repl(trns)
    c = JOIN(over = tail, joinee = joinee, on = on, left = n.left, right = n.right, lateral = lateral)
    Assemblage(c, cols = cols, repl = repl)
end

function assemble(n::FromNode, refs, ctx)
    output_columns = Set{Symbol}()
    for ref in refs
        match = @dissect ref (nothing |> Get(name = name))
        @assert match && name in n.table.column_set
        if !(name in output_columns)
            push!(output_columns, name)
        end
    end
    as = allocate_alias(ctx, n.table.name)
    cols = OrderedDict{Symbol, SQLClause}()
    for col in n.table.columns
        col in output_columns || continue
        cols[col] = ID(over = as, name = col)
    end
    tbl = ID(over = n.table.schema, name = n.table.name)
    c = FROM(AS(over = tbl, name = as))
    repl = Dict{SQLNode, Symbol}()
    for ref in refs
        if @dissect ref (nothing |> Get(name = name))
            repl[ref] = name
        end
    end
    Assemblage(c, cols = cols, repl = repl)
end

function assemble(n::GroupNode, refs, ctx)
    has_aggregates = any(ref -> (@dissect ref Agg()), refs)
    if isempty(n.by) && !has_aggregates
        return assemble(nothing, refs, ctx)
    end
    base = assemble(n.over, ctx)
    if @dissect base.clause (tail := nothing || FROM() || JOIN() || WHERE())
        base_as = nothing
    else
        base_as = allocate_alias(ctx, n.over)
        tail = FROM(AS(over = complete(base), name = base_as))
    end
    subs = make_subs(base, base_as)
    by = translate(n.by, ctx, subs)
    trns = Pair{SQLNode, SQLClause}[]
    for ref in refs
        if @dissect ref (nothing |> Get(name = name))
            @assert name in keys(n.label_map)
            ckey = by[n.label_map[name]]
            push!(trns, ref => ckey)
        elseif @dissect ref (nothing |> Agg(name = name))
            c = translate(ref, ctx, subs)
            push!(trns, ref => c)
        end
    end
    repl, cols = make_repl(trns)
    @assert !isempty(cols)
    if has_aggregates
        c = GROUP(over = tail, by = by)
        return Assemblage(c, cols = cols, repl = repl)
    else
        list = complete(cols)
        c = SELECT(over = tail, distinct = true, list = list)
        return Assemblage(c, repl = repl)
    end
end

assemble(n::HighlightNode, refs, ctx) =
    assemble(n.over, ctx)

function assemble(n::LimitNode, refs, ctx)
    base = assemble(n.over, ctx)
    if n.offset === nothing && n.limit === nothing
        return base
    end
    if @dissect base.clause (tail := nothing || FROM() || JOIN() || WHERE() || GROUP() || HAVING() || ORDER())
        base_as = nothing
    else
        base_as = allocate_alias(ctx, n.over)
        tail = FROM(AS(over = complete(base), name = base_as))
    end
    subs = make_subs(base, base_as)
    trns = Pair{SQLNode, SQLClause}[]
    for ref in refs
        push!(trns, ref => subs[ref])
    end
    repl, cols = make_repl(trns)
    c = LIMIT(over = tail, offset = n.offset, limit = n.limit)
    Assemblage(c, cols = cols, repl = repl)
end

function assemble(n::OrderNode, refs, ctx)
    base = assemble(n.over, ctx)
    if isempty(n.by)
        return base
    end
    if @dissect base.clause (tail := nothing || FROM() || JOIN() || WHERE() || GROUP() || HAVING())
        base_as = nothing
    else
        base_as = allocate_alias(ctx, n.over)
        tail = FROM(AS(over = complete(base), name = base_as))
    end
    subs = make_subs(base, base_as)
    by = translate(n.by, ctx, subs)
    trns = Pair{SQLNode, SQLClause}[]
    for ref in refs
        push!(trns, ref => subs[ref])
    end
    repl, cols = make_repl(trns)
    c = ORDER(over = tail, by = by)
    Assemblage(c, cols = cols, repl = repl)
end

function assemble(n::PartitionNode, refs, ctx)
    base = assemble(n.over, ctx)
    if !any(ref -> (@dissect ref Agg()), refs)
        return base
    end
    if @dissect base.clause (tail := nothing || FROM() || JOIN() || WHERE() || GROUP() || HAVING())
        base_as = nothing
    else
        base_as = allocate_alias(ctx, n.over)
        tail = FROM(AS(over = complete(base), name = base_as))
    end
    subs = make_subs(base, base_as)
    ctx′ = TranslateContext(ctx, subs = subs)
    by = translate(n.by, ctx′)
    order_by = translate(n.order_by, ctx′)
    partition = PARTITION(by = by, order_by = order_by, frame = n.frame)
    trns = Pair{SQLNode, SQLClause}[]
    for ref in refs
        if @dissect ref (nothing |> Agg(name = name))
            c = partition |> translate(ref, ctx′)
            push!(trns, ref => c)
        else
            push!(trns, ref => subs[ref])
        end
    end
    repl, cols = make_repl(trns)
    c = WINDOW(over = tail, list = [])
    Assemblage(c, cols = cols, repl = repl)
end

function assemble(n::SelectNode, refs, ctx)
    base = assemble(n.over, ctx)
    if !(@dissect base.clause (SELECT() || UNION()))
        base_as = nothing
        tail = base.clause
    else
        base_as = allocate_alias(ctx, n.over)
        tail = FROM(AS(over = complete(base), name = base_as))
    end
    subs = make_subs(base, base_as)
    cols = OrderedDict{Symbol, SQLClause}()
    for (name, i) in n.label_map
        col = n.list[i]
        cols[name] = translate(col, ctx, subs)
    end
    c = SELECT(over = tail, list = complete(cols))
    repl = Dict{SQLNode, Symbol}()
    for ref in refs
        ref_name = nothing
        @dissect ref (nothing |> Get(name = name))
        @assert name !== nothing
        repl[ref] = name
    end
    Assemblage(c, repl = repl)
end

function merge_conditions(c1, c2)
    if @dissect c1 OP(name = :AND, args = args1)
        if @dissect c2 OP(name = :AND, args = args2)
            return OP(:AND, args1..., args2...)
        else
            return OP(:AND, args1..., c2)
        end
    elseif @dissect c2 OP(name = :AND, args = args2)
        return OP(:AND, c1, args2...)
    else
        return OP(:AND, c1, c2)
    end
end

function assemble(n::WhereNode, refs, ctx)
    base = assemble(n.over, ctx)
    if (@dissect base.clause (nothing || FROM() || JOIN() || WHERE() || HAVING())) ||
       (@dissect base.clause GROUP(by = by)) && !isempty(by)
        subs = make_subs(base, nothing)
        condition = translate(n.condition, ctx, subs)
        if (@dissect condition LIT(val = val)) && val === true
            return base
        end
        if @dissect base.clause tail |> WHERE(condition = tail_condition)
            condition = merge_conditions(tail_condition, condition)
            c = WHERE(over = tail, condition = condition)
        elseif @dissect base.clause GROUP()
            c = HAVING(over = base.clause, condition = condition)
        elseif @dissect base.clause tail |> HAVING(condition = tail_condition)
            condition = merge_conditions(tail_condition, condition)
            c = HAVING(over = tail, condition = condition)
        else
            c = WHERE(over = base.clause, condition = condition)
        end
    else
        base_as = allocate_alias(ctx, n.over)
        tail = FROM(AS(over = complete(base), name = base_as))
        subs = make_subs(base, base_as)
        condition = translate(n.condition, ctx, subs)
        if (@dissect condition LIT(val = val)) && val === true
            return base
        end
        c = WHERE(over = tail, condition = condition)
    end
    trns = Pair{SQLNode, SQLClause}[]
    for ref in refs
        push!(trns, ref => subs[ref])
    end
    repl, cols = make_repl(trns)
    return Assemblage(c, cols = cols, repl = repl)
end


text:
#############################################################
# Fornberg algorithm

# This implements the Fornberg (1988) algorithm (https://doi.org/10.1090/S0025-5718-1988-0935077-0)
# to obtain Finite Difference weights over arbitrary points to arbitrary order

function calculate_weights(order::Int, x0::T, x::AbstractVector) where T<:Real
    #=
        order: The derivative order for which we need the coefficients
        x0   : The point in the array 'x' for which we need the coefficients
        x    : A dummy array with relative coordinates, eg. central differences
               need coordinates centred at 0 while those at boundaries need
               coordinates starting from 0 to the end point

        The approximation order of the stencil is automatically determined from
        the number of requested stencil points.
    =#
    N = length(x)
    @assert order < N "Not enough points for the requested order."
    M = order
    c1 = one(T)
    c4 = x[1] - x0
    C = zeros(T, N, M+1)
    C[1,1] = 1
    @inbounds for i in 1 : N-1
        i1 = i + 1
        mn = min(i, M)
        c2 = one(T)
        c5 = c4
        c4 = x[i1] - x0
        for j in 0 : i-1
            j1 = j + 1
            c3 = x[i1] - x[j1]
            c2 *= c3
            if j == i-1
                for s in mn : -1 : 1
                    s1 = s + 1
                    C[i1,s1] = c1*(s*C[i,s] - c5*C[i,s1]) / c2
                end
                C[i1,1] = -c1*c5*C[i,1] / c2
           end
            for s in mn : -1 : 1
                s1 = s + 1
                C[j1,s1] = (c4*C[j1,s1] - s*C[j1,s]) / c3
            end
            C[j1,1] = c4 * C[j1,1] / c3
        end
        c1 = c2
    end
    #=
        This is to fix the problem of numerical instability which occurs when the sum of the stencil_coefficients is not
        exactly 0.
        https://scicomp.stackexchange.com/questions/11249/numerical-derivative-and-finite-difference-coefficients-any-update-of-the-fornb
        Stack Overflow answer on this issue.
        http://epubs.siam.org/doi/pdf/10.1137/S0036144596322507 - Modified Fornberg Algorithm
    =#
    _C = C[:,end]
    _C[div(N,2)+1] -= sum(_C)
    return _C
end


text:
@testset "minimal_cancel_plus_test" begin
    @test cancelplus(interval(-Inf, -1.0), emptyinterval()) == entireinterval()
    @test cancelplus(interval(-1.0, Inf), emptyinterval()) == entireinterval()
    @test cancelplus(entireinterval(), emptyinterval()) == entireinterval()
    @test cancelplus(interval(-Inf, -1.0), interval(-5.0,1.0)) == entireinterval()
    @test cancelplus(interval(-1.0, Inf), interval(-5.0,1.0)) == entireinterval()
    @test cancelplus(entireinterval(), interval(-5.0,1.0)) == entireinterval()
    @test cancelplus(interval(-Inf, -1.0), entireinterval()) == entireinterval()
    @test cancelplus(interval(-1.0, Inf), entireinterval()) == entireinterval()
    @test cancelplus(emptyinterval(), interval(1.0, Inf)) == entireinterval()
    @test cancelplus(emptyinterval(), interval(-Inf,1.0)) == entireinterval()
    @test cancelplus(emptyinterval(), entireinterval()) == entireinterval()
    @test cancelplus(interval(-1.0,5.0), interval(1.0,Inf)) == entireinterval()
    @test cancelplus(interval(-1.0,5.0), interval(-Inf,1.0)) == entireinterval()
    @test cancelplus(interval(-1.0,5.0), entireinterval()) == entireinterval()
    @test cancelplus(entireinterval(), interval(1.0,Inf)) == entireinterval()
    @test cancelplus(entireinterval(), interval(-Inf,1.0)) == entireinterval()
    @test cancelplus(entireinterval(), entireinterval()) == entireinterval()
    @test cancelplus(interval(-5.0, -1.0), interval(1.0,5.1)) == entireinterval()
    @test cancelplus(interval(-5.0, -1.0), interval(0.9,5.0)) == entireinterval()
    @test cancelplus(interval(-5.0, -1.0), interval(0.9,5.1)) == entireinterval()
    @test cancelplus(interval(-10.0, 5.0), interval(-5.0,10.1)) == entireinterval()
    @test cancelplus(interval(-10.0, 5.0), interval(-5.1,10.0)) == entireinterval()
    @test cancelplus(interval(-10.0, 5.0), interval(-5.1,10.1)) == entireinterval()
    @test cancelplus(interval(1.0, 5.0), interval(-5.0,-0.9)) == entireinterval()
    @test cancelplus(interval(1.0, 5.0), interval(-5.1,-1.0)) == entireinterval()
    @test cancelplus(interval(1.0, 5.0), interval(-5.1,-0.9)) == entireinterval()
    @test cancelplus(interval(-10.0, -1.0), emptyinterval()) == entireinterval()
    @test cancelplus(interval(-10.0, 5.0), emptyinterval()) == entireinterval()
    @test cancelplus(interval(1.0, 5.0), emptyinterval()) == entireinterval()
    @test cancelplus(emptyinterval(), emptyinterval()) == emptyinterval()
    @test cancelplus(emptyinterval(), interval(1.0,10.0)) == emptyinterval()
    @test cancelplus(emptyinterval(), interval(-5.0,10.0)) == emptyinterval()
    @test cancelplus(emptyinterval(), interval(-5.0,-1.0)) == emptyinterval()
    @test cancelplus(interval(-5.1,-0.0), interval(0.0,5.0)) == interval(-0x1.999999999998P-4,0.0)
    @test cancelplus(interval(-5.1,-1.0), interval(1.0,5.0)) == interval(-0x1.999999999998P-4,0.0)
    @test cancelplus(interval(-5.0,-0.9), interval(1.0,5.0)) == interval(0.0, 0x1.9999999999998P-4)
    @test cancelplus(interval(-5.1,-0.9), interval(1.0,5.0)) == interval(-0x1.999999999998P-4,0x1.9999999999998P-4)
    @test cancelplus(interval(-5.0,-1.0), interval(1.0,5.0)) == interval(0.0,0.0)
    @test cancelplus(interval(-10.1, 5.0), interval(-5.0,10.0)) == interval(-0x1.999999999998P-4,0.0)
    @test cancelplus(interval(-10.0, 5.1), interval(-5.0,10.0)) == interval(0.0,0x1.999999999998P-4)
    @test cancelplus(interval(-10.1, 5.1), interval(-5.0,10.0)) == interval(-0x1.999999999998P-4,0x1.999999999998P-4)
    @test cancelplus(interval(-10.0, 5.0), interval(-5.0,10.0)) == interval(0.0,0.0)
    @test cancelplus(interval(0.9, 5.0), interval(-5.0,-1.0)) == interval(-0x1.9999999999998P-4,0.0)
    @test cancelplus(interval(1.0, 5.1), interval(-5.0,-1.0)) == interval(0.0,0x1.999999999998P-4)
    @test cancelplus(interval(0.0, 5.1), interval(-5.0,-0.0)) == interval(0.0,0x1.999999999998P-4)
    @test cancelplus(interval(0.9, 5.1), interval(-5.0,-1.0)) == interval(-0x1.9999999999998P-4,0x1.999999999998P-4)
    @test cancelplus(interval(1.0, 5.0), interval(-5.0,-1.0)) == interval(0.0,0.0)
    @test cancelplus(interval(0.0, 5.0), interval(-5.0,-0.0)) == interval(0.0,0.0)
    @test cancelplus(interval(0x1.FFFFFFFFFFFFP+0,0x1.FFFFFFFFFFFFP+0), interval(-0x1.999999999999AP-4,-0x1.999999999999AP-4)) == interval(0x1.E666666666656P+0,0x1.E666666666657P+0)
    @test cancelplus(interval(-0x1.999999999999AP-4,0x1.FFFFFFFFFFFFP+0), interval(-0x1.999999999999AP-4,0.01)) == interval(-0x1.70A3D70A3D70BP-4,0x1.E666666666657P+0)
    @test cancelplus(interval(0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFFp1023), interval(0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFFp1023)) == interval(0x1.FFFFFFFFFFFFFp1023,Inf)
    @test cancelplus(interval(-0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFFp1023), interval(-0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFFp1023)) == interval(0.0,0.0)
    @test cancelplus(interval(-0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFFp1023), interval(-0x1.FFFFFFFFFFFFEP+1023,0x1.FFFFFFFFFFFFFp1023)) == interval(0.0,0x1P+971)
    @test cancelplus(interval(-0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFFp1023), interval(-0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFEP+1023)) == interval(-0x1P+971,0.0)
    @test cancelplus(interval(-0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFEP+1023), interval(-0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFFp1023)) == entireinterval()
    @test cancelplus(interval(-0x1.FFFFFFFFFFFFEP+1023,0x1.FFFFFFFFFFFFFp1023), interval(-0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFFp1023)) == entireinterval()
    @test cancelplus(interval(-0x1P+0,0x1.FFFFFFFFFFFFFP-53), interval(-0x1P+0,0x1.FFFFFFFFFFFFEP-53)) == interval(-0x1.FFFFFFFFFFFFFP-1,-0x1.FFFFFFFFFFFFEP-1)
    @test cancelplus(interval(-0x1P+0,0x1.FFFFFFFFFFFFEP-53), interval(-0x1P+0,0x1.FFFFFFFFFFFFFP-53)) == entireinterval()
end

@testset "minimal_cancel_minus_test" begin
    @test cancelminus(interval(-Inf, -1.0), emptyinterval()) == entireinterval()
    @test cancelminus(interval(-1.0, Inf), emptyinterval()) == entireinterval()
    @test cancelminus(entireinterval(), emptyinterval()) == entireinterval()
    @test cancelminus(interval(-Inf, -1.0), interval(-1.0,5.0)) == entireinterval()
    @test cancelminus(interval(-1.0, Inf), interval(-1.0,5.0)) == entireinterval()
    @test cancelminus(entireinterval(), interval(-1.0,5.0)) == entireinterval()
    @test cancelminus(interval(-Inf, -1.0), entireinterval()) == entireinterval()
    @test cancelminus(interval(-1.0, Inf), entireinterval()) == entireinterval()
    @test cancelminus(emptyinterval(), interval(-Inf, -1.0)) == entireinterval()
    @test cancelminus(emptyinterval(), interval(-1.0, Inf)) == entireinterval()
    @test cancelminus(emptyinterval(), entireinterval()) == entireinterval()
    @test cancelminus(interval(-1.0,5.0), interval(-Inf, -1.0)) == entireinterval()
    @test cancelminus(interval(-1.0,5.0), interval(-1.0, Inf)) == entireinterval()
    @test cancelminus(interval(-1.0,5.0), entireinterval()) == entireinterval()
    @test cancelminus(entireinterval(), interval(-Inf, -1.0)) == entireinterval()
    @test cancelminus(entireinterval(), interval(-1.0, Inf)) == entireinterval()
    @test cancelminus(entireinterval(), entireinterval()) == entireinterval()
    @test cancelminus(interval(-5.0, -1.0), interval(-5.1,-1.0)) == entireinterval()
    @test cancelminus(interval(-5.0, -1.0), interval(-5.0,-0.9)) == entireinterval()
    @test cancelminus(interval(-5.0, -1.0), interval(-5.1,-0.9)) == entireinterval()
    @test cancelminus(interval(-10.0, 5.0), interval(-10.1, 5.0)) == entireinterval()
    @test cancelminus(interval(-10.0, 5.0), interval(-10.0, 5.1)) == entireinterval()
    @test cancelminus(interval(-10.0, 5.0), interval(-10.1, 5.1)) == entireinterval()
    @test cancelminus(interval(1.0, 5.0), interval(0.9, 5.0)) == entireinterval()
    @test cancelminus(interval(1.0, 5.0), interval(1.0, 5.1)) == entireinterval()
    @test cancelminus(interval(1.0, 5.0), interval(0.9, 5.1)) == entireinterval()
    @test cancelminus(interval(-10.0, -1.0), emptyinterval()) == entireinterval()
    @test cancelminus(interval(-10.0, 5.0), emptyinterval()) == entireinterval()
    @test cancelminus(interval(1.0, 5.0), emptyinterval()) == entireinterval()
    @test cancelminus(emptyinterval(), emptyinterval()) == emptyinterval()
    @test cancelminus(emptyinterval(), interval(-10.0, -1.0)) == emptyinterval()
    @test cancelminus(emptyinterval(), interval(-10.0, 5.0)) == emptyinterval()
    @test cancelminus(emptyinterval(), interval(1.0, 5.0)) == emptyinterval()
    @test cancelminus(interval(-5.1,-0.0), interval(-5.0, 0.0)) == interval(-0x1.999999999998P-4,0.0)
    @test cancelminus(interval(-5.1,-1.0), interval(-5.0, -1.0)) == interval(-0x1.999999999998P-4,0.0)
    @test cancelminus(interval(-5.0,-0.9), interval(-5.0, -1.0)) == interval(0.0, 0x1.9999999999998P-4)
    @test cancelminus(interval(-5.1,-0.9), interval(-5.0, -1.0)) == interval(-0x1.999999999998P-4,0x1.9999999999998P-4)
    @test cancelminus(interval(-5.0,-1.0), interval(-5.0, -1.0)) == interval(0.0,0.0)
    @test cancelminus(interval(-10.1, 5.0), interval(-10.0, 5.0)) == interval(-0x1.999999999998P-4,0.0)
    @test cancelminus(interval(-10.0, 5.1), interval(-10.0, 5.0)) == interval(0.0,0x1.999999999998P-4)
    @test cancelminus(interval(-10.1, 5.1), interval(-10.0, 5.0)) == interval(-0x1.999999999998P-4,0x1.999999999998P-4)
    @test cancelminus(interval(-10.0, 5.0), interval(-10.0, 5.0)) == interval(0.0,0.0)
    @test cancelminus(interval(0.9, 5.0), interval(1.0, 5.0)) == interval(-0x1.9999999999998P-4,0.0)
    @test cancelminus(interval(-0.0, 5.1), interval(0.0, 5.0)) == interval(0.0,0x1.999999999998P-4)
    @test cancelminus(interval(1.0, 5.1), interval(1.0, 5.0)) == interval(0.0,0x1.999999999998P-4)
    @test cancelminus(interval(0.9, 5.1), interval(1.0, 5.0)) == interval(-0x1.9999999999998P-4,0x1.999999999998P-4)
    @test cancelminus(interval(1.0, 5.0), interval(1.0, 5.0)) == interval(0.0,0.0)
    @test cancelminus(interval(-5.0, 1.0), interval(-1.0, 5.0)) == interval(-4.0,-4.0)
    @test cancelminus(interval(-5.0, 0.0), interval(-0.0, 5.0)) == interval(-5.0,-5.0)
    @test cancelminus(interval(0x1.FFFFFFFFFFFFP+0,0x1.FFFFFFFFFFFFP+0), interval(0x1.999999999999AP-4,0x1.999999999999AP-4)) == interval(0x1.E666666666656P+0,0x1.E666666666657P+0)
    @test cancelminus(interval(-0x1.999999999999AP-4,0x1.FFFFFFFFFFFFP+0), interval(-0.01,0x1.999999999999AP-4)) == interval(-0x1.70A3D70A3D70BP-4,0x1.E666666666657P+0)
    @test cancelminus(interval(0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFFp1023), interval(-0x1.FFFFFFFFFFFFFp1023,-0x1.FFFFFFFFFFFFFp1023)) == interval(0x1.FFFFFFFFFFFFFp1023,Inf)
    @test cancelminus(interval(-0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFFp1023), interval(-0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFFp1023)) == interval(0.0,0.0)
    @test cancelminus(interval(-0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFFp1023), interval(-0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFEP+1023)) == interval(0.0,0x1P+971)
    @test cancelminus(interval(-0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFFp1023), interval(-0x1.FFFFFFFFFFFFEP+1023,0x1.FFFFFFFFFFFFFp1023)) == interval(-0x1P+971,0.0)
    @test cancelminus(interval(-0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFEP+1023), interval(-0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFFp1023)) == entireinterval()
    @test cancelminus(interval(-0x1.FFFFFFFFFFFFEP+1023,0x1.FFFFFFFFFFFFFp1023), interval(-0x1.FFFFFFFFFFFFFp1023,0x1.FFFFFFFFFFFFFp1023)) == entireinterval()
    @test cancelminus(interval(0x0.0000000000001p-1022,0x0.0000000000001p-1022), interval(0x0.0000000000001p-1022,0x0.0000000000001p-1022)) == interval(0.0,0.0)
    @test cancelminus(interval(0x0.0000000000001p-1022,0x0.0000000000001p-1022), interval(-0x0.0000000000001p-1022,-0x0.0000000000001p-1022)) == interval(0x0.0000000000002p-1022,0x0.0000000000002p-1022)
    @test cancelminus(interval(0x1P-1022,0x1.0000000000002P-1022), interval(0x1P-1022,0x1.0000000000001P-1022)) == interval(0.0,0x0.0000000000001P-1022)
    @test cancelminus(interval(0x1P-1022,0x1.0000000000001P-1022), interval(0x1P-1022,0x1.0000000000002P-1022)) == entireinterval()
    @test cancelminus(interval(-0x1P+0,0x1.FFFFFFFFFFFFFP-53), interval(-0x1.FFFFFFFFFFFFEP-53,0x1P+0)) == interval(-0x1.FFFFFFFFFFFFFP-1,-0x1.FFFFFFFFFFFFEP-1)
    @test cancelminus(interval(-0x1P+0,0x1.FFFFFFFFFFFFEP-53), interval(-0x1.FFFFFFFFFFFFFP-53,0x1P+0)) == entireinterval()
end


text:
"""
    MeshLineModule

Module  for generation of meshes composed of line (curve) elements.
"""
module MeshLineModule

__precompile__(true)

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import ..FENodeSetModule: FENodeSet
import ..FESetModule: FESetL2, FESetL3
import ..MeshUtilModule: linearspace
import Statistics: mean

"""
L2block(Length::FFlt, nL::FInt)

Mesh of a 1-D block of L2 finite elements.
"""
function L2block(Length::FFlt, nL::FInt)
    fens,fes = L2blockx(collect(dropdims(linearspace(0.0,Length,nL+1)', dims=1)));
end

"""
L2blockx(xs::FFltVec)

Graded mesh of a 1-D block, L2 finite elements.
"""
function L2blockx(xs::FFltVec)
    xyz =reshape(sort(xs),length(xs),1);
    ncells=length(xs)-1;

    # create the nodes
    fens = FENodeSet(xyz);
    # Create the finite elements
    fes = FESetL2([(1:ncells) (2:ncells+1)]);

    return fens,fes;
end

"""
    L3blockx(xs::FFltVec)

Graded mesh of a 1-D block, L2 finite elements.
"""
function L3blockx(xs::FFltVec)
    fens,fes = L2blockx(xs)
    nxyz = zeros(count(fes), size(fens.xyz, 2))
    nconn = zeros(FInt, count(fes), 3)
    N = count(fens)
    for i = 1:count(fes)
        N = N + 1
        nxyz[i, :] = mean(fens.xyz[[k for k in fes.conn[i]], :], dims = 1)
        nconn[i, :] = vcat([k for k in fes.conn[i]], [N])
    end
    fens = FENodeSet([fens.xyz; nxyz]);
    # Create the finite elements
    fes = FESetL3(nconn);

    return fens,fes;
end

end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to IAU-2006 precession-nutation algorithm.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# File: ./src/transformations/precession_nutation_iau2006.jl
# ==========================================================

# Function precession_nutation_iau2006
# ------------------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example, using JD_TT = 2453101.828154745, one gets:
#
#   X = 80.531880"
#   Y =  7.273921"
#   s = -0.003027"
#
################################################################################

@testset "Function precession_nutation_iau2006" begin
    JD_TT = 2453101.828_154_745

    X,Y,s = precession_nutation_iau2006(JD_TT)

    @test X*180/pi*3600 ≈ 80.531880 atol = 5e-5
    @test Y*180/pi*3600 ≈  7.273921 atol = 5e-5
    @test s*180/pi*3600 ≈ -0.003027 atol = 1e-6
end


text:
module GoogleTrans
import HTTP
import JSON
export translate
#export batchtranslate

function package_rpc(text::AbstractString, lang_src::String, lang_tgt::String)::String
    parameter = [[text, lang_src, lang_tgt, true], [1]]
    escaped_parameter = JSON.json(parameter)
    rpc = [[["MkEWBc", escaped_parameter, nothing, "generic"]]]
    escaped_rpc = JSON.json(rpc)
    freq = "f.req=" * HTTP.escapeuri(escaped_rpc) * "&"
    return freq
end

function translate(text::AbstractString, lang_tgt="auto"::String, lang_src="auto"::String):String
    headers = Dict("Referer"=>"http://translate.google.com", "User-Agent"=>"Mozilla/5.0 (Windows NT 10.0; WOW64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/47.0.2526.106 Safari/537.36", "Content-type"=>"application/x-www-form-urlencoded;charset=utf-8")
    freq = package_rpc(text, lang_src, lang_tgt)
    r = HTTP.request("POST", "https://translate.google.com/_/TranslateWebserverUi/data/batchexecute", headers, freq)
    resp = String(r.body)
    for line in split(resp, "\n")
        if occursin("MkEWBc", line)
            response = line * "]"
            parsed = JSON.parse(response)
            #println(parsed)
            goodparsed = JSON.parse(parsed[1][3])[2][1][1][6]
            translated = ""
            for sentencelist in goodparsed
                translated = translated * " " * sentencelist[1]
            end
            return translated[2:end]
        end
    end
end

translate(inputstrings::AbstractVector{S}, lang_tgt="auto", lang_src="auto") where {S <: AbstractString} = translate.(inputstrings, lang_tgt, lang_src)

#function batchtranslate(inputstrings::Vector{AbstractString}, lang_tgt="auto"::String, lang_src="auto"::String)::Vector{String}
#    returnparts = Vector{String}()
#    for part in inputstrings
#        append!(returnparts, translate(part))
#    end
#    return returnparts
#end
end


text:
using HttpServer, Lazy

import Base.Meta.isexpr

export @app, serve

# `App` is just a box which allows the server to be
# redefined on the fly.
# In general these methods provide a simple way to
# get up and running, but aren't meant to be comprehensive.

type App
  warez
end

macro app(def)
  @assert isexpr(def, :(=))
  name, warez = def.args
  warez = isexpr(warez, :tuple) ? Expr(:call, :mux, map(esc, warez.args)...) : warez
  quote
    if isdefined($(Expr(:quote, name)))
      $(esc(name)).warez = $warez
    else
      const $(esc(name)) = App($warez)
    end
    nothing
  end
end

function http_handler(app::App)
  handler = HttpHandler((req, res) -> app.warez(req))
  handler.events["error"]  = (client, error) -> println(error)
  handler.events["listen"] = (port)          -> println("Listening on $port...")
  return handler
end

function ws_handler(app::App)
  handler = WebSocketHandler((req, client) -> app.warez((req, client)))
  return handler
end

const default_port = 8000

function serve(s::Server, port = default_port)
  @async @errs run(s, port)
  return
end

serve(h::App, port = default_port) =
  serve(Server(http_handler(h)), port)

serve(h::App, w::App, port = default_port) =
  serve(Server(http_handler(h), ws_handler(w)), port)


text:
function Beta_tr(mu, sd)
    a = ((1 - mu) / sd ^ 2 - 1 / mu) * mu ^ 2
    b = a * (1 / mu - 1)
    return a, b
end

function Gamma_tr(mu, sd)
    b = sd ^ 2 / mu
    a = mu / b
    return a, b
end

function InvGamma_tr(mu, sd)
    a = mu ^ 2 / sd ^ 2 + 2
    b = mu * (a - 1)
    return a, b
end

@model function rbc_kalman(z, m, p_f, α_prior, β_prior, ρ_prior, cache)
    α ~ truncated(Normal(α_prior[1], α_prior[2]), α_prior[3], α_prior[4])
    β_draw ~ Gamma(β_prior[1], β_prior[2])
    ρ ~ Beta(ρ_prior[1], ρ_prior[2])
    β = 1 / (β_draw / 100 + 1)
    p = [α, β, ρ]
    sol = generate_perturbation(m, p; p_f, cache)
    if !(sol.retcode == :Success)
        Turing.@addlogprob! -Inf
        return
    end
    Turing.@addlogprob! solve(sol, sol.x_ergodic, (0, length(z)); observables = z).logpdf
end

@model function rbc_second(z, m, p_f, α_prior, β_prior, ρ_prior, cache, x0 = zeros(m.n_x))
    α ~ truncated(Normal(α_prior[1], α_prior[2]), α_prior[3], α_prior[4])
    β_draw ~ Gamma(β_prior[1], β_prior[2])
    ρ ~ Beta(ρ_prior[1], ρ_prior[2])
    β = 1 / (β_draw / 100 + 1)
    p = [α, β, ρ]
    T = length(z)
    ϵ_draw ~ MvNormal(m.n_ϵ * T, 1.0)
    ϵ = map(i -> ϵ_draw[((i - 1) * m.n_ϵ + 1):(i * m.n_ϵ)], 1:T)
    sol = generate_perturbation(m, p; p_f, cache)
    if !(sol.retcode == :Success)
        Turing.@addlogprob! -Inf
        return
    end
    Turing.@addlogprob! solve(sol, x0, (0, T); noise = ϵ, observables = z).logpdf

end

@model function FVGQ20_kalman(z, m, p_f, params, cache)
    # Priors
    β_draw ~ Gamma(params.β[1], params.β[2])
    β = 1 / (β_draw / 100 + 1)
    h ~ Beta(params.h[1], params.h[2])
    ϑ = 1.0
    κ ~ truncated(Normal(params.κ[1], params.κ[2]), params.κ[3], params.κ[4])
    α ~ truncated(Normal(params.α[1], params.α[2]), params.α[3], params.α[4])
    θp ~ Beta(params.θp[1], params.θp[2])
    χ ~ Beta(params.χ[1], params.χ[2])
    γR ~ Beta(params.γR[1], params.γR[2])
    γy ~ truncated(Normal(params.γy[1], params.γy[2]), params.γy[3], params.γy[4])
    γΠ ~ truncated(Normal(params.γΠ[1], params.γΠ[2]), params.γΠ[3], params.γΠ[4])
    Πbar_draw ~ Gamma(params.Πbar[1], params.Πbar[2])
    Πbar = Πbar_draw / 100 + 1
    ρd ~ Beta(params.ρd[1], params.ρd[2])
    ρφ ~ Beta(params.ρφ[1], params.ρφ[2])
    ρg ~ Beta(params.ρg[1], params.ρg[2])
    g_bar ~ Beta(params.g_bar[1], params.g_bar[2])
    σ_A ~ InverseGamma(params.σ_A[1], params.σ_A[2])
    σ_d ~ InverseGamma(params.σ_d[1], params.σ_d[2])
    σ_φ ~ InverseGamma(params.σ_φ[1], params.σ_φ[2])
    σ_μ ~ InverseGamma(params.σ_μ[1], params.σ_μ[2])
    σ_m ~ InverseGamma(params.σ_m[1], params.σ_m[2])
    σ_g ~ InverseGamma(params.σ_g[1], params.σ_g[2])
    Λμ ~ Gamma(params.Λμ[1], params.Λμ[2])
    ΛA ~ Gamma(params.ΛA[1], params.ΛA[2])
    # Likelihood
    θ = [β, h, ϑ, κ, α, θp, χ, γR, γy, γΠ, Πbar, ρd, ρφ, ρg, g_bar, σ_A, σ_d, σ_φ, σ_μ, σ_m, σ_g, Λμ, ΛA]
    println(θ)
    sol = generate_perturbation(m, θ; p_f, cache)
    if !(sol.retcode == :Success)
        Turing.@addlogprob! -Inf
    else
        z_trend = params.Hx * sol.x + params.Hy * sol.y
        z_detrended = map(i -> z[i] - z_trend, eachindex(z))
        Turing.@addlogprob! solve(sol, sol.x_ergodic, (0, length(z_detrended)); observables = z_detrended).logpdf
    end
end

# Joint likelihood
@model function FVGQ20_joint(z, m, p_f, params, cache, x0 = zeros(m.n_x))
    T = length(z)
    # Priors
    β_draw ~ Gamma(params.β[1], params.β[2])
    β = 1 / (β_draw / 100 + 1)
    h ~ Beta(params.h[1], params.h[2])
    ϑ = 1.0
    κ ~ truncated(Normal(params.κ[1], params.κ[2]), params.κ[3], params.κ[4])
    α ~ truncated(Normal(params.α[1], params.α[2]), params.α[3], params.α[4])
    θp ~ Beta(params.θp[1], params.θp[2])
    χ ~ Beta(params.χ[1], params.χ[2])
    γR ~ Beta(params.γR[1], params.γR[2])
    γy ~ truncated(Normal(params.γy[1], params.γy[2]), params.γy[3], params.γy[4])
    γΠ ~ truncated(Normal(params.γΠ[1], params.γΠ[2]), params.γΠ[3], params.γΠ[4])
    Πbar_draw ~ Gamma(params.Πbar[1], params.Πbar[2])
    Πbar = Πbar_draw / 100 + 1
    ρd ~ Beta(params.ρd[1], params.ρd[2])
    ρφ ~ Beta(params.ρφ[1], params.ρφ[2])
    ρg ~ Beta(params.ρg[1], params.ρg[2])
    g_bar ~ Beta(params.g_bar[1], params.g_bar[2])
    σ_A ~ InverseGamma(params.σ_A[1], params.σ_A[2])
    σ_d ~ InverseGamma(params.σ_d[1], params.σ_d[2])
    σ_φ ~ InverseGamma(params.σ_φ[1], params.σ_φ[2])
    σ_μ ~ InverseGamma(params.σ_μ[1], params.σ_μ[2])
    σ_m ~ InverseGamma(params.σ_m[1], params.σ_m[2])
    σ_g ~ InverseGamma(params.σ_g[1], params.σ_g[2])
    Λμ ~ Gamma(params.Λμ[1], params.Λμ[2])
    ΛA ~ Gamma(params.ΛA[1], params.ΛA[2])
    ϵ_draw ~ MvNormal(m.n_ϵ * T, 1.0)
    ϵ = map(i -> ϵ_draw[((i-1)*m.n_ϵ+1):(i*m.n_ϵ)], 1:T)
    # Likelihood
    θ = [
        β,
        h,
        ϑ,
        κ,
        α,
        θp,
        χ,
        γR,
        γy,
        γΠ,
        Πbar,
        ρd,
        ρφ,
        ρg,
        g_bar,
        σ_A,
        σ_d,
        σ_φ,
        σ_μ,
        σ_m,
        σ_g,
        Λμ,
        ΛA,
    ]
    # println(θ)
    sol = generate_perturbation(m, θ; p_f, cache)
    if !(sol.retcode == :Success)
        Turing.@addlogprob! -Inf
        return
    end
    z_trend = params.Hx * sol.x + params.Hy * sol.y
    z_detrended = map(i -> z[i] - z_trend, eachindex(z))
    Turing.@addlogprob! solve(sol, x0, (0, T); noise = ϵ, observables = z_detrended).logpdf
end


text:
struct Timing
  ns::UInt
end

function Base.:-(a::Timing,b::Timing)
  ns = a.ns-b.ns
  Timing(ns)
end

@inline Timing() = Timing(time_ns())

mutable struct PTimer{A,B<:Dict}
  parts::A
  timings::B
  current::Timing
  verbose::Bool
end

function PTimer(parts::AbstractPData{<:Integer};verbose::Bool=false)
  current = Timing()
  timing = map_parts(p->current,parts)
  T = typeof(timing)
  timings = Dict{String,T}()
  PTimer(parts,timings,current,verbose)
end

function Base.getproperty(t::PTimer, sym::Symbol)
  if sym == :data
    data = map_main(t.parts) do part
      d = (min=0.0,max=0.0,avg=0.0)
      T = typeof(d)
      Dict{String}{T}()
    end
    for (name,timing) in t.timings
      timing_main = gather(timing)
      map_main(timing_main,data) do part_to_timing,data
        ns = map(i->Float64(i.ns)/1.0e9,part_to_timing)
        d = (min=minimum(ns),max=maximum(ns),avg=sum(ns)/length(ns))
        data[name] = d
      end
    end
    return data
  else
    getfield(t, sym)
  end
end

function Base.propertynames(x::PTimer, private=false)
  (fieldnames(typeof(x))...,:data)
end

function tic!(t::PTimer)
  t.current = Timing()
end

function toc!(t::PTimer,name::String)
  current = Timing()
  dt = current-t.current
  timing = map_parts(p->dt,t.parts)
  t.timings[name] = timing
  if t.verbose == true
    map_main(timing) do i
      println("[$(lstrip(_nice_time(Float64(i.ns)/1.0e9))) s in MAIN] $name")
    end
  end
  t.current = Timing()
end

function Base.show(io::IO,k::MIME"text/plain",t::PTimer)
  print_timer(io,t)
end

function print_timer(
  io::IO,
  t::PTimer;
  linechars::Symbol=:unicode)

  map_main(t.data) do data
    kv = collect(data)
    v = map(i->i[2].max,kv)
    i = reverse(sortperm(v))
    sorteddata = kv[i]
    _print_on_main(io,sorteddata,linechars)
  end
end

print_timer(t::PTimer; kwargs...) = print_timer(stdout,t;kwargs...)

function _print_on_main(io,data,linechars)
  w = length(_nice_time(eps()))
  longest_name = maximum(map(i->length(i[1]),data))
  longest_name = max(longest_name,length("Section"))
  _print_header(io,longest_name,w,linechars)
  for (name,d) in data
    _print_section(io,longest_name,name,d)
  end
  _print_footer(io,longest_name,w,linechars)
end

_nice_time(t) = @sprintf("%12.3e",t)

function _print_header(io,longest_name,w,linechars)
  rule = linechars == :unicode ? "─" : "-"
  sec_header = rpad("Section",longest_name)
  max_header = lpad("max",w)
  min_header = lpad("min",w)
  avg_header = lpad("avg",w)
  header = sec_header*max_header*min_header*avg_header
  header_w = length(header)
  println(io,rule^header_w)
  println(io,header)
  println(io,rule^header_w)
end

function _print_section(io,longest_name,name,d)
  sec_header = rpad(name,longest_name)
  max_header = _nice_time(d.max)
  min_header = _nice_time(d.min)
  avg_header = _nice_time(d.avg)
  header = sec_header*max_header*min_header*avg_header
  println(io,header)
end

function _print_footer(io,longest_name,w,linechars)
  rule = linechars == :unicode ? "─" : "-"
  header_w = longest_name+3*w
  println(io,rule^header_w)
end


text:
# from DifferentialEquations.jl
macro CSI_str(str)
    return :(string("\x1b[", $(esc(str)), "m"))
end
const TYPE_COLOR = CSI"36"
const NO_COLOR = CSI"0"

inbounds(x::T1, xmin::T2, xmax::T3) where {T1<:Real, T2<:Real, T3<:Real} =
	max(xmin, min(x, xmax))

getiterator(x::NTuple{N,AbstractVector}) where N = Iterators.product(x...)
getiterator(x::NTuple{1,AbstractVector}) = x[1]

const FuncOrNothing = Union{Function,Nothing}

"""Returns a modified version of the function f that if forced to stay between
the bounds defined by the tuple of vectors tVectors"""
function wrapinbounds(f, tVectors::NTuple{N,Vector{C}}) where {N,C}
	lowerbounds::NTuple{N,C} = minimum.(tVectors)
	upperbounds::NTuple{N,C} = maximum.(tVectors)
	return function f2(args...)
		inbounds.(f(args...), lowerbounds, upperbounds)
	end
end
function wrapinbounds(f, tVectors::NTuple{1,Vector{C}}) where {C}
	lowerbounds::C = minimum(tVectors[1])
	upperbounds::C = maximum(tVectors[1])
	return function f1(args...)
		inbounds(f(args...), lowerbounds, upperbounds)
	end
end


text:
defaultOutputFormat(plt::AbstractPlot) = "png"

function png(plt::AbstractPlot, fn::@compat(AbstractString))
  fn = addExtension(fn, "png")
  io = open(fn, "w")
  writemime(io, MIME("image/png"), plt)
  close(io)
end
png(fn::@compat(AbstractString)) = png(current(), fn)

function svg(plt::AbstractPlot, fn::@compat(AbstractString))
  fn = addExtension(fn, "svg")
  io = open(fn, "w")
  writemime(io, MIME("image/svg+xml"), plt)
  close(io)
end
svg(fn::@compat(AbstractString)) = svg(current(), fn)

function pdf(plt::AbstractPlot, fn::@compat(AbstractString))
  fn = addExtension(fn, "pdf")
  io = open(fn, "w")
  writemime(io, MIME("application/pdf"), plt)
  close(io)
end
pdf(fn::@compat(AbstractString)) = pdf(current(), fn)

function ps(plt::AbstractPlot, fn::@compat(AbstractString))
  fn = addExtension(fn, "ps")
  io = open(fn, "w")
  writemime(io, MIME("application/postscript"), plt)
  close(io)
end
ps(fn::@compat(AbstractString)) = ps(current(), fn)

function tex(plt::AbstractPlot, fn::@compat(AbstractString))
  fn = addExtension(fn, "tex")
  io = open(fn, "w")
  writemime(io, MIME("application/x-tex"), plt)
  close(io)
end
tex(fn::@compat(AbstractString)) = tex(current(), fn)

# ----------------------------------------------------------------

@compat const _savemap = Dict(
    "png" => png,
    "svg" => svg,
    "pdf" => pdf,
    "ps"  => ps,
    "tex" => tex,
  )

function getExtension(fn::@compat(AbstractString))
  pieces = split(fn, ".")
  length(pieces) > 1 || error("Can't extract file extension: ", fn)
  ext = pieces[end]
  haskey(_savemap, ext) || error("Invalid file extension: ", fn)
  ext
end

function addExtension(fn::@compat(AbstractString), ext::@compat(AbstractString))
  try
    oldext = getExtension(fn)
    if oldext == ext
      return fn
    else
      return "$fn.$ext"
    end
  catch
    return "$fn.$ext"
  end
end

function savefig(plt::AbstractPlot, fn::@compat(AbstractString))

  # get the extension
  local ext
  try
    ext = getExtension(fn)
  catch
    # if we couldn't extract the extension, add the default
    ext = defaultOutputFormat(plt)
    fn = addExtension(fn, ext)
  end

  # save it
  func = get(_savemap, ext) do
    error("Unsupported extension $ext with filename ", fn)
  end
  func(plt, fn)
end
savefig(fn::@compat(AbstractString)) = savefig(current(), fn)

# ---------------------------------------------------------

gui(plt::AbstractPlot = current()) = display(PlotsDisplay(), plt)

# override the REPL display to open a gui window
Base.display(::Base.REPL.REPLDisplay, ::MIME"text/plain", plt::AbstractPlot) = gui(plt)

# a backup for html... passes to svg
function Base.writemime(io::IO, ::MIME"text/html", plt::AbstractPlot)
    writemime(io, MIME("image/svg+xml"), plt)
end

# ---------------------------------------------------------
# IJulia
# ---------------------------------------------------------

const _ijulia_output = ASCIIString["text/html"]

function setup_ijulia()
    # override IJulia inline display
    if isijulia()
        @eval begin
            import IJulia
            export set_ijulia_output
            function set_ijulia_output(mimestr::ASCIIString)
                # info("Setting IJulia output format to $mimestr")
                global _ijulia_output
                _ijulia_output[1] = mimestr
            end
            function IJulia.display_dict(plt::AbstractPlot)
                global _ijulia_output
                Dict{ASCIIString, ByteString}(_ijulia_output[1] => sprint(writemime, _ijulia_output[1], plt))
            end
        end

        # IJulia.display_dict(plt::AbstractPlot) = Dict{ASCIIString, ByteString}("text/html" => sprint(writemime, "text/html", plt))
        set_ijulia_output("text/html")
    end
end

# ---------------------------------------------------------
# Atom PlotPane
# ---------------------------------------------------------

function setup_atom()
    # @require Atom begin
    if isatom()
        # @eval import Atom, Media
        @eval import Atom

        # connects the render function
        for T in (GadflyBackend,ImmerseBackend,PyPlotBackend,GRBackend)
            Atom.Media.media(AbstractPlot{T}, Atom.Media.Plot)
        end
        # Atom.Media.media{T <: Union{GadflyBackend,ImmerseBackend,PyPlotBackend,GRBackend}}(Plot{T}, Atom.Media.Plot)

        # Atom.displaysize(::AbstractPlot) = (535, 379)
        # Atom.displaytitle(plt::AbstractPlot) = "Plots.jl (backend: $(backend(plt)))"

        # this is like "display"... sends an html div with the plot to the PlotPane
        function Atom.Media.render(pane::Atom.PlotPane, plt::AbstractPlot)
            Atom.Media.render(pane, Atom.div(Atom.d(), Atom.HTML(stringmime(MIME("text/html"), plt))))
        end

        # function Atom.Media.render(pane::Atom.PlotPane, plt::Plot{PlotlyBackend})
        #     html = Media.render(pane, Atom.div(Atom.d(), Atom.HTML(stringmime(MIME("text/html"), plt))))
        # end
    end
end


text:
using StatsBase:sample

"""
    struct TabularQ <: AbstractQApproximator{Int, Int}
        table::Array{Float64, 2}
    end

Using a `table` of type `Array{Float64,2}` to record the action value of each state.
"""
struct TabularQ <: AbstractQApproximator{Int, Int}
    table::Array{Float64, 2}
end

"""
    TabularQ(ns::Int, na::Int=1, init::Float64=0.)

Initial a table of size `(ns, na)` filled with value of `init`.
"""
TabularQ(ns::Int, na::Int=1, init::Float64=0.) = TabularQ(fill(init, ns, na))

(Q::TabularQ)(s::Int, ::Val{:dist}) = Q.table[s, :]
(Q::TabularQ)(s::Int, ::Val{:max}) = @views findmax(Q.table[s, :])[1]
(Q::TabularQ)(s::Int, ::Val{:argmax}) = @views sample(findallmax(Q.table[s, :])[2])
(Q::TabularQ)(s::Int, a::Int) = Q.table[s, a]
(Q::TabularQ)(s::Int) = Q(s, Val(:dist))

"""
!!! warning
    Without bound check! Unsafe!
"""
function update!(Q::TabularQ, s::Int, errors::Vector{Float64})
    for (i, e) in enumerate(errors)
        Q.table[s, i] += e
    end
end

update!(Q::TabularQ, s::Int, a::Int, e::Float64) = Q.table[s, a] += e